Survey on Adsorption of Low Molecular Weight Compounds in Cu-BTC Metal–Organic Framework: Experimental Results and Thermodynamic Modeling

This contribution aims at providing a critical overview of experimental results for the sorption of low molecular weight compounds in the Cu-BTC Metal–Organic Framework (MOF) and of their interpretation using available and new, specifically developed, theoretical approaches. First, a literature review of experimental results for the sorption of gases and vapors is presented, with particular focus on the results obtained from vibrational spectroscopy techniques. Then, an overview of theoretical models available in the literature is presented starting from semiempirical theoretical approaches suitable to interpret the adsorption thermodynamics of gases and vapors in Cu-BTC. A more detailed description is provided of a recently proposed Lattice Fluid approach, the Rigid Adsorbent Lattice Fluid (RALF) model. In addition, to deal with the cases where specific self- and cross-interactions (e.g., H-bonding, Lewis acid/Lewis base interactions) play a role, a modification of the RALF model, i.e., the RALFHB model, is introduced here for the first time. An extension of both RALF and RALFHB is also presented to cope with the cases in which the heterogeneity of the rigid adsorbent displaying a different kind of adsorbent cages is of relevance, as it occurs for the adsorption of some low molecular weight substances in Cu-BTC MOF.


Introduction
Metal-organic frameworks (MOFs) are a class of hybrid crystalline materials composed of transition metal divalent (Zn 2+ , Cu 2+ , Mg 2+ , etc.) or trivalent (Cr 2+ , Al 2+ , Fe 2+ , etc.) cations joined together by organic ligands (phosphonates, carboxylates, or sulfonates) via strong coordinative bonds. Due to their structures, MOFs are characterized, as compared to activated carbon and zeolites, by significantly higher surface areas and a perfectly ordered molecular arrangement with pore sizes between 3 and 35 Å [1]. Different distributions of pore size and shape as well as functionalities can be achieved by simply selecting the metal center and/or the ligand [1]. This unique feature allows these materials to be used in numerous applications [2], ranging from storage media [3] and adsorbents for separation processes [4,5] to drug delivery carriers and catalysts [6][7][8][9]. In this perspective, H 2 , CH 4 , and CO 2 are gases of particular interest due to their environmental and economic importance. H 2 is a promising energy carrier for the substitution of liquid fuel resources applied in the automobile sector, being environmentally friendly and a fully renewable energy source [10]. CO 2 is emitted from the combustion of fossil fuels and its emissions give a major contribution to the greenhouse effect and thus to climate change [11]. In  [30]. Copyright 2013 American Chemical Society.
In the present survey, the attention is focused on the adsorption of low molecular weight compounds in Cu-BTC. The objectives are to critically review available experimental results from the literature with a special consideration for approaches exploiting the wealth of information obtainable from vibrational spectroscopy and to extensively discuss the theoretical models developed in the literature to interpret sorption data in MOFs. In addition, a modified version of an available model based on a lattice fluid approach is introduced and tested against experimental sorption isotherms of several low molecular weight compounds in Cu-BTC. Table 1 reports a compendium of the most significant experimental works available in the literature on this topic. It can be observed that, in the past two decades, various reports have been focused on the adsorption capacities of small apolar gases in Cu-BTC, such as noble gases, N 2 , H 2 , CH 4 , and CO 2 . Only a few reports regarding the adsorption of small polar alcohols have been published. This is particularly notable for methanol since it is one of the most abundant organic compounds in the atmosphere and, therefore, its removal in the atmosphere is of great significance [37].  [72] a adsorption temperature. b adsorption pressure range. c values specified as P/P 0 .

Vibrational Spectroscopy and Other Spectroscopic Techniques
One of the main tools to investigate MOFs at the molecular level is vibrational spectroscopy (e.g., FTIR, Raman, ultrafast 2D-IR) [73]. Regarding MOFs, several vibration frequencies can be assigned to the different constituents of the framework or the guest molecules that interact with the sites belonging to the different cages. Understanding the nature of these interactions is crucial to tailor and improve MOFs properties. Structural changes of the host and/or guest molecules (e.g., adsorbed gas molecules or catalytic reactions) can be attributed and assigned to specific signal intensity changes and/or observed frequency shifts [74]. These occurrences indicate possible interactions or changes within the MOF framework. However, a correct and accurate assignment of the vibrational modes represents a challenge in vibrational spectroscopy [74]. In fact, several experimental studies making use of vibrational spectroscopy are often accompanied by theoretical interpretations, employing first-principles simulations [74][75][76] since combining theory with experiments provides complementary information for a deeper molecular characterization. IR spectroscopy is processed through Fourier transform (FTIR) and, as its main advantage, allows for a combined solid-state technique and sampling flexibility. For instance, collection modes in FTIR such as Fourier Transform Attenuated Total Reflection (FT-ATR) and Diffuse Reflection Infrared Fourier Transform (DRIFT) do not require any sample preparation. Transmission FTIR spectra can be collected on a range of solid-state samples such as pressed pellets and thin films, while DRIFT only allows for the measurement of fine powders.
These techniques are very sensitive to intermolecular interactions, such as hydrogenbonding (H-bonding) between water and host materials, allowing for the identification and quantification of the different molecular aggregates that are being formed [77,78]. The formation of strong H-bonds can lead to extended networks that severely affect fundamental processes such as hydration processes [79], chemical synthesis and reactions [80,81], heat dissipation [82,83], and macroscopic structure formations [83,84]. To this regard, it is essential to provide an accurate description regarding the formation of the water H-bonding network in MOFs, that can unambiguously discriminate between the framework-water and water-water interactions. For instance, it is noteworthy to mention the study conducted by Singha et al. [85] in which the authors highlighted the complex mechanism that regulates the water adsorption on MIL-53(Al) by means of DRIFT. Characteristic peaks of the OHstretching vibration were related to water interactions with the adsorbent sites: at lower hydration levels the water interacted tightly with multiple sites, while at higher hydration levels the water interacted with fewer sites. Another study conducted on a Co2Cl2BTDD by DRIFT [78] elucidated the adsorption mechanism as a function of relative humidity (RH). It was found that water strongly binds with the open Co 2+ sites of the framework to subsequently form one-dimensional chains of H-bonded molecules that develop bridges between the Co 2+ sites. Upon an increase in RH, the water chains filled the pores of the network and occupied the entire pore volume before RH attained 30%.
In the context of Cu-BTC MOF, different works have been performed to elucidate the interactions of guest species with the framework. It is noteworthy to note the spectroscopic study regarding the adsorbate-adsorbent interaction in HKUST-1 performed by Bordiga et al. [23,40]. In their first work regarding this topic [40], the authors investigated the dehydration process (activation) by means of the XRD, UV−Vis, EXAFS, XANES, and Raman spectroscopies. They experimentally showed that the removal of coordinated water molecules, chemically bound to the Cu 2+ sites, led to an unchanged oxidation state of copper, a preserved crystalline nature of the material, and promoted the reduction in the cell volume due to the shrinking of the [Cu 2 C 4 O 8 ] cage. In the dehydrated state, they observed the formation of labile Cu 2+ ···CO and Cu 2+ ···H 2 adducts, detecting for the first time the signal of Cu(II) carbonyl and dihydrogen complexes formed inside a crystalline microporous hosting matrix. In a follow-up study, Bordiga et al. [23] improved the preparation method for the Cu-BTC synthesis and performed IR spectroscopic measurements in transmission mode using a properly designed cryogenic cell, assessing the interaction of HKUST-1 sites with several adsorbates, such as NO, CO 2 , CO, N 2 , and H 2 . Interestingly, the interaction of CO 2 , CO, and N 2 allowed to distinguish between a first type of Cu 2+ sites located at the external faces of the crystals and a second type of Cu 2+ sites regularly contained within the cages of the framework.
IR and Raman spectroscopy are complementary techniques with different selection rules and are often implemented jointly to study the adsorption of gases into MOFs. Raman spectroscopy is based on the detection of photons that are inelastically scattered from the sample under observation when it interacts with the radiation of a single frequency laser [85]. As the laser frequency is shifted up or down due to the interaction of the molecular vibrations, the produced spectral lines (fingerprint) correspond to the different vibrational modes of the sample material. Raman spectroscopy presents the main advantage, as compared to FTIR spectroscopy, of being able to collect high-resolution vibrational spectra at very low wavenumbers (ca. 10 cm −1 ). To this regard, metal containing materials are characterized by low frequency modes, such as metal-ligand stretches [85]. In general, Raman spectroscopy methods allow for the monitoring of solid-state samples such as pressed powders, thin films, and suspensions. However, the major limit for these techniques when studying the interactions between MOFs and adsorbates is the sample fluorescence as even weak fluorescent backgrounds can overcome Raman signals. In fact, several MOF samples that contain organic building blocks can lead to fluorescence phenomena that can mask the Raman signal [85]. In this sense, Resonance Raman (RR) may not be suitable since it relies on sources at frequencies near to those of a molecule's electronic transition [86]. For this reason, Raman spectroscopy is suitable for MOF materials with no characteristic fluorescence [85].
Raman spectroscopy has been widely used to investigate the hydrolyzation and water stability of MOF networks [40,87,88]. Notably, for Cu-BTC frameworks, in [88], the authors reported an experimental investigation by means of Raman spectroscopy on the decomposition process of Cu-BTC exposed to air moisture at 300 K and 70% of RH. Raman measurements indicated structural deterioration of the framework due to hydrolysis which affected a significant fraction of the Cu-O bonds of the crystal. This occurrence led to an irreversible process for exposure times longer than 20 days. Raman spectra revealed a shifting of the peak positions and variable intensities of the main Raman bands of the material, attributed to Cu-Cu, Cu-O, and O-C-O stretching modes. Interestingly, the coexistence of two types of paddle-wheels with different structures was observed, corresponding to hydrolyzed and nonhydrolyzed paddle-wheels, as confirmed by the detection of the Raman peak split attributed to the Cu-Cu vibration in two well distinguishable components.
Other significant techniques suitable to monitor and characterize the MOF structure are optical electronic spectroscopy and X-ray spectroscopy. Optical electronic spectroscopy allows for the detection of electronic energy levels and bonding features of different molecules and materials [89,90]. Structural frameworks formed by transition metal-based complexes and coordination polymers, such as MOFs, can exhibit a large range of electronic behaviors (from semiconductor to conductor depending on the framework structure) and, therefore, their electron transitions fall in the UV-visible and Near Infrared (NIR) wavelength range [89,90]. Therefore, optical spectroscopy methods, such as absorption and emission spectroscopy, are suitable to investigate the electronic structure of MOFs. These techniques allow for the treating of homogenous microcrystalline powders and are easily tunable for in situ measurements at controlled environmental conditions. X-ray spectroscopy allows for the analysis of electron transitions upon the absorption or emission of X-ray photons [91,92]. Based on an excited state induced by the energy of a photon, an atom moves to a higher energy level and then returns to its "unexcited" energy level. These energy transitions translate to the emission of photons with a wavelength characteristic of the sample material under observation [91,92]. In recent years, regarding MOFs, the use of techniques based on the X-ray absorption spectroscopy (XAS) has spread. That is, X-ray absorption near edge structure (XANES) and extended X-ray absorption fine structure (EXAFS) allow for the analysis of atomic distances, the coordination geometry, and oxidation state of a specific metal element, making them suitable to collect data regarding the MOF structural changes and the host-guest interactions. Similarly, in this case, these techniques deal with solid-state samples and allow in situ monitoring of the adsorption processes.
Reference to several significant works based on vibrational spectroscopy conducted on Cu-BTC networks and aimed at elucidating their interactions with adsorbates is provided in Table 2.

Theoretical Modeling and Simulation Approaches for Sorption of Low Molecular Weight Compounds in MOFs
Vibrational spectroscopy is often assisted by theoretically calculated frequencies that can help in assigning the spectroscopic bands of interest. One simple method consists of the normal coordinate analysis that is usually used to provide a harmonic vibrational frequency [95][96][97][98]. However, this method has often been inconsistent with the assignment of the spectral bands of complex materials such as MOFs [98]. In this respect, a first alternative, preferred by many researchers, consists of the empirical modeling that provides a quick and effortless evaluation of the adsorption phenomena.
The Langmuir-Freundlich (LF) isotherm, widely known also as Sips isotherms [99], is a semiempirical three-parameter model that includes mathematical features of both the Langmuir and Freundlich isotherms. Although the thermodynamic consistency of this model exhibits flaws in the regions of very low pressure (since it does not recover the Henry's law limit), the simple form of the equation allows for the modeling of either subcritical or supercritical isotherms, without requiring for the definition of the saturation pressure for the adsorbate. The equation is the following [50]: where q represents the adsorbed amount, p is the pressure, q max is the maximum adsorption capacity, b is the affinity constant, and n is the heterogeneity coefficient, respectively. In particular, for n = 1, the Sips isotherm reduces to the classic Langmuir isotherm, applicable to homogeneous adsorbent-adsorbate systems. Sips parameters are dependent on the temperature [100] but, typically, q max and n are considered independent from the temperature to keep the model application simpler. In the work performed by Aprea et al. [50], the Sips model was adopted to model the CO 2 adsorption isotherms on a laboratory-synthesized Cu-BTC framework at several temperatures (283, 293, 318, and 343 K) and for pressures up to 0.1 MPa. Although it was possible to gather some preliminary data regarding the Cu-BTC saturation capacity and a homogenous-type (Langmuir) adsorption system, this simple model did not allow for the retrieval of information regarding the nature of the adsorbate-adsorbent interaction. Under this perspective, the Virial-Langmuir (VL) model allows for the gathering of some information regarding the nature of the adsorbate-adsorbent interaction. This model consists of the addition of two virial coefficients to the Henry constant and naturally recovers Henry's law in the low concentration limit [101]. The equation is expressed as follows [45,54]: where β is the Henry constant, and b and c are the second and third virial coefficients, respectively. For instance, in [54], authors conducted a comparative adsorption study of three gases, such as CO, CO 2 , and CH 4 , on two adsorbents, namely Cu-BTC (or HKUST-1) and Cr-BDC (or MIL-101). The gravimetric adsorption equilibrium measurements on the samples were performed using a magnetic suspension balance at three different temperatures (295, 318, and 353 K) and pressures ranging from 0 to 10 MPa. In this instance, the use of the model allowed the authors to firstly evaluate the enthalpy of adsorption at zero loading and the enthalpy variation with loading at 295 K. In the former case, the enthalpy of adsorption of the three gases only resulted in small differences, due to coordinatively unsaturated metal centers present in the Cu-BTC framework [49] that were either not open or hindered by the presence of solvent molecules left over from the synthesis procedure; in the latter case, both CH 4 and CO 2 showed a slight variation in enthalpy with loading, while CO only showed a considerable decrease in enthalpy of adsorption. This occurrence was attributed to electrostatic interactions that dominate only the low loading region, while as the sites available for the interaction are progressively filled up, the enthalpy of adsorption drops down sharply. From the adsorption isotherm data, it was also possible to observe that CO 2 exhibited the highest capacity, while CH 4 had a lower capacity than CO in the low-pressure region and then progressively exceeded it in the high-pressure region. The initial behavior of the comparison between CH 4 and CO was also confirmed by Henry's constant (evaluated from the VL model), that is higher for CO than for CH 4 at 295 K due to its dipole moment. However, as interaction sites are progressively filled up, initial electrostatic interactions are overcome by dispersion interactions and the larger polarizability of CH 4 results in its higher capacity. Empirical modeling has provided a theoretical framework not only for the guest-host interaction, but for studying multicomponent adsorption mechanisms as well. To this aim, the Dubinin-Astakhov (DA) model is widely used for modeling the single-component isotherms and provides a fair prediction of the multicomponent isotherms [102]. As opposed to the Langmuir adsorption isotherm, this model is related to the micropore volume filling without the formation of successive surface layers [102]. The equation is expressed as follows [102]: where R is the ideal gas constant, T is the temperature, ε is the characteristic energy of adsorption, p s is the saturation pressure, and n is the heterogeneity parameter, respectively. The introduction of the parameter n allows us to account for the surface heterogeneity typical of microporous adsorbents such as MOFs. In their work [63], Gomez et al. used the DA model for a non-linear regression of single-component adsorption isotherms measured experimentally at 298 K and in a pressure range of 0-5 MPa. Based on that, the authors were then able to predict the adsorption isotherms of binary and ternary mixtures containing CO 2 , CH 4 , and N 2 within Cu-BTC frameworks. Although empirical models are readily implemented to retrieve data on the adsorbentadsorbate interaction, such as the isosteric heat of adsorption, results are not always consistent with experimental observations [35,43,49,51,57,65,69,70,93,103]. To provide a robust background on the mechanisms taking place during the gas adsorption on MOF materials and an interpretation of their relative vibrational spectra, over the last decade different ab initio computational methods have been implemented. For Cu-BTC frameworks, large efforts have been made in this direction using both density functional theory (DFT) [49,57,65,103] and Grand Canonical Monte Carlo (GCMC) calcula-tions [35,43,51,69,70,93]. Typically, a large computational effort is required for these calculations as the structure of MOFs comprise large crystal cells. Therefore, the use of these approaches is generally limited to small areas, such as clusters, or is combined with the use of molecular dynamics (MD) as well.
DFT methods are usually preferred as they possess a relatively higher computational efficiency and present an acceptable accuracy in predicting the vibrational spectra with their relative band intensities and frequencies [96]. A major drawback of DFT methods is the incapability of calculating weak van der Waals forces and, hence, methods such as dispersion-corrected DFT (DFT-D) and van der Waals-DFT are used to account for this discrepancy [95]. Conversely, MD simulations of MOFs rely on the development of force fields (FFs) that can describe the molecular interaction between the framework and the guest molecules. Commonly used Force Fields (FFs) for MOFs are the general Amber force field (GAFF) [104] and the universal force field (UFF) [105] that allow us to accurately model the organic links, but are less effective in describing the coordination and the geometry of the metal center [106]. To address this issue, FF extension for MOFs has also been developed, i.e., the MOF-FF [107] and the Quick-FF [108].
Concerning the Cu-BTC frameworks, some works warrant mention. For instance, Supronowicz et al. [74] studied the interactions of CO, CO 2 , OCS, SO 2 , NO, NO 2 , N 2 O, NH 3 , PH 3 , and other small molecules with the undercoordinated metal centers of the HKUST-1, by means of the DFT. Authors retrieved the adsorption energies on the Cu 2+ sites of the paddle-wheel and found the following ranking: They classified the observed interactions into three categories: (1) weak physisorption, (2) polarization and electrostatics, and (3) strong acid-base.
García-Pérez et al. [47] investigated the adsorption of several quadrupolar and nonpolar gases on Cu-BTC by means of combined adsorption isotherms and Monte Carlo (MC) simulations. In this study, the authors could identify four main adsorption sites: site I close to the copper atoms, site I within the larger cavities, site II located in the center of the smaller octahedral cages, and site III at the windows of the four open faces of the octahedral cage. Monte Carlo simulations allowed us to detect the octahedral cages (sites II and III) and the big cages (site I ) as the preferred positions for adsorption, while in the case of site I (near the copper atoms), sites remain empty over the entire range of analyzed pressures, possibly due to the reduced accessibility of these sites. Interestingly, the occupation of the different sites by CH 4 and C 2 H 6 exhibited small differences as compared to O 2 and N 2 ; this finding being attributed to the quadrupole moment of the polar molecules. While CH 4 molecules predominantly occupied the sites of type II, the N 2 occupied both I and II an equivalent amount. The molecular sitting for O 2 showed an intermediate behavior between those observed for CH 4 and N 2 .
In their work, Van Assche et al. [58] studied the adsorption of various polar adsorbates (such as methanol, ethanol, 1-propanol, 2-propanol, 1-butanol, 1-hexanol, water, acetone, acetonitrile, tetrahydrofuran, and N,N-dimethylformamide) as well as apolar adsorbates (such as n-hexane, n-heptane, and m-xylene cyclohexane) on Cu-BTC frameworks. The authors observed that alcohols characterized by a longer carbon chain (and thus less polar) have higher uptakes at lower vapor pressures. Interestingly, a significant two-step uptake was noticed for smaller alcohols in the measured vapor pressure range, the effect being more remarkable for methanol and ethanol. Regarding the alkanes, these molecules filled up the sites of the Cu-BTC structure at low vapor pressures due to their favorable interaction with the host structure. Despite the strong interaction of polar adsorbates with the Cu-BTC structure, the material also showed a quite apolar nature due to the presence of the aromatic counterpart. This behavior was also observed when performing GCMC simulations at 313-343 K for methanol using the Cu-BTC crystal structure proposed by Chui et al. [22].
Listed in Table 2 is a series of combined experimental and theoretical approaches used to characterize the interaction between the adsorbent and adsorbate, retrieved from the literature.
All the described semiempirical thermodynamics models, as well as other phenomenological approaches, which have been extensively reviewed by Brandani [109], suffer from a lack of thermodynamic consistency in dealing with adsorbates molecules with an appreciable difference in size and do not account in a full predictive fashion for the adsorbateadsorbate and/or adsorbate-adsorbent interactions on the basis of the pure component properties [109]. Despite these semiempirical models being able to exhibit a good fitting capability in view of the large number of phenomenological fitting parameters, they are not suitable for a full predictive approach since these adjustable parameters are not rooted in a rigorous physical background, so their safe use is limited to the condition of the experimental data adopted for the non-linear regression.
To overcome this drawback, a lattice fluid equation of state theory firmly rooted on a statistical thermodynamics background [109,110], aimed at modeling the adsorption thermodynamics of multicomponent fluid mixtures within a rigid adsorbent has been recently proposed in the literature. This approach, known as the Rigid Adsorbent Lattice Fluid model (RALF), has been successfully applied to mixtures of gases and/or vapors within rigid zeolites and MOF systems. To this regard, in the present contribution, we have implemented for the first time the RALF model to investigate the adsorption of pure CO 2 in Cu-BTC. Moreover, in the present investigation we propose an extension of the original RALF model (which is intrinsically a pure mean-field theory), accounting for specific adsorbates-adsorbates and adsorbates-adsorbent interactions. This model, named RALFHB model, has been applied to investigate the adsorption of CH 3 OH in Cu-BTC. In Section 2, we report the fundamentals of the RALF and RALFHB models.

RALF Model
The Rigid Adsorbent Lattice Fluid (RALF) model has been developed by Brandani [109,110] with the aim of describing sorption isotherms of low molecular weight fluid mixtures within a rigid solid adsorbent (such as zeolites and several types of MOFs). Brandani has recently proposed [111] an extension of this model that also accounts for the flexibility of the solid, which has been implemented in the case of a MOF adsorbent that undergoes structural changes in the presence of an adsorbate (breathing transition). In the present contribution, since such kind of structural transitions are not present in the case of binary Cu-BTC/penetrant systems, we take the assumption of rigid solid as reasonable for a quantitative description of the sorption thermodynamics. Therefore, the original version of the Brandani model (indicated hereafter as RALF) has been considered. For the sake of brevity, in the following, we only report the basic equations of this model, referring the interested reader to the original literature for the details regarding the derivation of the equations and the meaning of all the involved variables [109,110].
The RALF model represents an ad hoc extension of the original Sanchez and Lacombe (SL) multicomponent Lattice Fluid model [112][113][114], originally developed to deal with compressible fluid mixtures. The SL model accounts for only self-and cross-mean-field pair interactions and describes the compressibility of the system assuming the presence of empty sites within the lattice. To re-adapt this LF model to the case of a solid adsorbent, Brandani [109,110] assumes that the volume of this mixture V is identified as the apparent volume of the solid V S (i.e., one of the solid including its micropores) in the mixture so that the following relationships hold: where m S , ρ S , ρ, m j , and t represent the mass of the solid, the apparent density of the solid, the density of the adsorbent phase, the mass of j-th component, and the total number of components in the adsorbent phase, respectively. Commonly, the Equation of State, EoS, for the mixtures formed by the adsorbent phase is not available, but it could be possible to use the EoS and/or dilatometric equilibrium data of the pure solid to circumvent the problem. Brandani [109] proposed the following relationship, valid at any temperature, T: where V 0 S is the equilibrium solid volume in vacuum and V ∞ S is the equilibrium solid volume approached at infinite pressure P, in absence of the adsorbates (i.e., pure solid apparent volume) at the given T. In principle, these two values can be provided by an EoS of the pure solid. β T represents the isothermal equilibrium compressibility factor at the given T and in many solid phases, which do not undergo allotropic transformations, can be assumed to be quite independent upon P (this result should be consistently confirmed by the pure solid EoS adopted). Such an approximation is applied in Equation (5). More in general, to avoid the described simplification on β T , the first two terms can be lumped in a single expression V pure S (T, P, N) directly provided by the adopted pure solid EoS (where the i-th component of the vector N represents the number of molecules of component i). Finally, ∆V S represents a correction term accounting for possible solid volume rearrangements induced by the adsorbates and it is, therefore, a function of P, T, and concentration. This term is significant only in the case of flexible solid structures, so that the RALF model assumes that both ∆V S and β T are equal to zero. Therefore, V S (T) is neither a function of adsorbate concentration nor of P; then, according to Equation (4), the volume of the mixture V is not provided by any EoS for the mixture and it is approximated by V S (T), whose value is commonly supplied by a preliminary experimental investigation of the pure solid. To this regard, a common further reasonable assumption is that the value of V = V S is taken also to be independent of T.
In conclusion, based on the rigid network of the solid, it descends the main assumption that the density of the adsorbent mixture is not dictated by the EoS expression provided by the SL theory. However, it is postulated that the Gibbs energy still follows the functional dependence on the state variables of the general out-of-equilibrium constitutive class provided by the SL model. For instance, in an N, P, T ensemble, these are given by P, T, V, and by the vector of the number of molecules of each component in the phase considered, N (hereafter, the symbol N stands for the scalar variable representing the total number of molecules of the phase of interest). Brandani recognizes that this approach represents formally the basis of the Non-Equilibrium Lattice Fluid (NELF) model previously proposed by Sarti and Doghieri [115][116][117]. In fact, the NELF model is an extension of the SL model, in the framework of thermodynamics with internal state variables [118][119][120], to deal with sorption thermodynamics in polymers kinetically frozen in an out-of-equilibrium glassy state. It is worth noting that, in the case of glassy polymers, the frozen value of the mixture volume depends on the non-equilibrium thermomechanical history of the pure polymer sample up to the start of the sorption test and, hence, it can significantly differ (it is commonly higher) from the corresponding equilibrium value at the same P and T [115][116][117]. Conversely, the described simplifications in Equation (5) of the RALF model are rooted in classical equilibrium thermodynamics and the non-equilibrium value of V S = V must be properly intended as an "average" value of the true equilibrium value in the range of T, P, and concentration of interest so that, under the rigid solid assumption, it is expected to be close to the true equilibrium value of the mixture. However, in the development of both the RALF and NELF models, the different rationale that stays behind the assumption of a frozen V is not relevant, since the only significant hypothesis is given by the use of the same functional form of the Gibbs energy provided by the SL framework.
In the SL model and, hence, in the RALF model, N j and r j represent the number of molecules and the lattice sites occupied by a molecule of species j, respectively. Moreover, N 0 represents the number of vacancies (the vacancies are assumed to simply occupy one empty site, so that r 0 is assumed to be equal to 1). Therefore, if v * is the cell volume within the lattice, the total volume V of the mixture formed by t species (including the solid in the case of the adsorbent phase of RALF model) is given by: The close-packed volume of the mixture (V * ) is obtained by setting N 0 = 0 in Equation (6). Therefore, the reduced volume ( v) of the mixture is given by: In Equation (7), ρ represents the actual density of the mixture and ρ * the density of the mixture in the close-packed condition. The volumetric fraction of the j-th species and of the LF empty sites are, respectively (hereafter, superscript L stands for lattice fluid including the empty sites): A first issue related to the use of the SL arises in the need to define a mixing rule for v * . To this regard, in the original SL model, the authors [112][113][114] proposed a mixing rule which allows us to retain the close-packed molecular volume proper of the pure state also within the mixture for each component: In Equation (9), r 0 j represents the number of mers (i.e., of cells) occupied by a molecule of the species j-th in its pure phase (superscript 0 refers to the pure component). By summing upon the t species, it is evident that Equation (9) implies that the close-packed volume of the mixture is additive in terms of the close-packed volume of the pure species. Equation (9) introduce an r j term for each species so that there are t equations but t + 1 mixture variables (including v * ) need to be determined.
To close the problem, in the same series of papers [112][113][114], the authors assumed that the total number of LF cells occupied within the mixture is additive with respect to the LF representation of pure components so that the following equation holds: where x j represents the molar fraction of j-th component. From Equation (10) it follows that: Equation (11) represents the searched operative expressions for the "average" mixture parameter r only as a function of the corresponding pure component parameters and concentration. According to Equations (9) and (10), the close-packed volume is given by: Equation (12) provides the additional function which allows us to close the mixing rules conditions along with Equation (9). From Equations (9)-(11), the following expression is obtained [112][113][114]: where: Equations (13) and (14) represent the needed operative expressions for the v * of the mixture as a function only of pure component parameters and concentration. Finally, the close-packed mass density of the mixture, according to Equation (12), is provided by the following expression: where the last equality is provided by the following relationship for each pure component [115][116][117]: where M w,j represents the molecular weight of j-th species. From Equation (15), the following equation is easily derived: where w j represents the mass fraction of j-th component. Equation (17) represents the operative expression which provides the mixing rule for ρ * as a function of concentration and only of pure component parameters. We remark that the same set of mixing rules provided by Equations (9) and (10) appear both in the RALF and in the NELF models. In a full predictive framework, according to the original SL model, the close-packed volume of each pure penetrant to be adopted in the case of the adsorbent phase can be assumed to be equal to the one adopted in describing the external fluid phase. However, in the case of adsorbed molecules in a rigid solid, due to confinement constraints, it is reasonable (still retaining as mixing rules the Equations (9) and (10) also in the adsorbent phase) to assume that in the solid-penetrant phase the close-packed volume of each pure component, to be adopted for the generic adsorbate, is larger than that in the pure bulk fluid phase [121]. To this regard, Brandani [109,110], in the development of the RALF model, introduced a pair mixture parameter for each penetrant ξ jA , to properly re-scale the v * j to be adopted within the adsorbent phase: In Equation (18), the subscript jA stands for the couple made by the j-th component (different from the solid adsorbent) and by the solid itself (A stands for the adsorbent) and, according to this approach, the parameter ξ jA is expected to be non-negative and, for a given solid adsorbent, to be larger for penetrant molecules with a larger kinetic diameter. Along with the described procedure for v * j , Brandani proposed, according to Equation (18), to re-scale the close-packed density of the component j, ρ * j , whose value is consistently expected to assume a lower value in the case of the adsorbent phase. In particular, the shape of Equation (18) allows us to retain the same value of r 0 j that is used for each adsorbed molecule both in the external fluid phase and in the adsorbent phase.
This represents the main departure of the RALF model from the NELF approach but, since it does not introduce any additional concentration dependence for the mixture parameters, it does not affect the formal derivation of the RALF expressions for the chemical potentials. To this regard, the subscript jA can be formally replaced by the subscript j, just recognizing that in the adsorbent phase, v * j and ρ * j are intended to be the concentration independent quantities provided by Equation (12). This formal simplification is adopted by the expressions of the RALF model reported in the following discussion. By inspecting the SL expressions for the equilibrium chemical potentials, Brandani recognized that the SL model suffers from an intrinsic thermodynamic inconsistency: in multiphase phases differently from the pure case, when the volume (at a given T and concentration) diverges positively, the model does not recover the expression of the chemical potentials of the ideal gas state. This is a well-documented issue [122][123][124][125][126][127][128] which arises from the common set of mixing rules of v * adopted in the SL literature. To this regard, different approaches to correct this inconsistency have been proposed [122][123][124][125][126]. Moreover, it has been pointed out that a more complex lattice fluid model, the so-called Non-Random Hydrogen-Bonding (NRHB) model [129,130], is conceived in a way that this drawback is overcome [127,128].
The first step in the development of the RALF model, working in an N, V, T ensemble, is the introduction of a posteriori correction of the expression of the combinatorial term of the SL Helmholtz energy which allows us to disgard of the mentioned inconsistency. Starting from this corrected expression, it was found [109] that (within an additive constant and terms that are only dependent upon T and that are not involved in the determination of the chemical potential expressions as well as of the EoS) the Gibbs energy of the SL model can be written as: Starting from Equation (19), the corresponding residual Gibbs energy takes the form: (20) where: represents the close-packed volumetric fraction of the j-th component and n represents the total number of moles. We remark here that expressions (19) and (20) represent only a first step in the development of the RALF model. At the end of this section, we will discuss a further correction of the Gibbs energy expression introduced by Brandani to properly re-adapt the SL model approach to the case of an adsorbent solid phase (see Equation (30)).
In Equations (19) and (20), the reduced temperature T and pressure P are, respectively, given by: Moreover, z is the non-equilibrium compressibility factor that, in the framework of the SL model, is given by [112][113][114]: where n t represents the total number of moles. If the mixture density would be allowed to attain its equilibrium value, then one would have: where: In the SL model, RT * and P * represent the characteristic energy and the characteristic energy density of the mixture which, in this model, are only ascribed to "mean-field" interactions. The following expression relates these two characteristic parameters: Such an expression holds for any multicomponent phase including the pure component. Based on Equations (13), (14) and (26), only one mixing rule concerning T * or, alternatively, P * needs to be specified. In the framework of the version of SL on which the RALF model is rooted, the following mixing rule is adopted for P * [109,110]: where: In Equations (27) and (28), when j = k, the corresponding P * jj represents the characteristic pressure P * j of the component j and the corresponding k jj is equal to zero. Moreover, from Equation (28) it follows that k kj = k jk .
In conclusion, the proposed mixing rule for P * introduces a set of dimensionless parameters, k jk , each one defined for a given couple of components of the multicomponent phase of interest. Therefore, for the case of t components, there is a set of k jk parameters (with j = k) in a number of t(t − 1)/2, which represent, in general, a set of optimization parameters of the model for the system of interest along with the set formed by the (t − 1) ξ jA parameters. According to the proposed mixing rule, any k kj (with j = k) can be considered as characteristic of the specific couple of components involved and, regarding the penetrant components, in the RALF model it is still assumed to be the same both in the adsorbent phase and in the external fluid phase, as it occurs in the SL model.
In addition to the correction proposed by Brandani for the characteristic parameters v * j of the adsorbate, P * j and T * j should also be corrected. However, to minimize the number of fitting parameters and based on the physical meaning of the term RT * j , in the original implementation of the RALF model it has been assumed that T * j takes the same value associated with the bulk fluid phase, where no geometrical constraints are present. Moreover, following the line of thought of the original SL model, it has been also assumed that the relationship (28) still holds in the adsorbent phase for each adsorbent. Consequently, when the correction given by Equation (18) for v * j in the adsorbent phase is implemented, the corresponding P * j will be consistently lower by the factor 1 1+ξ jA .
In conclusion, the RALF model introduces (t − 1)/ t 2 + 1 optimization parameters (the set of k jk with j = k and the set of ξ jA ). In order to adopt the model in a predictive fashion for the case of adsorption in Cu-BTC, each one of the described optimization parameters should be retrieved through a non-linear regression procedure to fit solubility data of any other system containing the couple of components associated.
We remark here that, since the RALF model has been developed for multicomponent systems, the sub-case of a pure fluid component is included in all its expressions. In particular, the expressions for G of multicomponent fluid mixtures, for the chemical potentials, and for the compressibility factor z naturally recover the case of pure component in the limit of its molar fraction approaching to 1, and all the mixing rules adopted converge consistently to the corresponding pure component value.
In the expressions reported above, each summation includes all the t components of the phase of interest. It is worth noting that the RALF model does not consider for any possible adsorption process occurring on the surfaces of the solid crystals. In fact, the latter contribution is expected to be, as a first approximation, negligible in comparison to the largely prevailing contribution of adsorption within the bulk phase of the solid crystals. Hence, in the case of the adsorbent phase, Brandani [109] assumes that, according to the homogeneous approach of the LF model proposed, the solid is formed by a single huge molecule (so-called 'single bulk monocrystal' assumption). Operatively, in using the model for data correlation purposes, the surface adsorption contribution as well as the polycrystallinity dispersion are lumped, in the form of an effective "averaged" additional contribution, with the bulk adsorption contribution. In certain cases, this approximation could result in some inconsistent values of best fitting parameters and, consequently, could determine some deviations from the experimental data when the model is used for a fully predictive approach, as it will be discussed in Section 4.
Having assumed that the solid is formed by a single molecule, the number of solid molecules N s in the continuous LF approach tends to zero. Consequently, its molecular mass diverges. Finally, since the cell volume v * s , the mixture cell volume v * , and the closepacked pure solid volume V * 0 S are described by real numbers, the number of lattice sites occupied by the solid in the pure lattice (r 0 s ) as well as in the mixtures (r s ) diverges. In fact, from Equation (9), the following relationship is derived: where ρ * s represents the close-packed density of the solid component. Equation (29) allows us to determine the limit condition for the term r s N s from which, in principle, it is possible to re-express also the rN term in the limit condition. Alternatively, the limit expression of rN can be straightforwardly obtained by observing that Equation (15) can be recast as: Then, from Equations (4), (7) and (30), rN can be expressed by the finite positive quantity given by: Moreover, as indicated by Brandani, in the limit condition of N s → 0 , the SL expression of the residual Gibbs energy given by Equation (20) needs to be properly re-adapted for the adsorbent phase. In fact, in the derivation of the RALF model [109,110], G R,SL has been modified considering that, since N s → 0 and the solid is assumed perfectly rigid, the adsorbent component cannot contribute to the combinatorial term but has the effect of reducing the volume available to the adsorbent molecules within the lattice, as also expressed by the proposed scaling of the v * j of the penetrants. To properly account for this issue, in the RALF model the combinatorial part of the G R is ad hoc modified with respect to the SL expression adopted for the external fluid phase (given by Equation (20)), according to the following expression [109]: In Equation (32), the index i replaces the index j of Equation (20). In fact, i refers only to the t − 1 adsorbate components while the solid component is explicitly labelled with subscript s. It is worth noting that the expression of the compressibility factor, z, in the case of the adsorbent phase, is still provided by Equation (23).
Equation (32), under the hypothesis that the solid "frozen" volume is coincident with the adsorbent phase volume, represents the starting point of the RALF model to consistently derive the expression of the residual chemical potential contribution of any penetrant present in the adsorbent phase, which is the only contribution required for the calculation of the solubility isotherms. As for the penetrant chemical potential in the external fluid phase, the expressions of the residual penetrant chemical potentials are provided by performing the derivative of Equation (20). In the following, we briefly discuss the operative expressions of the penetrant chemical potential contribution required for sorption thermodynamics calculations. Regarding the external fluid phase, the expression of the equilibrium residual chemical potentials on a molecular basis for the component k, µ R,eq k , can be calculated by determining the derivative of the equilibrium expression of residual G in an N, P, T ensemble, i.e., where V(T, P, N) is given by the EoS (Equation (24)) and G R (T, P, N, V(T, P, N)) is given by Equation (20). Hereafter, the superscript eq stands for equilibrium. The functional dependence V(T, P, N) in the model under consideration is not available in a closed form. However, we recall that the equilibrium volume V of the mixture in a T, P, N ensemble is provided by the minimization condition of G towards the phase volume (i.e., the EoS) [112][113][114] and that the dependence of G from the mixture volume is only present in the G R term [109,131,132]. Therefore, according to the derivative chain rules, the following relationship is obtained: In Equation (34) which is obtained by coupling the following general non-equilibrium expression (which can be calculated in closed form): with the EoS (Equation (24)). In particular, the following expression holds in the case of the external fluid phase: where the subscript F stands for fluid phase and z andρ F are obtained by solving Equations (24) and (25). In the case that the external phase is mono-component, Equation (36a) becomes: Regarding the adsorbent phase, the volume of the mixture has been fixed to be in an out-of-equilibrium state and the solid is assumed insoluble in the external fluid phase. These are the same conditions hypothesized in the case of the NELF model. Consequently, Brandani [115][116][117]127] infers that in the RALF model, the phase equilibrium conditions between the external fluid phase and the adsorbate phase are dictated by the same set of equations determined for the case of the NELF model when dealing with the phase equilibrium between a multicomponent fluid phase and the glassy polymer-penetrant phase kinetically locked in an out-of-equilibrium state.
The NELF phase equilibrium conditions impose that for each adsorbate, the equilibrium chemical potentials in the external phase are equal to the non-equilibrium (NE) chemical potentials in the polymer-penetrant phase. In particular, the appropriate nonequilibrium chemical potential required in the NELF model is provided for the k-th penetrant by [115][116][117]127]: where in this case, V is the assumed out-of-equilibrium fixed value of mixture volume, taken coincident with the out-of-equilibrium value of the pure polymer right before the sorption test.
Invoking the same approach of the NELF, the RALF model phase equilibrium condition is dictated by the following set of equations expressed in terms of residual chemical potentials [109,110]: where y k and x k refers to the molar fraction of the k-th adsorbate in the external fluid phase and in the adsorbent phase, respectively. In Equation (38), µ R,eq k,F is given by Equation (36a,b) coupled with Equation (2). µ R k is referred to the adsorbent phase and it is still provided by Equation (35) but, differently from the equilibrium external phase case, the G R adopted for its calculation is provided by Equation (32) coupled with the non-equilibrium expression of z (Equation (23)). The associated value of volume to be used in these expressions is given by Equation (4).
We remark that in the case of the NELF model, the condition of phase equilibrium based upon the non-equilibrium expression of the chemical potentials provided by Equation (37) has been proven to be thermodynamically consistent with the constraints imposed by the second law of the thermodynamics on the whole biphasic system considered based on the thermodynamics with internal state variables [115][116][117]. In this respect, a well-established constitutive class for the function expressing the kinetic evolution of the volume mixture, provided by a viscoelastic model, was adopted [115][116][117]. Conversely, in the RALF approach the equivalent condition given by Equation (38) does not rely upon such kind of analysis, but it is inferred based on the same formal conditions used in the NELF model regarding the "frozen" mixture volume.
In order to close the phase equilibrium problem when using the RALF model, in the following we need to focus on the operative calculation of µ R k . As a consequence of the discussed limit condition N s → 0 , some comments are in order regarding the calculation of the derivative in Equation (35) in the case of the adsorbent phase. In principle, starting from Equation (32) coupled with Equations (23) and (24), the limit condition just imposes that rN is expressed by Equation (31) and that N approaches the total number of adsorbate molecules N p . In this way, the functional dependence G R (T, P, N, V) can be formally re-expressed as G R T, P, N p , m s , ρ s , V = G R T, P, N p , V , where N p is the subvector of N referred only to the adsorbate components, and the last identity of Equation (31) has been applied to lump m s and ρ s in terms of V. The expression of G R T, P, N p , m s , ρ s , V can be, therefore, equivalently adopted in the derivative of Equation (35) in order to obtain µ R k , since neither m s nor ρ s depends on N p . In this way, the operative expression of µ R k is provided by: with: Incidentally, we observe that, in the case of a pure fluid, Equation (40) collapses into the original version of SL EoS [112][113][114]. We recall that in Equations (39) and (40), the volume V is assigned according to Equation (4) so that ρ is directly provided as: In Equations (39) and (40), the close-packed volumetric fraction of the solid is given by: where the limit condition N s → 0 imposed to obtain the last member is equivalent to divide for N and N p . Finally, according to Equation (29), the term r s x s can be expressed as: where the last identity follows from the recalled limit condition. From Equation (31), in the limit condition, r can be expressed as: Equations (43) and (44) introduce the ratio N p /m s : an intensive quantity expressing the total number of molecules of adsorbate per mass of solid. This is in line with the proposed approach in which, being that the solid mass is not soluble in the external phase, the adsorbate concentration vector can be expressed per solid mass.
In the case of a single adsorbate, in Equations (39) and (40), the conditions x 1 = 1 and r 1 r = ϕ 1 = 1 − ϕ s hold, so that Equations (39) and (40) result in the single equation:

RALFHB Model
The RALF model is based upon an LF statistics which accounts only for "mean-field" interactions between the component molecules, assuming that the empty sites do not contribute to the interaction energy. In the following, we propose an extension of the original RALF model aimed at describing the possible presence of self-and cross-specific interactions (i.e., Hydrogen-Bonding and/or Lewis acid-bases interactions) between the components, referred in the following as the RALFHB model. On the basis of the discussed formal correspondence between the NELF and the RALF models, the proposed extension follows the same procedure adopted by Mensitieri et al. [127,133] in extending the original "mean-field" version of the NELF model to the case of sorption thermodynamics of glassy polymer-penetrant mixtures displaying possible self-and cross-HB interactions. This latter extension is referred to as the NELFHB model [133]. In the following, we briefly report the fundamentals of the RALFHB approach. Full details on the procedure for the derivation of NELFHB model (and consequently, of the RALFHB one), which is rooted in thermodynamics with internal state variables, are reported in [127,133].
The first Step of the NELFHB, as well as of the RALFHB model, requires a general non-equilibrium expression of G accounting for both mean-field and specific interactions. This expression is provided by Panayiotou and Sanchez (PS) [134] who proposed a model based upon the factorization of the partition function between a "mean-field" contribution and an excess contribution related to specific interactions. In fact, this factorization results in two additive contributions that provide the expression for G: one "mean -field" LF term (G LF ), which is similar to the one provided by the SL model but that can be in principle provided by any "mean-field" model, and an additional term (G HB ) accounting for specific interactions, which is provided by the Veytsman [135] statistics [134,135]. Its extension to deal with sorption thermodynamics of polymer-penetrant mixtures characterized by the presence of strong interactions and kinetically locked in an out-of-equilibrium glassy state is, hereafter, referred to as the NELFHB model. A similar approach can be used to develop a model for adsorption thermodynamics in a rigid adsorbent phase where strong interactions occur. In this way, the RALFHB model is derived from the RALF approach in the same way in which the NELFHB model is obtained from the NELF approach.
Regarding the RALFHB model, consistently with the RALF model case, we have adopted for G LF the expression given by Equations (20) and (32), respectively, for the external fluid phase and for the adsorbent phase, taking the original G HB term of the PS model for both the phases. The expression of G HB can be found in the original literature of PS model [134] and reads: In Equation ( Hereafter, G 0 αβ represents the molar Gibbs energy of formation of the specific interaction between a proton donor of kind α and a proton acceptor of kind β. We remark that the energies of formation of HB must be considered as an excess energy contribution with respect to the mean-field energy which is described by the LF contribution.
Therefore, the HB contribution to G introduces in principle m × n HB formation energies, G 0 αβ , where each one is defined by: With U 0 αβ , V 0 αβ , and S 0 αβ , respectively, representing the molar internal energy, the molar volume, and the molar entropy of formation of the specific interaction between a proton donor of kind α and a proton acceptor of kind β. According to a well-established approach [127], the assumption V 0 αβ = 0 is commonly adopted, so that, according to Equation (47), the HB contribution introduces 2(m × n) energetic parameters, provided by the two sets of U 0 αβ and S 0 αβ terms. The values of such parameters regarding self-and cross-HB interactions between the adsorbates, in principle, can be obtained by regressions of VLE equilibrium data according to the version of PS model implemented for the external phase. Therefore, in an N, P, T ensemble the general non-equilibrium expression of G can be formally expressed as [133]: where: G R,eq ≡ G R N, P, T, V(N, P, T), N H αβ (N, P, T) (51) and V(N, P, T) and N H αβ (N, P, T) are, in principle, provided by the solutions of Equation (49a,b). Equation (50) follows from Equation (49a,b) since the dependence of G from V and N H αβ is only involved in the calculation of the G R term [132][133][134]. The last member of Equation (50), which can be expressed in a closed form, coupled with Equation (49a,b), therefore, provides the operative expression of µ R,eq k for the equilib-rium external fluid phase, thus circumventing the problem derived from the fact that the solution of Equation (49a,b) is available only numerically.
To close the sorption thermodynamics modeling, we need to focus now on the adsorbent phase. As discussed, in framework of the RALF model as well as of the NELF model, it is assumed that the adsorbent phase does not follow the equilibrium condition of the LF model adopted. In the case of the NELFHB model and, consequently, of the RALFHB model, this results in assuming that Equation (49a,b) does not hold. In principle, under the assumption of a "frozen" adsorbent phase, by following the procedure used to develop the NELFHB model in the case where the specific interactions are accounted for by the additional set of internal variables N H αβ , one needs to provide the fixed out-of-equilibrium values of V and of N H αβ . Commonly, the latter information is not available so that an "instantaneous equilibrium" assumption occurs [127,133]. In fact, in the NELFHB model and, consequently, in the proposed RALFHB model, in view of the microscopic scale regarding the nature of the HB contacts that allows for a very fast rearrangement of the interactions, it is assumed that N H αβ still follows the equilibrium condition dictated by Equation (49b), but in correspondence with the fixed out-of-equilibrium mixture volume V (i.e., it is assumed that Equation (49a) does not hold) [127,133].
For the adsorbent phase, under the N s → 0 assumption, the general out-of-equilibrium residual Gibbs energy of the model, G R is now given by: where N H αβ (N P , P, T, V) is, in principle, provided by the solutions of Equation (49b), V is the fixed mixture volume, and G R,LFHB (N P , P, T, V, N H αβ ) is the sum of the G LF term given by Equation (32) and of the G HB term of the original PS and NELFHB models.
The expression for G R (N P , P, T, V) obtained under the HB instantaneous equilibrium assumption, therefore, represents the general out-of-equilibrium expression of a model which displays the volume V as the only internal state variable. Consequently, the problem has been traced back to the well-established NELF framework already adopted in the RALF model, and the out-of-equilibrium residual chemical potential of k-th adsorbate, µ R k , is provided by derivation of G R (N P , P, T, V), according to the expression (35).
However, the solution of the sub-set (49b) in correspondence with the fixed V is in general available only numerically so that G R (N P , P, T, V) is not available in a closed form. To circumvent this problem, the calculation of µ R k can be performed as: that has been obtained by applying Equation (49b). Equation (53) results from the fact that the dependence of G upon N H αβ involves only the residual term G R . The last term in Equation (53) can be now expressed in a closed form. Operatively, the value of µ R k is provided by Equation (53) coupled with Equation (49b).
Again, since the instantaneous equilibrium assumption formally brings back the formulation of the model to the original NELF framework, the set of Equation (38) properly re-casted in terms of the RALFHB model, still dictates the phase equilibrium state (i.e., the sorption thermodynamics of the adsorbates within the rigid adsorbent phase).
In the following, we show the operative expression of the residual chemical potentials. To this aim, it is worth noting that: where G id represents the ideal gas term of the Gibbs energy and G (R,LF) represents the residual Gibbs energy of the mean-field LF contribution, whose expression is provided by Equations (20) and (32), respectively, in the case of the external fluid phase and the adsorbent phase. Therefore, in the calculation of the residual equilibrium chemical potentials of the k-th penetrant in the external fluid phase, from Equation (50) it follows: where: and Equations (56) and (57) need to be coupled with Equation (49a,b) and the subscript F has been added to underline that this is referred to as the fluid phase. Analogously, for the adsorbent phase, from Equation (53) it follows that: where: As before, Equations (59) and (60) need to be coupled with Equation (49b) with the assigned volume V provided by Equation (4). It is worth noting that the last identities in Equations (56) and (59) follow from the fact that the general non-equilibrium expressions of G LF in each phase do not depend upon N H αβ [127,133,134]. Consequently, the formal expressions of µ R,LF k,F and of µ R,LF k coincide with the corresponding ones of the RALF model reported in the previous section, i.e., Equations (36a) and (39) for the fluid phase and the adsorbent phase, respectively. In the following, we report the operative expressions of µ R,HB k,F and µ R,HB k along with the minimization conditions represented by Equation (49) which are required to close the phase equilibrium problem. We observe that Equations (57) and (59) are formally given by the same expression. Hence, by following the same procedure of refs. [133,134], one obtains: where d k a and a k β represent, respectively, the number of donors of kind α and of acceptors of kind β present on a molecule of species k-th, and we have also that: The minimization conditions represented by Equation (49a) imply that: where, in obtaining the HB contribution it can be followed by the procedure reported in [127,133].
In the case of a pure fluid, Equation (63) collapses into the equation of state of the PS model [134]. Finally, Equation (49b), according to [127,133], implies that: for each α = 1, 2, . . . , m and β = 1, 2 . . . , n On the basis of the previous discussion, for the external fluid phase, Equation (61) is intended to be coupled with Equations (63) and (64), and, for the adsorbent phase, Equation (61) is intended to be coupled only with Equation (64) while the fixed volume V is provided by Equation (4).
We conclude this section addressing an issue which concerns the application of both the RALF model and its extension, the RALFHB model. In fact, from inspection of the shape of the expressions of the adsorbate chemical potentials in the adsorbent phase (i.e., Equation (39) in the case of the RALF model and Equation (61) and Equation (39) coupled with Equation (64) in the case of the RALFHB model), it is evident that the simultaneous dependence on the out-of-equilibrium volume, V = V S , and on the mass of solid, m s , can be expressed as a dependence only on the solid mass density, ρ S . This result is expected since the chemical potentials are intensive properties. Operatively, the following relationship can be used: In Equation (65), ρ * s and V m represent, in the framework of the RALF approach, the solid skeletal density and the pore volume of the solid per mass of solid itself, respectively. Both the adsorbent parameters can be measured experimentally, as it will be discussed in Section 4. Finally, we observe that the dependence on ρ s is, in turn, only expressed as a dependence on the ρ of the adsorbent phase mixture, which can be simply calculated as:

Dual-Site Extension of the RALF Approach: The RALF-DS and RALFHB-DS Models
In this section, we deal with the extension of the original RALF approach to the case of a heterogeneous rigid adsorbent displaying different kinds of adsorbent cages. Indeed, as detailed in the Introduction section, this is the case of the Cu-BTC MOF.
Firstly, we briefly examine the extension of the RALF approach proposed by Brandani [136] to address the case of a dual cage system, which results in the RALF dual solid model (RALF-DS). As it will be evident in the following, this approach can be easily extended to any multi-cage case once the required information on the heterogenous structure of the rigid adsorbent as well of the corresponding set of model parameters are available.
The fundamental assumption is that the two kinds of adsorbent cages contribute in parallel to the adsorption process. In the DS approach, these two different families of cages are ascribed to two different "fictive" solids. The whole adsorption isotherm of the actual rigid adsorbent for a given external fluid phase is then calculated as the averaged sum of the contributions associated to each fictive solid. To this regard, each fictive solid can be, in turn, modeled by using the RALF model or the RALFHB model in the case that specific interactions are expected. Following the approach of Brandani [136] to develop the RALF-DS model, we propose here the RALFHB-DS model in which the behavior of each fictive solid is described using the RALFHB model.
The implementation of the RALF-DS model introduces two additional parameters as compared to the homogeneous RALF model, namely, the mass fraction of fictive solid I with respect to the total solid mass of the actual solid (indicated as ω I ) and the fraction of pore volume of fictive solid I with respect to the total pore volume of the actual solid (indicated as f vI ). The density of each fictive solid is given by [136]: In Equation (67) where m k,tot refers to the total mass of adsorbate k-th in the actual solid and m s,tot is the total mass of the actual solid. Finally, the subscripts I and II refer to fictive solids I and II, respectively.

Summary of the Theoretical Approaches
Empirical and semiempirical modeling have been proposed in the literature to describe sorption in porous materials in general, and in MOFs in particular. These approaches provide a quick and effortless evaluation of the adsorption phenomena. Relevant examples are the Langmuir-Freundlich (LF) isotherm, characterized however by a thermodynamic inconsistency in the low-pressure region, the Virial-Langmuir (VL) model which allows us to gather some information regarding the nature of the adsorbate-adsorbent interaction, and the Dubinin-Astakhov (DA) model which provides a fair prediction of multicomponent isotherms Semiempirical thermodynamics models, despite their good fitting capability, suffer from a lack of thermodynamic consistency in dealing with adsorbates molecules with an appreciable difference in size and do not account in a full predictive fashion for the adsorbate-adsorbate and/or adsorbate-adsorbent interactions. Moreover, they are not suitable for a fully predictive approach since their adjustable parameters are not rooted in a rigorous physical background.
To provide a robust background on the mechanisms taking place during the gas adsorption on MOF materials and an interpretation of their relative vibrational spectra, over the last decade different ab initio computational methods have been implemented which allowed for the identification of the adsorption sites. For Cu-BTC frameworks, many contributions have been produced in this direction using DFT, GCMC, and MD calculations, but a large computational effort is required for these calculations with MOFs.
To overcome this drawback, a lattice fluid equation of state theory firmly rooted on a statistical thermodynamics background aimed at modeling the adsorption thermodynamics of multicomponent fluid mixtures within a rigid adsorbent has been recently proposed in the literature. This approach, known as RALF, has been successfully applied to mixtures of gases and/or vapors within rigid zeolites and MOF systems. The RALF model is an ad hoc extension of the original SL multicomponent Lattice Fluid model but, differently from SL theory, in the RALF approach, the volume of the adsorbent-adsorbate system is not provided by an EoS and is assumed to be constant. In view of its structure, the RALF model is intrinsically a pure mean-field theory. Then, based the on RALF approach, we propose here an extension of this model to deal with cases where specific adsorbates-adsorbates and adsorbates-adsorbent interactions need to be accounted for, the RALFHB model.
Finally, to deal with the case of heterogeneity of adsorption sites, we have shown how one can easily extend both the RALF and the RALFHB models by considering different kinds of adsorbent cages and introducing the so-called 'Dual-Solid' or 'Multi-Solid' models.

Experimental Section
In this section, we summarize the information on the material and procedures we adopted to obtain in our lab, some experimental data that are discussed in Section 4, and the literature data.

Materials
Basolite ® C300 [Copper benzene-1,3,5-tricarboxylate in powder form was provided by Sigma-Aldrich (Milano, Italy). According to the supplier, the specific surface area was in the range of 1500-2100 m 2 /g, the bulk density was 0.35 g/cm −3 , and the average particle size of the crystalline powder was 15.96 µm. Carbon dioxide with 99.99% purity was supplied by SOL S.p.A. (Monza, Italy). Chloroform with 99.8% purity and Potassium Bromide (KBr) windows 2.0 mm thick with a diameter of 13 mm were supplied by Sigma-Aldrich (Milano, Italy).
A Cu-BTC dispersion was prepared from a 5.0 wt% mixture of Cu-BTC in chloroform stirred for 2 h and further sonicated for 30 min. Samples for spectroscopy measurements were prepared by casting a few drops of the Cu-BTC dispersion on a KBr window allowing solvent evaporation for 2 h at room temperature. Afterwards, the window was slightly pressed at 5 bar for compacting the MOF layer, thus optimizing the FTIR transmission measurement.

FTIR Spectroscopy
The FTIR measurements were performed under gas flowing using a modified Linkam cell, THMS350V (Surrey, UK), equipped with temperature control (83-623 K) and a vacuum system. The cell was connected through service lines to a mass-flow-controller [MKS Type GM50A (Andover, MA, USA)] to set the CO 2 flux, while a solenoid valve regulated the downstream pressure. The system was equipped with a Pirani vacuometer and an MKS Baratron 121 pressure transducer (full scale 1000 Torr, resolution 0.01 Torr, accuracy ±0.5% of the reading) (Andover, MA, USA). Further details on the experimental apparatus are reported in [137]. The equipment allowed for the in situ activation of the sample at 150 • C under a vacuum and the collection of isothermal data at different temperatures. The diffusion cell was coupled to a Spectrum-100 FTIR spectrometer (Perkin-Elmer, Norwalk, CT, USA), equipped with a beam splitter made of a thin Ge film supported on KBr plates and a wide-band DTGS detector.

Analysis of Cu-BTC/CO 2 Systems by FTIR Spectroscopy
A comparison between the spectra of Cu-BTC after activation (red trace, see Figure 2) and after equilibration at 35 • C with an external CO 2 pressure of 150 Torr (blue trace, see Figure 2) highlights the signals produced by the guest molecule at 2338 cm −1 (O=C=O antisymmetric stretching, ν 3 ) and at 656 cm −1 (O=C=O bending, ν 2 ).
The absorbance of the ν 3 band is linearly dependent on the amount of sorbed CO 2 , as demonstrated by correlating the sorption isotherm monitored spectroscopically with that from an independent volumetric measurement (see Figure S1, Supplementary Materials [137]). This correlation allows for the calibration of the photometric observable according to the Beer-Lambert relationship and validates the quantitative analysis by spectroscopic means. The ν 3 bandshape provides further information on the molecular environment. A pronounced shoulder is detected at about 2324 cm −1 ; previous studies [138,139] demonstrated that this feature is a (ν 3 + ν 2 ) − ν 2 hot-band enhanced by Fermi resonance with the nearby fundamental. It does not originate from a CO 2 population distinct from that producing the main signal and, after suitable resolution, can be neglected in the bandshape analysis. The main component at 2338 cm −1 displays a lower resolution and a full-width-at-halfheight (FWHH) considerably larger than those observed in various host/guest systems of the same kind (see, for instance, a comparison with the polyetherimide/CO 2 system represented in Figure S2, Supplementary Materials [140]). This effect is indicative of a multicomponent bandshape with an unresolved fine structure. Resolution-enhancement approaches are needed to analyze such complex profiles: we adopted two-dimensional correlation spectroscopy (2D-COS), a well-established tool very effective in diffusion studies for investigating weakly interacting systems [141][142][143]. The resolution enhancement brought about by 2D-COS originates from the spreading of the spectral data over a second frequency axis, coupled with the vanishing of the asynchronous correlation intensity for signals evolving at the same rate. In addition, the technique provides information on the evolution of the system as a function of the perturbing variable (CO 2 concentration, in the present case) [144].
The asynchronous map relative to isothermal measurements at 0 • C is represented in Figure 3. It displays a well-resolved four-peak pattern that reveals the presence of three components at 2337, 2346, and 2333 cm −1 . These correspond to probe populations that interact with different chemical environments in the nanostructure; the diverse interactions promote the separation of the ν 3 signal in three components and induce a difference in terms of molecular mobility of the guest molecules that grants the detection by 2D-COS. This conclusion has been confirmed by the analysis of the bending mode in the 670-640 cm −1 range [137] and has been associated with the occurrence of three specific adsorption sites, theoretically predicted by first-principles studies, that are characterized by distinct values of the interaction energy [144].
The FTIR analysis also reveals that the CO 2 /HKUST-1 interactions are weak. In fact, the perturbation induced to the ν 3 mode is unable to produce a sizable resolution of the components, not even in the form of shoulders or bandshape asymmetry. The occurrence of multiple components can only be evidenced by applying 2D-COS spectroscopy. Fur-thermore, the ν 3 peak position is close to that observed in weakly interacting systems (see Figure S2, Supplementary Materials).

Modeling Sorption Isotherms of Low Molecular Weight Compounds in Cu-BTC
In this section several experimental results available in the literature for sorption isotherms of gases and vapors in Cu-BTC and their interpretation with RALF, RALF-DS, RALFHB, and RALFHH-DS are presented and reviewed.
To the aim of modeling the adsorption of CO 2 and CH 3 OH in Cu-BTC by using the RALF approach, one needs to provide structural and energetic parameters of the model. In particular, the RALF model requires pure LF parameters of the adsorbent solid: the close-packed density (ρ * s ), T * s , and P * s . Since this set of parameters cannot be found through the fitting of dilatometric equilibrium data of the solid, the common procedure consists of estimating the pure solid parameters through the simultaneous non-linear regression of solubility data of several penetrants [137]. According to this procedure, in the present contribution we have directly implemented the model to fit simultaneously the solubility of light weakly interacting gases (O 2 , CH 4 , N 2 , and CO 2 ) and specific interacting vapor (CH 3 OH).
Based on the Cu-BTC structure and the kinetic diameters of the adsorbates considered (reported in Table 3), three different cages are, in principle, involved in the adsorption process of the investigated light gases and vapor.  [146] In the case of weakly interacting gases, the observed type I adsorption behavior [147] suggests that the penetrants independently explore and access the different cages of the framework under observation and are not affected by the presence of the different guesthost interactions related to the double Cu 2+ paddle-wheel. Consequently, a 'single solid' RALF model can be applied to describe the mean-field adsorption process within the whole solid. Conversely, in the case of CH 3 OH, according to [148], specific interactions are expected to occur between Cu 2+ and hydroxyl group of methanol. Due to heterogenous spatial distribution of the copper doublet (L2 with respect to L3 and S1), the methanol molecules probe different energetic environments when exploring the polar cages (L2) and the apolar cages (L3 and S1) [58]; indeed, this is confirmed by the observed type IV adsorption behavior [147]. To account for this heterogenous solid behavior, in the case of methanol, the 'dual solid' approach has been implemented, combining the RALFHB model, used for describing the polar cages contribution, with the pure mean-field RALF model, used for describing the apolar cages contribution. Cross-and self-specific interactions have been modeled by assuming that each hydroxyl group of methanol has 1 proton donor d 1 1 and 1 proton acceptor a 1 1 ; each copper of Cu-BTC has 2 proton donors d S 2 . For the sake of comparison, we have also implemented, in the case of methanol, the RALF-DS model in which both solids have been modeled disregarding HB contribution.
The set of optimization parameters adopted in the simultaneous fitting procedure is formed by: , ω I , and f V,I . The corresponding optimized values of these parameters are reported in Table 4. We remark that, to reduce the number of fitting parameters, a single value of T * s , P * s , and ρ * s have been adopted for both the 'single solid' and the 'two fictive solids' (I and II) models. Therefore, the whole fitting procedure is coupled for all investigated adsorbates. The energetic heterogeneity of the system in the case of methanol is taken into account for both the specific interaction contribution (present only in the solid I) and different mean-field interaction parameters of the 'two fictive solids' (k I CH 3 OH−S and k I I CH 3 OH−S ). In addition, the cage's heterogeneity, which is significant in the case of methanol, is accounted for through the two volume correction parameters (ζ I CH 3 OH−S and ζ I I CH 3 OH−S ) as well as through the intrinsic double solid parameters ω I and f V,I . In order to close the modeling, one needs to provide the skeletal solid density which in the RALF model framework corresponds to ρ * s , the pore volume, V m , of the whole MOF expressed as pore volume per mass of the whole solid, and the RALF lattice fluid characteristic parameters, T * , P * , and ρ * , of each penetrant which actually correspond to the Sanchez-Lacombe parameters, as discussed in Section 2.1. In addition, the self-specific interaction parameters U 0 CH 3 OH−CH 3 OH and S 0 CH 3 OH−CH 3 OH for methanol are also required. The list of the adopted values of these parameters is reported in Table 3. Regarding the methanol, two sets of parameters are reported in Table 5: the first one [112] refers to a pure mean-field RALF model and, consequently, the characteristic parameters coincide with the Sanchez-Lacombe ones, while the second set is associated with the RALFHB model and it has been retrieved in the present work by a non-linear regression of liquid-vapor equilibrium data.  Figure 4 shows the comparison between the experimental data for vapor pressure of methanol as a function of temperature and of density of methanol at liquid-vapor equilibrium taken from [149] and the model fitting using the optimized values of fitting parameters of the RALFHB model. It is evident that the model exhibits an excellent fitting capability in the whole range of the temperature investigated. In fact, the RALFHB model reproduces simultaneously the saturation densities of liquid and vapor phase along with the saturation pressure, provided that the data are sufficiently far away from the critical transition. We remark that the only fitting parameters are provided by the characteristic lattice fluid parameters, in fact the hydrogen-bonding formation parameters have been fixed according to the literature of PS [134]. To this regard, as discussed in Section 2.2, the RALFHB parameters for the external fluid phase coincide with the ones of the PS model in the case of the single penetrant system. The last required input parameters for the modeling are ρ * s and ρ s . The first one is identified with the skeletal density of the solid and, according to its crystalline structure, is taken in the literature to be equal to 2.685 g cm −3 [150]; the second parameter, ρ s , usually is not provided directly by the experiments but it can be calculated in the case of a single solid by Equation (65) while, in the case of the RALF-DS model, is provided by Equation (67) for each fictive solid. To this regard, the parameter V m (usually determined experimentally by pycnometry) is required to close the problem. V m depends on the way the Cu-BTC sample is prepared. For the sets of data regarding the weakly interacting gases, V m is equal to 0.79 cm 3 g −1 [48] and for Cu-BTC/CH 3 OH sets of data, V m is equal to 0.731 cm 3 g −1 [58].
Once all the input parameters have been provided, simultaneous fitting of all the sets of adsorption data considered can be performed. In Figures 5-9, a comparison between the experimental adsorption data for the weakly interacting gases and the optimized predictions provided by the simultaneous fitting procedure is reported.     It is evident that the modeling approach exhibits an excellent fitting capability both in the case of weakly interacting gases (in particular, reproducing well the type I behavior) and in the case of specific interacting Cu-BTC/CH 3 OH systems (in particular, reproducing well the stepwise type IV behavior commonly observed in a semi-log scale). To better elucidate the DS approach, in Figure 10, the two additive contributions associated to the fictive solids I and II at 70 • C (similar results, omitted for brevity, have been found also at the other investigated temperatures) are explicitly reported. It is worth noting that the polar solid I dominates the adsorption process at low pressures approaching a saturation plateau. This is expected in view of the strong cross-specific interactions occurring between the hydroxyl group of methanol and Cu 2+ , whose total number is limited by the stoichiometry of the MOF. Once the solid I contribution approaches its saturation limit, the stepwise contribution of the apolar solid II becomes significant allowing the overall solubility curve of type IV of the whole solid to reproduce. It is worth noting that the optimized value of S 0 CH 3 OH−S is significantly low in comparison to the values commonly observed for hydrogen-bonding-specific interactions (see for instance the value of S 0 CH 3 OH−CH 3 OH in Table 3). Indeed, this result is in line with the observation that the interaction between OH − of methanol and Cu 2+ is not associated to an angular constraint of the rotational mobility, differently from the case of a classic hydrogen-bonding. In the framework of Veystman [26,27] statistics, the correlated reduction in rotational mobility is associated with the entropy of formation of the specific interaction. Then, in the case of non-directional-specific interactions, as in the case of interest here, the value of the entropy of formation of the interaction is still expected to be negative but negligible. In conclusion, the values of U 0 CH 3 OH−S and S 0 CH 3 OH−S point to a quite athermal Gibbs energy of formation of cross-specific interactions, which is more negative as compared with the common values of hydrogen-bonding Gibbs energy of formation (as reported in the specific table in ref. [127]). Indeed, this is reasonable in Lewis' acid-base involving ions and polar groups.
To better understand RALFHB model results, one can analyze the self-and crossspecific interactions occurring in the system. For instance, regarding the isotherm reported in Figure 10, at 70 • C for the solid associated with the polar cage (L2) at the lowest pressure investigated (1.148·10 −4 MPa), the model predicts that around 77.5% of the adsorbentspecific interacting sites are involved in the cross interaction with the methanol hydroxyl, while at the highest pressure investigated (0.02486 MPa), this value raises to around 91.2%. This result is consistent with the shape of the adsorption contribution of solid I which represents the polar cages (see Figure 10). Regarding the ratio between the methanol selfand cross-specific interactions in the solid I, its value exhibits a crossover in the range of pressures investigated spanning from 0.47 to 1.49 MPa; this result is in line with the observed trend of the saturation value of the MOF-specific interacting sites and, in addition, it highlights that more non-cross interacting methanol molecules are located in the polar cages as the external pressure is increased. The modeling approach also predicts that the same qualitative trend of the methanol self-specific interacting molecules as a function of the pressure is observed in the solid II, which represents the apolar cages. Similar results are observed at the other temperatures investigated and, in particular, by reducing the temperature, an increase in both self-and cross-specific interactions of the methanol in each solid occurs and this effect is more prominent for the interactions characterized by the lower entropy of formation. In fact, in the case of cross-specific interactions the behavior is quite athermal.
For the sake of comparison, we have also implemented a simultaneous fitting procedure of the whole set of data disregarding the contribution of specific interactions, i.e., by applying the pure mean-field RALF model also in describing the polar cage in the case of the Cu-BTC/CH 3 OH system, but still retaining the DS approach. Figure 11 reports the comparison between the experimental data and the optimized fitting results for the Cu-BTC/CH 3 OH systems. It is evident that disregarding the HB contribution does not allow us to properly reproduce quantitatively the experimental data, particularly at low pressure, where the contribution of Cu-BTC-CH 3 OH cross-HB interaction is expected to be significant. Regarding the weakly interacting gases, the correlated predictions are quite like the ones reported for the previous simultaneous fitting case (i.e., the ones reported in Figures 5-8). For the sake of brevity, the table reporting the optimized parameters for the simultaneous fitting performed with the pure mean-field RALF model and the figures showing the comparison of the experimental data of the weakly interacting gases with the optimized prediction curves are reported in the Supplementary Materials (see Figures S3-S7  and Table S1).
To further test the robustness of the RALF type of approach, based on the values of the optimized fitting parameters reported in Table 4, we have implemented the model in a fully predictive fashion to reproduce five CO 2 adsorption isotherms of two Cu-BTC systems which differ from the ones adopted in the previous simultaneous fitting procedure. These two Cu-BTC structures are characterized by different specific volume porosities, V m , still retaining the same ρ * s implemented in the case of sets of data of the fitting procedure, so that ρ s assumes different values.
For the CO 2 adsorption isotherms at 10 • C, 20 • C, 45 • C, and 70 • C, V m is equal to 0.57 g cm −3 [50] and ρ s is calculated according to Equation (65). In the case of the CO 2 adsorption isotherms at 35 • C, a commercial Cu-BTC (Basolite ® C300) which displays ρ s = 1.26 g cm −3 [56] has been used. Figure 12 proves the excellent agreement of RALF model predictions (implemented with the single solid framework) with the experimental data of CO 2 adsorption isotherms at 10 • C, 20 • C, 45 • C, 70 • C (data taken from [50]), and 35 • C (data taken from [137]). In conclusion, by properly accounting for the actual specific porosity of the different Cu-BTC samples, the model can consistently predict the adsorption isotherms in a wide range of temperatures and pressures. To better elucidate the model performance regarding the type IV behavior exhibited by the methanol isotherms, we have also investigated the fitting capability of the RALFHB model assuming a simple 'single solid' framework. To this aim, we have performed a simultaneous non-linear regression of the adsorption isotherms at 40 • C, 50 • C, 60 • C, and 70 • C. As an example, in Figure 13 the results of the best fitting procedure of the experimental data at 40 • C are reported. It is evident that accounting for specific interactions in the 'single solid' approach is not enough to reproduce the type IV behavior, indicating that the combined energetic and structural heterogeneity of the adsorbent material should be properly considered using a multi-solid approach.

Conclusions
An overview of experimental adsorption data of low molecular weight compounds in Cu-BTC is provided, also analyzing some semiempirical and theoretical models used to interpret adsorption thermodynamics.
The theoretical approach based on the original RALF approach has been described in detail and extended to account for the contribution of specific self-and cross-interactions by introducing the new RALFHB model. This general framework allows us to account also for the heterogeneity of adsorption sites of the adsorbent materials by properly modifying the structure of the RALF and RALFHB models considering a 'multi solid' structure of the adsorbent material. The case of adsorption of weakly interacting gases in Cu-BTC can be dealt with simply by using the RALF theory while the thermodynamics of adsorption of methanol can be properly addressed only by considering a dual solid (DS) structure, using the RALF-DS model to account for sites heterogeneity or, better, using the RALF/RALFHB-DS model to account for both sites heterogeneity and the occurrence of specific interactions for polar sites. In this way, a theoretical 'platform' is made available based on the original RALF framework, well-suited for the interpretation of adsorption thermodynamics in MOFs and, in particular, in Cu-BTC.
In fact, this theoretical approach has been proven to be very effective both in fitting the adsorption data of some gases and of methanol vapor in a specific Cu-BTC sample and in providing excellent predictions for adsorption in other Cu-BTC specimens displaying a different specific volume porosity.