<figure><img src="image_1.png" /></figure>

# Efficient Modeling of Water Adsorption in MOFs Using Interpolated Transition Matrix Monte Carlo

Bartosz Mazur,* Lucyna Firlej, and Bogdan Kuchta*

<figure><img src="image_2.png" /></figure>

Cite This: ACS Appl. Mater. Interfaces 2024, 16, 25559–25567

<figure><img src="image_3.png" /></figure>

Read Online

## ACCESS

<figure><img src="image_4.png" /></figure>

Metrics & More

<figure><img src="image_5.png" /></figure>

Article Recommendations

<figure><img src="image_6.png" /></figure>

Supporting Information

**ABSTRACT:** With the specter of accelerating climate change, securing access to potable water has become a critical global challenge. Atmospheric water harvesting (AWH) through metal-organic frameworks (MOFs) emerges as one of the promising solutions. The standard numerical methods applied for rapid and efficient screening for optimal sorbents face significant limitations in the case of water adsorption (slow convergence and inability to overcome high energy barriers). To address these challenges, we employed grand canonical transition matrix Monte Carlo (GC-TMMC) methodology and proposed an efficient interpolation scheme that significantly reduces the number of required simulations while maintaining accuracy of the results. Through the example of water adsorption in three MOFs: MOF-303, MOF-LA2-1, and NU-1000, we show that the extrapolation of the free energy landscape allows for prediction of the adsorption properties over a continuous range of pressure and temperature.

This innovative and versatile method provides rich thermodynamic information, enabling rapid, large-scale computational screening of sorbents for adsorption, applicable for a variety of sorbents and gases. As the presented methodology holds strong applicative potential, we provide alongside this paper a modified version of the RASPA2 code with a ghost swap move implementation and a Python library designed to minimize the user's input for analyzing data derived from the TMMC simulations.

**KEYWORDS:** water adsorption, atmospheric water harvesting, computational screening, molecular modeling, metal-organic frameworks, transition matrix Monte Carlo

<figure><img src="image_7.png" /></figure>

## INTRODUCTION

In the coming years, as climate change accelerates, access to potable water may become a luxury. One of the solutions to this pressing issue lies in the extraction of drinking water directly from the atmosphere in a decentralized manner using atmospheric water harvesting (AWH).¹⁻⁷ Materials and even devices that enable this process, have already been demonstrated with metal-organic frameworks (MOFs) emerging as particularly promising sorbent candidates.⁸⁻¹¹ MOFs offer advantages such as relatively low synthesis costs, the ability to tailor their pores' shape and size to the desired profile of the water adsorption isotherm, and stability in an aqueous environment, allowing for multiple adsorption-desorption cycles.

Despite the continuous advancement of sophisticated experimental methods to study adsorption phenomena with ever deeper understanding at the microscopic/atomistic scale, numerical simulations offer a level of resolution, understanding, and control over the investigated system that the physical experiment have yet to achieve. Therefore, they remain a remarkably powerful tool for delving into the microscopic world, and for modern material design. This is particularly true when examining MOFs and processes in which they are involved.

Given the modular nature of MOFs, the number of potential structures is theoretically infinite and at least difficult to estimate. The CSD MOF database contains over 12 000 experimentally synthesized MOFs,¹² CoRE MOF database exceeds 14 000,¹³ and the hMOF database comprises more than 51 000 hypothetical MOFs.¹⁴ Conducting experimental investigation on such a vast array of materials and the selection of the most suitable one for a specific application is prohibitively costly in terms of time and resources. Therefore, the utilization of large-scale high-throughput numerical screening methods becomes indispensable.

Although water is one of the most common molecules, essential to many chemical, physical, biological, and technological processes, fully comprehending its interaction with surrounding heterogeneous structures remains challenging. A better understanding of the intricate mechanism of water adsorption in MOFs is essential for optimizing the adsorption

Received: February 15, 2024

Revised: April 25, 2024

Accepted: April 29, 2024

Published: May 6, 2024

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>
process in current materials and designing new ones, not only for water harvesting but also for processes where the presence of water influences the adsorption,¹⁵,¹⁶ separation,¹⁷ or catalysis¹⁸ of other molecules.

Among various computational techniques, Monte Carlo simulations in a grand canonical ensemble (GCMC) have proven to be particularly valuable for modeling the adsorption of small molecules on surfaces and in porous structures.¹⁹⁻²² However, challenges persist in water adsorption simulations, including long convergence times and frequently observed discrepancies between experimental and numerical results.²³⁻²⁷ We believe that they are mostly the consequence of two factors: (i) incomplete description of the water interaction with its heterogeneous surroundings, often neglecting strong polarization effects that significantly influence adsorption²⁸ (interactions are typically described using rigid, nonpolarizable models); and (ii) high energetic barriers between states, resulting from very strong interactions within water itself, primarily through hydrogen bonding; these barriers often cause the GCMC simulations to become trapped in local energy minima.²⁹ In this work, our aim is to address the latter issue; the former requires separate investigation.

To solve the problem of long-living metastable states, instead of standard GCMC algorithm, we employed the grand canonical transition matrix Monte Carlo (GC-TMMC) method using NVT + *ghost swap* approach. Recent studies have shown that this technique significantly enhances the computational efficiency of water adsorption simulations.²⁹,³⁰ It consists in conducting simulations in canonical ensemble in which only *ghost* trial insertion or deletion of a molecule (a swap move) is performed but never accepted; only the transition probability is recorded for further analysis. The NVT + *ghost swap* method has been successfully used in various simulations, including adsorption of rigid and flexible molecules under supercritical conditions,³¹ adsorption of small molecules such as methane in rigid³² and flexible³³ frameworks, and adsorption of water in rigid frameworks.²⁹

One of the key advantages of GC-TMMC lies in its ability to study the properties of the adsorbed gas, even when the gas is present in more than one phase.³⁴ This allows the analysis of metastable phases, which are responsible for the observation of hysteresis loops on the adsorption–desorption isotherms.³⁵,³⁶ As hysteresis is a nonequilibrium phenomenon, its accurate determination depends on the duration of experimental observation.³⁷ Consequently, the unambiguous experimental determination of hysteresis is challenging, complicating in turn verification of numerical results. However, GC-TMMC simulations allow direct determination of both the limit of stability of the adsorbed phase and pressure of equilibrium phase transition.³⁴

The major drawback of the NVT + *ghost swap* method is the necessity to conduct *N* independent simulations, covering the full spectrum of possible macrostates, ranging from an empty to fully saturated framework. Due to substantial porosity and volume available for fluid in MOFs, the required number *N* of NVT simulations remains high, often exceeding 1000.

Recently, a modification of the NVT + *ghost swap* procedure, referred to as the C-map method and dedicated to screening materials in adsorption applications, has been developed.²⁹ In this approach, at specified pressure and temperature (*pT*) conditions, a TMMC simulation is conducted at predetermined uptake, and the probabilities of molecule insertion and deletion are calculated and accumulated. Consequently, if the

system exhibits a higher probability of accepting insertion, it suggests that at the given *pT* conditions the equilibrium loading is higher than the one at which calculations were performed. Similarly, if the probability of particle removal exceeds that of insertion, the equilibrium loading has a lower value. By conducting such simulations for various loading values, one can determine the equilibrium loading range for a specific material under defined *pT* conditions.

A simplified version of this method was recently used³⁸,³⁹ to rapidly assess the hydrophobicity or hydrophilicity of MOFs. The simulation was performed only for a single arbitrarily chosen loading value. In such a situation, assuming a perfectly stepped water adsorption isotherm, if the probability of accepting the insertion of a molecule is greater than the probability of its removal, the material under given *pT* conditions is fully filled (hydrophilic); otherwise, it is considered empty (hydrophobic).

In this paper, we propose a methodology with a computational cost comparable to that of the C-map method, yet it also yields the complete water adsorption isotherm. We introduce an effective interpolation scheme wherein to obtain full adsorption isotherm, direct simulations are conducted only for selected macrostates. This approach reduces the total number (and time) of simulations by 2 orders of magnitude while maintaining good agreement with the isotherms calculated from simulations performed for all macrostates. This represents a significant improvement in methodology, especially for applications such as water harvesting processes and other temperature/pressure swing adsorption processes, where the value of interest is working capacity (i.e., the deliverable amount of water or other gas, in a single adsorption–desorption cycle).³,⁴⁰ This value can be determined from the isotherm or isobar of adsorption.

Then, we present, for the first time, an extrapolation of the free energy landscape using these reduced data sets, from which we estimate the adsorption uptake at temperatures that were not simulated directly. The use of the reweighting and extrapolation techniques results in the significant advantage of the GC-TMMC method over the GCMC method—the ability to determine the adsorption properties of a sorbent (in this case, three selected MOFs) across a continuous range of temperatures and pressures based on calculations conducted *at only one* temperature and pressure. This methodology represents a pioneering approach that enables simultaneously fast, accurate, and efficient high-throughput screening of MOFs for water adsorption applications.

<figure><img src="image_10.png" /></figure>

## METHODS

Here, we outline the fundamental assumptions of the grand canonical transition matrix Monte Carlo (GC-TMMC) procedure. For a detailed thermodynamic description of the method, we direct the reader to references.³⁴,⁴¹,⁴²

The main objective of the GC-TMMC method is to calculate the macrostate probability distribution (MPD), instead of direct calculation of the ensemble averages. To achieve this goal, the simulation is biased to sample all macrostates with equal probability. In this way, macrostates with relatively high free energy may be explored. When investigating adsorption phenomena, it is convenient to use the number of particles in the system, denoted as *N*, as a macrostate variable. The MPD in the grand canonical ensemble is then given by
$$ \Pi(N; \mu, V, \beta) = \frac{\exp(\beta\mu N)Q(N, V, \beta)}{\Xi(\mu, V, \beta)} \quad (1) $$

where $\mu$ is chemical potential, $V$ is system volume, $\beta = 1/k_B T$ ($T$ being temperature and $k_B$ Boltzmann constant), and $Q(N,V,\beta)$ and $\Xi(\mu,V,\beta)$ are respectively canonical and grand canonical partition functions. To determine MPD, each time the system visits a macrostate $N$, the so-called collection matrix (C-matrix) is updated with the unbiased probabilities of accepting a swap move ($p_{acc}$), according to

$$ C(N_o \rightarrow N_n) = C(N_o \rightarrow N_n) + p_{acc}(o \rightarrow n) \quad (2) $$

$$ C(N_o \rightarrow N_n) = C(N_o \rightarrow N_n) + [1 - p_{acc}(o \rightarrow n)] \quad (3) $$

where the labels $o$ and $n$ refers to old and new configurations for attempted move. Then, the probability transition matrix is computed by normalization of the C-matrix:

$$ P(N \rightarrow N + \delta) = \frac{C(N \rightarrow N + \delta)}{\sum_{\Delta \in \{-1,0,1\}} C(N \rightarrow N + \Delta)} \quad (4) $$

The $\Pi(N; \mu, V, \beta)$ can be then calculated using relation:

$$ \ln \Pi(N + 1; \mu, V, \beta) = \ln \Pi(N; \mu, V, \beta) + \ln \left[ \frac{P(N \rightarrow N + 1)}{P(N + 1 \rightarrow N)} \right] \quad (5) $$

Using the following protocol: first, an arbitrary value is given to $\ln \Pi(N_{min})$, then values of $\ln \Pi$ for next $N$ are calculated sequentially. The minimum and maximum particle numbers are set to ensure the sampling of the entire domain of interest, from empty ($N_{min} = 0$) to fully saturated system (with $N_{max}$ calculated from the volume of the system available for adsorption; for details read section Simulation Details in Supporting Information). To ensure that the entire domain of interest has been sufficiently sampled, one should always check whether $\Pi$ has a sufficiently low value in $N_{max}$:

$$ \max[\ln \Pi(N_{liquid})] - \ln \Pi(N_{max}) > \text{tolerance} \quad (6) $$

where $N_{liquid}$ corresponds to the macrostates in the liquid domain. In our simulations, we used the tolerance of 10 as the states of relative probability of $e^{-10}$ do not contribute significantly to the ensemble averages.

To calculate the ensemble average of any physical quantity $A$ in the grand canonical ensemble, we use the equation:⁴³

$$ \langle A_{\alpha} \rangle = \frac{\sum_{N \in \alpha} A(N, V, T) \Pi(N; \mu, V, T)}{\sum_{N \in \alpha} \Pi(N; \mu, V, T)} \quad (7) $$

where the sums are calculated over the macrostates that belong to phase $\alpha$ and $A$ denotes value of interest which is accumulated for each visited macrostate $N$. The GC-TMMC simulations performed at one $\mu$ value can be easily recalculated for any other $\mu'$ value by simple histogram reweighting using relation:⁴²

$$ \ln \Pi(N; \mu', V, T) = \ln \Pi(N; \mu, V, T) + \beta N(\mu' - \mu) \quad (8) $$

In this way, from a single set of calculations (performed at single $\mu$) we can obtain the full isotherm at infinite resolution. For a more detailed description of the identification of phases and the calculation of thermophysical properties from MPD, we refer the reader to the paper by Siderius et al.³⁴

To compute MPD, it is essential to sample all relevant macrostates to gather adequate transition statistics. Typically, this is achieved in a multiple-macrostate approach, where a single simulation sweeps through a range of adsorbed particles. To sample low-probability states, a biasing function is employed, such as in the Wang–Landau (WL) algorithm.⁴⁴,⁴⁵ It is also possible to use a single-macrostate approach, where $N$ simulations in NVT ensemble are conducted, which artificially imposes equal sampling of a domain and allows for the utilization of the TMMC principles without employing a biasing function.²⁹,³¹,³³ In this case, transition probabilities are calculated using the so-called *ghost swap* move. This move resembles a standard swap move, involving a trial insertion or deletion of a molecule, but it is never accepted: only transition probabilities are recorded. Consequently, the total number of molecules remains constant through simulation. Hatch et al.⁴⁶ demonstrated that single-macrostate simulations are less efficient than multiple-macrostate simulations, mainly due to the lack of microstate sampling. However, we show that with single-macrostate simulations, it is possible to skip sampling certain macrostates and interpolate transition probabilities to obtain the full MPD with reduced simulation cost. At the same time, in this study, we do not undertake a comprehensive analysis of the differences in efficiency among various approaches, which we discuss in detail in the last part of the *Results and Discussion* section.

In this work, NVT + *ghost swap* simulations were conducted for each $N$ within the macrostate range using modified version of RASPA2 code.⁴⁷ In the *Supporting Information*, we also compare the results obtained from these simulations with those obtained using the WL/TMMC method to eliminate the risk of referring to erroneous results. Following the simulation, we reduced the NVT + *ghost swap* simulation data sets by factors of 20, 50, and 100. Table 1 provides the exact numbers of

<table><caption>**Table 1.** Number of Macrostates for which Direct Transition Probability Calculations Were Performed</caption><thead><tr><th>reduction factor</th><th>MOF-303</th><th>MOF-LA2-1</th><th>NU-1000</th></tr></thead><tbody><tr><td>reference</td><td>745</td><td>889</td><td>1200</td></tr><tr><td>20</td><td>37</td><td>44</td><td>60</td></tr><tr><td>50</td><td>15</td><td>18</td><td>24</td></tr><tr><td>100</td><td>7</td><td>9</td><td>12</td></tr></tbody></table>

macrostates in the reduced data sets for each system studied: MOF-303,⁴⁸ MOF-LA2-1,⁹ and NU-1000.⁴⁹ Transition probability values for macrostates not included in the reduced data sets were calculated by linear interpolation between the nearest states in the simulation. A visualization of such an approach for MOF-303 is presented in Figure 1. The selection of these materials was based on their spectrum of isotherms, which closely approximates the typical isotherms of water adsorption in high-potential MOFs for use in AWH.⁴⁰ This encompasses adsorption by continuous pore filling (MOF-303), adsorption with stepped isotherms with minimal presence of metastable states (MOF-LA2-1), and adsorption with stepped isotherms with pronounced and long metastable states (NU-1000). The rationale behind the choice of materials is discussed in detail in the Choice of Porous Materials section in *Supporting Information*.

To extrapolate the free energy landscape obtained at single simulation temperature to other temperatures, we followed the
<figure><img src="image_11.png" /><figcaption>**Figure 1.** Macrostate representation for MOF-303. Macrostates for which the transition probabilities were explicitly calculated are represented by open symbols; macrostates with interpolated transition probabilities are marked with a dark line connecting open points. The light bold lines indicating the reference macrostates are provided for better comparison of the used macrostates reduction.</figcaption></figure>

method published by Mahynski et al.⁵⁰,⁵¹ The ln Π(N) is expanded in a Taylor series:

$$ \ln \Pi(N; \beta) = \ln \Pi(N; \beta^0) + \frac{\partial \ln \Pi(N; \beta)}{\partial \beta} |_{\beta^0} \Delta \beta + \frac{1}{2!} \frac{\partial^2 \ln \Pi(N; \beta)}{\partial \beta^2} |_{\beta^0} \Delta \beta^2 + \dots \quad (9) $$

where Δβ = β - β₀ and β₀ is the inverse simulation temperature. In the grand canonical ensemble, the first two terms of the Taylor expansion simplify to

$$ \frac{\partial \ln \Pi(N; \beta)}{\partial \beta} = \mu N - \langle E \rangle \quad (10) $$

$$ \frac{\partial^2 \ln \Pi(N; \beta)}{\partial \beta^2} = \langle E^2 \rangle - \langle E \rangle^2 \quad (11) $$

where E is the potential energy and ⟨E⟩ is the ensemble average potential energy in the canonical ensemble. In this study, all extrapolated isotherms are calculated using only the first order Taylor expansion, since the next terms converged much more slowly. In most cases, it caused significant noise, which considerably affected the obtained results, as discussed in detail in the Supporting Information.

## RESULTS AND DISCUSSION

**Definition of the Problem and Proposed Methodology.** We begin by justifying the choice of the grand canonical transition matrix Monte Carlo (GC-TMMC) method for simulations of water adsorption. Our case study is water adsorption in NU-1000 at temperature T = 298 K. First, we conducted two series of calculations using the GCMC method: (i) starting from an empty system, we calculated the adsorption isotherm, and (ii) starting with a pre-equilibrated system containing 528 molecules/unit cell (corresponding to a fully saturated system), we calculated the desorption isotherm (Figure 2, left panel). Both curves served as reference for subsequent analysis. Due to the slow convergence of GCMC water adsorption simulations, extensive calculations were conducted, comprising 2 × 10⁶ cycles for adsorption simulations and at least 1 × 10⁶ cycles for desorption, with the longest desorption simulation spanning 90 days. Despite the considerable duration of these simulations, the system has

<figure><img src="image_12.png" /><figcaption>**Figure 2.** (left) Isotherms of water adsorption in NU-1000 at 298 K simulated using the GCMC (2 × 10⁶ and 1 × 10⁶ cycles for adsorption and desorption, respectively) and GC-TMMC (10⁵ cycles) methods. The solid yellow line represents a track of stable states; the dashed yellow line represents metastable branches of the isotherm. GCMC adsorption simulation was started with an empty system. GCMC desorption simulation was started with system containing 528 molecules/unit cell that was already equilibrated in NVT simulation (for water cluster formation). In GCMC simulations, during both adsorption and desorption, insertion and removal of molecules were allowed. (right) Free energy profiles at relative pressure indicated by a vertical line of the corresponding color on the plot on the left.</figcaption></figure>

not reached the equilibrium state. In contrast, using the GC-TMMC method, we conducted only one set of calculations at T = 298 K and p = 3200 Pa, comprising only 10⁵ cycles in production runs. From these calculations, using reweighting, we were able to determine both the low- and high-density branches of the isotherm (referred to as gas and liquid branches, respectively), phase stability ranges, and equilibrium-transition pressure. Additionally, from the same simulation, we calculated the free energy of a system at various pressures (Figure 2, right panel). This allowed us to explain both the mechanism of water adsorption in NU-1000 and why GCMC simulations did not deliver an equilibrium isotherm.

As presented in Figure 2 (left), the GC-TMMC adsorption isotherm aligns with the GCMC gas phase branch across the entire pressure range, even as the latter becomes metastable. Similarly, the liquid phase branch of the GC-TMMC desorption isotherm coincides with the GCMC throughout the pressure range where the liquid phase occurs. As the metastable liquid phase vanishes, the system is emptying and the desorption branch merges with the gas phase branch, the only phase that exists at low pressure.

In principle, the GCMC isotherms should follow only stable states in both adsorption and desorption simulations. However, this is not the case because simulations often become trapped in local energy minima. To illustrate this phenomenon, we calculated (from GC-TMMC simulations) the free energy profiles at various pressures: (i) when only the gas phase is present, (ii) when the liquid phase is metastable, and (iii) when the gas phase is metastable (see Figure 2, right panel). We observe the emergence of two minima—one at low density (ld) and the other at high density (hd), corresponding to uptakes of ~30 and ~500 molecules per unit cell, respectively. Even at saturation pressure, starting the GCMC simulation from an empty system does not lead to the increase of the adsorbed amount. In fact, despite a free energy difference of approximately 400 kJ/mol favoring the hd state, the energy barrier (~200 kJ/mol) between the ld and hd states is roughly 80 times higher than the thermal energy k_B T at 298 K.
<figure><img src="image_13.png" /><figcaption>**Figure 3.** Isotherms of water adsorption at 298 K calculated for three MOFs: (top) MOF-303, (middle) MOF-LA2-1-ZUS, and (bottom) NU-1000, using data sets reduced (left) 20 times, (middle) 50 times, and (right) 100 times. The reference isotherms (calculated using full data sets) are colored with a slightly wider blue line for better visibility.</figcaption></figure>

K, causing the system to remain blocked in the local minimum in the *ld* state. Similarly, when the simulation begins with a fully saturated system, at *P*/*P*<sub>0</sub> = 0.43, it remains trapped at the local minimum in the *hd* state, despite a free energy difference of about 1200 kJ/mol in favor of the *ld* state. Again, the barrier between the *ld* and *hd* states exceeds the thermal energy. As the pressure decreases, the local minimum associated with the *hd* state disappears, and the simulation initiated from the *hd* state at *P*/*P*<sub>0</sub> = 0.26 easily reaches the equilibrium state at *ld* as there is no energy barrier hindering the process.

Consequently, we constate that, compared to standard GCMC simulations, the GC-TMMC method offers increased accuracy, reduced simulation time, and deeper understanding of the thermodynamic information of the system under investigation.

**Isotherms from Reduced Data Sets.** As previously mentioned, the main limitation of the NVT + ghost swap method is the number of simulations needed to calculate the macrostate probability distribution (MPD), which serves as the basis for extracting all thermodynamic properties of a system. We show that it is possible to conduct simulations at reduced number of macrostates and subsequently interpolate the missing ones from the directly calculated values.

To assess the accuracy of interpolation of transition probabilities, we first conducted simulations across all macrostates for three systems: MOF-303, MOF-LA2-1, and NU-1000. Then, during postprocessing, we reduced the number of macrostates for which transition probabilities are known directly 20, 50, and 100 times (Table 1). The resulting isotherms are presented in Figure 3. For all three MOFs, the isotherms calculated using data sets reduced 20 and 50 times

show good or near-perfect agreement with the reference isotherms. Upon further reduction in the number of simulations (by 100 times), we observe larger (albeit still acceptable) discrepancies: (i) for MOF-303, the uptake is slightly underestimated at lowest pressure and overestimated at highest pressure; (ii) for MOF-LA2-1 the saturation amount and location of the step on isotherm are accurately predicted, with minor differences in gas branch; however, these differences are not significant enough to impact the analysis of the results; (iii) for NU-1000, differences are more pronounced—while the saturation volume is reproduced correctly, the uptake in the gas phase and the step pressure slightly deviates from the reference values.

We observe a distinct trend in which, as the system size expands, and consequently the energy barrier between states increases, the reduction in direct simulations results in a noticeable degradation of the calculated adsorption isotherms. One factor contributing to this result may be that as the size of the system increases, the range of macrostates that separate the phases broadens. Consequently, more simulations fall within the range of low probability states compared to the number of simulations that sample high probability states.

However, when searching for optimal materials for atmospheric water harvesting applications, we prioritize specific characteristics, including (i) the relative humidity (or *P*/*P*<sub>0</sub>) at which the material adsorbs water, (ii) shape of the isotherm (a step-shaped isotherm being most desirable), (iii) the volume of water that can be collected in a single adsorption-desorption cycle, (iv) the hydrothermal stability of the material (to ensure long-term operation), and (v) components used for adsorbent synthesis (preferably, non-
toxic, and abundant). The first three features can be confidently obtained directly from short NVT + ghost swap simulations. According to Figure 3, even with data sets reduced 100 times (thus calculated using only 7, 9, and 12 macrostates for MOF-303, MOF-LA2-1, and NU-1000, respectively), these characteristics are in perfect or good agreement with reference values. In the Supporting Information, we present the calculation of working capacity for NU-1000, identified as the worst performing material in this study. We show that the error for the calculated values is only 1.16%, 3.60%, and 8.51% for simulations performed with 60, 24, and 12 macrostates, respectively (Table S2).

A direct comparison of the computational times of GCMC and GC-TMMC is not straightforward, given their different approaches of phase space sampling. In GCMC simulations, a single simulation provides an average uptake under specific *p*T conditions. In contrast, in GC-TMMC, a series of simulations for different uptake values are conducted, from which an isotherm with infinite resolution is derived through pressure reweighting in postprocessing. Consequently, the time required for a single simulation in GC-TMMC depends on the number of molecules present in the system (thus size of the system), rather than the selected pressure. For this reason, for some hydrophilic MOFs (such as MOF-303 or MOF-LA2-1) that do not present a challenge to GCMC simulations, the computational cost required to calculate a few isotherm points might be similar to performing a series of GC-TMMC simulations, assuming the interpolation scheme proposed in this work is used. However, predicting the simulation outcome prior to its execution is not possible. Furthermore, it is only after completing a GCMC simulation that one can determine whether the simulation will be challenging or not. Additionally, if the GCMC simulation becomes trapped at a local energy minimum, fluctuations in properties as a function of the number of MC cycles may falsely indicate equilibrium (Figure S7). In conclusion, both methods can achieve similar performance for the most hydrophilic MOFs; however, GC-TMMC offers better control over simulation convergence. Furthermore, other benefits of the GC-TMMC method, such as infinite isotherm resolution, access to the system's free energy, ability to explore metastable states, and temperature extrapolation, should also be taken into account.

**Temperature Extrapolation.** One of the major advantages of the TMMC method is its capability to extrapolate the free energy landscape computed at a single temperature to other temperatures through simple postprocessing of the data.<sup>50,51</sup> Here, for the first time, we evaluate the accuracy of this extrapolation technique while simultaneously reducing the number of simulations. We limited the analyzed temperatures to the 298 K–343 K range, as this represents an operating range for water adsorption in MOFs.

Figure 4 shows isotherms and isobars calculated using the NVT + ghost swap method and the extrapolation technique. The left panel presents water adsorption isotherms at 343 K, extrapolated from data collected at 298 K with a reduced number of macrostates. In all cases, the extrapolated isotherms agree well or nearly align with the reference isotherms computed directly at 343 K across the entire pressure range. When the number of macrostates decreases, some minor differences appear in the low-density range and in the metastable branches; however, the most important characteristics of the isotherm (such as the pressure of the *ld* to *hd* transition or total uptake) remain in good agreement with

<figure><img src="image_14.png" /><figcaption>**Figure 4.** Extrapolated isotherms (*T* = 343 K, left) and isobars (at 1.7 and 3.5 kPa for MOF-303 and MOF-LA2-1, and at 5.0 and 10.0 kPa for NU-1000, right). The reference isotherms (in blue) were calculated for each MOF in direct simulations at *T* = 343 K. All other curves were extrapolated from simulations performed at *T* = 298 K using the number of states indicated in the legend.</figcaption></figure>

simulation performed directly at 343 K. Considering the results obtained for 298 K with interpolation of transition probabilities (as shown in Figure 3), a similar degree of decrease in accuracy is observed during simultaneous extrapolation and interpolation. Consequently, within the studied range of temperatures, the observed error is mainly due to interpolation of the results rather than their extrapolation.

This observation represents a substantial advancement in the screening capabilities of simulations. So far, most screening procedures employing classical simulation methods have focused on extracting a single material characteristic (in most cases, the maximal loading at a selected pressure). With the NVT + ghost swap method, we can derive thermodynamic data for a diverse range of temperatures and pressures from just a few simulations with practically infinite resolution, given that extrapolation and pressure reweighting are obtained during data postprocessing.

Most industrial applications of adsorption phenomena use temperature/pressure swing adsorption (T/PSA) procedures, where cyclic changes in pressure and/or temperature prompt alternating adsorption and desorption of fluid. The working capacity of a particular sorbent is evaluated as the amount of fluid that can be recovered (or harvested, in the case of water adsorption from the atmosphere) in one adsorption-desorption cycle. This quantity can be easily calculated from the adsorption isotherms/isobars at various temperatures/pressures of interest as a difference in the average uptake between any two (*p*₁*T*₁), (*p*₂*T*₂) state points.
To illustrate this approach, the right panel of Figure 4 shows the adsorption isobars for the discussed systems. Let us consider determining the working capacity of the MOF for water adsorption at T = 298 K and p = 3.5 kPa and the desorption at T = 333 K and p = 1.7 kPa. By comparing MOF-303 and MOF-LA2–1, we can conclude that under identical operating conditions MOF-LA2–1 yields approximately twice as many molecules of water per unit cell as MOF-303 (60 molecules/unit cell vs 30 molecules/unit cell, respectively).

**NVT + Ghost Swap Convergence Time.** The main constrain in water adsorption screening studies using classical simulation methods is the extremely long simulation time. This issue is solved within the NVT + ghost swap method. Therefore, we also verified (following previously discussed interpolation technique) how the number of simulation cycles affects the isotherm profile. The results are presented in Figure 5: isotherms derived from data collected after 5 × 10³ cycles of

<figure><img src="image_15.png" /><figcaption>**Figure 5.** Isotherm of water adsorption in NU-1000 at 298 K calculated with 5 × 10³ production cycles and variable number N of macrostates used in the calculation of direct transition probabilities. Solid and dashed lines correspond to stable and metastable states, respectively.</figcaption></figure>

simulations closely reproduce the reference isotherm, even using the minimal number of simulations. As presented in Figure 2, the isotherm derived from the GCMC simulations failed to converge to equilibrium values despite conducting simulation with a number of cycles 3 orders of magnitude higher than those presented in Figure 5. This suggests that the presented approach (short, low-demand simulations) can serve for preliminary materials selection in large-scale screening studies, allowing more accurate (longer) calculations to be performed on already narrowed selection of sorbents.

The efficiency of the presented NVT + ghost swap method compared to the multiple-macrostate TMMC method remains an open question, particularly when considering the advantages of probability interpolation. Results presented by Hatch et al.⁴⁶ suggest that for adsorption in a repulsive porous network, single-macrostate simulations might be 3 orders of magnitude less efficient than double-macrostate simulations, due to the lack of microstate sampling. However, their study focused on gas adsorption in a purely repulsive porous network, potentially corresponding to a homogeneous hydrophobic material. In contrast, MOFs are purely heterogeneous materials with diverse interactions with adsorbing molecules, from strongly repulsive to strongly attractive. Moreover, the authors agreed that the choice of sets of Monte Carlo moves significantly impacts the difference in performance between

single- and double-macrostate simulations. Furthermore, the object of comparison was the natural logarithm of macrostate probability distribution, which converges relatively slowly, while the resulting average number of adsorbed molecules reaches the equilibrium value more rapidly. This is a consequence of the fact that the MPD quickly reaches values close to equilibrium values, while exhibiting significant noise that slowly disappears with increasing number of MC moves (Figure 7 in cited paper⁴⁶). This noise, however, has little effect on the average number of molecules, as the average is strongly dependent on the location of the peak on the MPD profile because values from this region have largest weight, which follows directly from eq 7.)

## CONCLUSIONS

In this study, we presented and explored the practical application of the TMMC method in the NVT + ghost swap approach for computation of water adsorption isotherms in MOFs. This method uses the interpolation of transition probabilities to reduce computational cost while maintaining the accuracy of the results. Importantly, the TMMC method provides considerably more thermodynamic insight into adsorbing systems compared to GCMC simulations. In the case of water adsorption, it significantly reduces the calculation time and improves the precision due to overcoming the issues associated with simulation getting trapped in local energy minima. Moreover, for the first time, we combined interpolation of probabilities with temperature extrapolation method, obtaining high-quality results at temperatures beyond those directly simulated.

In the context of water harvesting, our study demonstrated that the working capacity of the sorbent, crucial for T/PSA processes, can be calculated with very high accuracy across an extensive range of pressures and temperatures using exclusively a single set of simulations conducted at single temperature and pressure. This allows for a rapid and effective material preselection for a given application and facilitate the optimization of the adsorption process itself. In particular, this methodology marks a significant step toward the advancement of large-scale computational screening of MOFs for water (and other gases) adsorption applications. Moreover, MPD data can be stored and easily recalculated for different thermodynamic conditions, aligning the described methodology with the current trend of working with reusable data.⁵² Furthermore, since MPD and free energy are directly related, these data can potentially be used for training machine learning algorithm for a more efficient search for new materials.⁵³,⁵⁴

To enhance the accessibility of the presented methodology, which holds strong applicative potential, we provide alongside this paper a modified version of the RASPA2 code⁴⁷ with a ghost swap move implementation. Additionally, we provide a Python library designed to minimize the user's input for analyzing data derived from the TMMC simulations. This library simplifies the interpolation of probabilities and extrapolation of MPD and allows for automatic calculation of isotherms (including writing to the AIF file⁵⁵ for data standardization purposes).

Future research could be dedicated to better understanding of the influence of all factors affecting the course of TMMC simulations of gas adsorption in MOFs. Furthermore, given the low acceptance rate of MC moves when simulating water adsorption and the formation of highly ordered water network
during adsorption, to further accelerate water adsorption
simulations advanced MC moves must be explored.

■ ASSOCIATED CONTENT

Data Availability Statement

Modified RASPA2 code for GC-TMMC simulations, python
code for TMMC data analysis with example Jupyter Notebook
and all exemplary simulation input files are available at https://
github.com/b-mazur.

⊕ Supporting Information

The Supporting Information is available free of charge at
https://pubs.acs.org/doi/10.1021/acsami.4c02616.

The following files are available free of charge. Details on
simulation methodology used in this work, additional
simulation results of adsorption isotherms, and tables
and figures referred in the main text (PDF)

■ AUTHOR INFORMATION

Corresponding Authors

**Bartosz Mazur** – Department of Micro, Nano, and Bioprocess Engineering, Faculty of Chemistry, Wroclaw University of Science and Technology, Wroclaw 50-370, Poland;
● orcid.org/0000-0002-6267-7135;

Email: bartosz.mazur@pwr.edu.pl

**Bogdan Kuchta** – Department of Micro, Nano, and Bioprocess Engineering, Faculty of Chemistry, Wroclaw University of Science and Technology, Wroclaw 50-370, Poland; MADIREL, CNRS, Aix-Marseille University, Marseille 13013, France; ● orcid.org/0000-0002-8635-4154;

Email: bogdan.kuchta@pwr.edu.pl

Author

**Lucyna Firlej** – Department of Micro, Nano, and Bioprocess Engineering, Faculty of Chemistry, Wroclaw University of Science and Technology, Wroclaw 50-370, Poland; Laboratoire Charles Coulomb (L2C), Université de Montpellier - CNRS, Montpellier 34095, France;

● orcid.org/0000-0002-8205-3522

Complete contact information is available at:
https://pubs.acs.org/10.1021/acsami.4c02616

Author Contributions

The manuscript was written through contributions of all
authors. All authors have given approval to the final version of
the manuscript.

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

B.M. is supported by the National Science Centre (NCN), Poland (grant no. 2022/45/N/ST4/03507). B.K. is supported by the National Science Centre (NCN), Poland (grant no. 2022/45/B/ST8/02028). We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2023/016797 which we used for DFT and WL/TMMC calculations. NVT + ghost swap calculations were carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (https://wcss.pl), grant no. 033.

■ REFERENCES

(1) Kalmutzki, M. J.; Diercks, C. S.; Yaghi, O. M. Metal-Organic Frameworks for Water Harvesting from Air. *Adv. Mater.* **2018**, *30* (37), 1704304.

(2) Hanikel, N.; Prévot, M. S.; Yaghi, O. M. MOF Water Harvesters. *Nat. Nanotechnol.* **2020**, *15* (5), 348–355.

(3) Zhou, X.; Lu, H.; Zhao, F.; Yu, G. Atmospheric Water Harvesting: A Review of Material and Structural Designs. ACS Mater. Lett. **2020**, *2* (7), 671–684.

(4) Lu, H.; Shi, W.; Guo, Y.; Guan, W.; Lei, C.; Yu, G. Materials Engineering for Atmospheric Water Harvesting: Progress and Perspectives. *Adv. Mater.* **2022**, *34* (12), 2110079.

(5) Hanikel, N.; Pei, X.; Chheda, S.; Lyu, H.; Jeong, W.; Sauer, J.; Gagliardi, L.; Yaghi, O. M. Evolution of Water Structures in Metal-Organic Frameworks for Improved Atmospheric Water Harvesting. *Science* **2021**, *374* (6566), 454–459.

(6) LaPotin, A.; Kim, H.; Rao, S. R.; Wang, E. N. Adsorption-Based Atmospheric Water Harvesting: Impact of Material and Component Properties on System-Level Performance. *Acc. Chem. Res.* **2019**, *52* (6), 1588–1597.

(7) Feng, A.; Akther, N.; Duan, X.; Peng, S.; Onggowarsito, C.; Mao, S.; Fu, Q.; Kolev, S. D. Recent Development of Atmospheric Water Harvesting Materials: A Review. ACS Mater. Au **2022**, *2* (5), 576–595.

(8) Hanikel, N.; Prévot, M. S.; Fathieh, F.; Kapustin, E. A.; Lyu, H.; Wang, H.; Diercks, N. J.; Glover, T. G.; Yaghi, O. M. Rapid Cycling and Exceptional Yield in a Metal-Organic Framework Water Harvester. ACS Cent. Sci. **2019**, *5* (10), 1699–1706.

(9) Hanikel, N.; Kurandina, D.; Chheda, S.; Zheng, Z.; Rong, Z.; Neumann, S. E.; Sauer, J.; Siepmann, J. I.; Gagliardi, L.; Yaghi, O. M. MOF Linker Extension Strategy for Enhanced Atmospheric Water Harvesting. ACS Cent. Sci. **2023**, *9* (3), 551–557.

(10) Ejeian, M.; Wang, R. Z. Adsorption-Based Atmospheric Water Harvesting. *Joule* **2021**, *5* (7), 1678–1703.

(11) Kim, H.; Rao, S. R.; Kapustin, E. A.; Zhao, L.; Yang, S.; Yaghi, O. M.; Wang, E. N. Adsorption-Based Atmospheric Water Harvesting Device for Arid Climates. *Nat. Commun.* **2018**, *9* (1), 1191.

(12) The CSD MOF collection contains just over 12,000 MOF crystal structures available to use for free in academic research. https://www.ccdc.cam.ac.uk/free-products/csd-mof-collection/ (accessed 2023-07-24).

(13) Chung, Y. G.; Haldoupis, E.; Bucior, B. J.; Haranczyk, M.; Lee, S.; Zhang, H.; Vogiatzis, K. D.; Milisavljevic, M.; Ling, S.; Camp, J. S.; et al. et al. Advances, Updates, and Analytics for the Computation-Ready, Experimental Metal–Organic Framework Database: CoRE MOF 2019. *J. Chem. Eng. Data* **2019**, *64* (12), 5985–5998.

(14) Chung, Y. G.; Gómez-Gualdrón, D. A.; Li, P.; Leperi, K. T.; Deria, P.; Zhang, H.; Vermeulen, N. A.; Stoddart, J. F.; You, F.; Hupp, J. T.; et al. In Silico Discovery of Metal-Organic Frameworks for Precombustion CO₂ Capture Using a Genetic Algorithm. *Sci. Adv.* **2016**, *2* (10), No. e1600909.

(15) Li, G.; Xiao, P.; Webley, P. A.; Zhang, J.; Singh, R. Competition of CO2/H2O in Adsorption Based CO2 Capture. *Energy Procedia.* **2009**, *1* (1), 1123–1130.

(16) Ho, C.-H.; Paesani, F. Elucidating the Competitive Adsorption of H₂ O and CO₂ in CALF-20: New Insights for Enhanced Carbon Capture Metal–Organic Frameworks. ACS Appl. Mater. Interfaces **2023**, *15* (41), 48287–48295.

(17) Erucar, I.; Keskin, S. Unlocking the Effect of H₂ O on CO₂ Separation Performance of Promising MOFs Using Atomically Detailed Simulations. *Ind. Eng. Chem. Res.* **2020**, *59* (7), 3141–3152.

(18) Caratelli, C.; Hajek, J.; Cirujano, F. G.; Waroquier, M.; Llabrés i Xamena, F. X.; Van Speybroeck, V. Nature of Active Sites on UiO-66 and Beneficial Influence of Water in the Catalysis of Fischer Esterification. *J. Catal.* **2017**, *352*, 401–414.

(19) Daglar, H.; Keskin, S. Recent Advances, Opportunities, and Challenges in High-Throughput Computational Screening of MOFs for Gas Separations. *Coord. Chem. Rev.* **2020**, *422*, 213470.
(20) Düren, T.; Bae, Y.-S.; Snurr, R. Q. Using Molecular Simulation to Characterise Metal–Organic Frameworks for Adsorption Applications. *Chem. Soc. Rev.* **2009**, *38* (5), 1237.

(21) Getman, R. B.; Bae, Y.-S.; Wilmer, C. E.; Snurr, R. Q. Review and Analysis of Molecular Simulations of Methane, Hydrogen, and Acetylene Storage in Metal–Organic Frameworks. *Chem. Rev.* **2012**, *112* (2), 703–723.

(22) Han, S. S.; Mendoza-Cortés, J. L.; Goddard III, W. A. Recent Advances on Simulation and Theory of Hydrogen Storage in Metal–Organic Frameworks and Covalent Organic Frameworks. *Chem. Soc. Rev.* **2009**, *38* (5), 1460.

(23) Paranthaman, S.; Coudert, F. X.; Fuchs, A. H. Water Adsorption in Hydrophobic MOF Channels. *Phys. Chem. Chem. Phys.* **2010**, *12* (28), 8123.

(24) Sarkisov, L.; Centineo, A.; Brandani, S. Molecular Simulation and Experiments of Water Adsorption in a High Surface Area Activated Carbon: Hysteresis, Scanning Curves and Spatial Organization of Water Clusters. *Carbon* **2017**, *118*, 127–138.

(25) Zhang, H.; Snurr, R. Q. Computational Study of Water Adsorption in the Hydrophobic Metal–Organic Framework ZIF-8: Adsorption Mechanism and Acceleration of the Simulations. *J. Phys. Chem. C* **2017**, *121* (43), 24000–24010.

(26) Choi, J.; Lin, L.-C.; Grossman, J. C. Role of Structural Defects in the Water Adsorption Properties of MOF-801. *J. Phys. Chem. C* **2018**, *122* (10), 5545–5552.

(27) Ghosh, P.; Kim, K. C.; Snurr, R. Q. Modeling Water and Ammonia Adsorption in Hydrophobic Metal–Organic Frameworks: Single Components and Mixtures. *J. Phys. Chem. C* **2014**, *118* (2), 1102–1110.

(28) Cirera, J.; Sung, J. C.; Howland, P. B.; Paesani, F. The Effects of Electronic Polarization on Water Adsorption in Metal-Organic Frameworks: H2O in MIL-53(Cr). *J. Chem. Phys.* **2012**, *137*(5)

(29) Datar, A.; Witman, M.; Lin, L.-C. Improving Computational Assessment of Porous Materials for Water Adsorption Applications via Flat Histogram Methods. *J. Phys. Chem. C* **2021**, *125* (7), 4253–4266.

(30) Datar, A.; Witman, M.; Lin, L. Monte Carlo Simulations for Water Adsorption in Porous Materials: Best Practices and New Insights. *AichE J.* **2021**, *67* (12), No. e17447.

(31) Witman, M.; Mahynski, N. A.; Smit, B. Flat-Histogram Monte Carlo as an Efficient Tool to Evaluate Adsorption Processes Involving Rigid and Deformable Molecules. *J. Chem. Theory Comput.* **2018**, *14* (12), 6149.

(32) Mazur, B.; Formalik, F.; Roztocki, K.; Bon, V.; Kaskel, S.; Neimark, A. V.; Firlej, L.; Kuchta, B. Quasicontinuous Cooperative Adsorption Mechanism in Crystalline Nanoporous Materials. *J. Phys. Chem. Lett.* **2022**, *13* (30), 6961–6965.

(33) Witman, M.; Wright, B.; Smit, B. Simulating Enhanced Methane Deliverable Capacity of Guest Responsive Pores in Intrinsically Flexible MOFs. *J. Phys. Chem. Lett.* **2019**, *10* (19), 5929–5934.

(34) Siderius, D. W.; Shen, V. K. Use of the Grand Canonical Transition-Matrix Monte Carlo Method to Model Gas Adsorption in Porous Materials. *J. Phys. Chem. C* **2013**, *117* (11), 5861–5872.

(35) Neimark, A. V.; Ravikovitch, P. I.; Vishnyakov, A. Adsorption Hysteresis in Nanopores. *Phys. Rev. E* **2000**, *62* (2), R1493–R1496.

(36) Ravikovitch, P. I.; Domhnaill, S. C. O.; Neimark, A. V.; Schueth, F.; Unger, K. K. Capillary Hysteresis in Nanopores: Theoretical and Experimental Studies of Nitrogen Adsorption on MCM-41. *Langmuir* **1995**, *11* (12), 4765–4772.

(37) Nakamura, M.; Ohba, T.; Branton, P.; Kanoh, H.; Kaneko, K. Equilibration-Time and Pore-Width Dependent Hysteresis of Water Adsorption Isotherm on Hydrophobic Microporous Carbons. *Carbon* **2010**, *48* (1), 305–308.

(38) Wang, Y.-M.; Datar, A.; Xu, Z.-X.; Lin, L.-C. *In Silico* Screening of Metal–Organic Frameworks for Water Harvesting. *J. Phys. Chem. C* **2024**, *128* (1), 384–395.

(39) Xu, Z.-X.; Wang, Y.-M.; Lin, L.-C. Connectivity Analysis of Adsorption Sites in Metal–Organic Frameworks for Facilitated Water Adsorption. *ACS Appl. Mater. Interfaces* **2023**, *15* (40), 47081–47093.

(40) Liu, X.; Wang, X.; Kapteijn, F. Water and Metal–Organic Frameworks: From Interaction toward Utilization. *Chem. Rev.* **2020**, *120* (16), 8303–8377.

(41) Shen, V. K.; Errington, J. R. Metastability and Instability in the Lennard-Jones Fluid Investigated by Transition-Matrix Monte Carlo. *J. Phys. Chem. B* **2004**, *108* (51), 19595–19606.

(42) Errington, J. R. Direct Calculation of Liquid–Vapor Phase Equilibria from Transition Matrix Monte Carlo Simulation. *J. Chem. Phys.* **2003**, *118* (22), 9915–9925.

(43) Siderius, D. W.; Hatch, H. W.; Errington, J. R.; Shen, V. K. Comments on “Monte Carlo Simulations for Water Adsorption in Porous Materials: Best Practices and New Insights. *AichE J* **2022**, *68* (8), No. e17686.

(44) Wang, F.; Landau, D. P. Determining the Density of States for Classical Statistical Models: A Random Walk Algorithm to Produce a Flat Histogram. *Phys. Rev. E* **2001**, *64* (5), 056101.

(45) Landau, D. P.; Tsai, S.-H.; Exler, M. A New Approach to Monte Carlo Simulations in Statistical Physics: Wang-Landau Sampling. *Am. J. Phys.* **2004**, *72* (10), 1294–1302.

(46) Hatch, H. W.; Siderius, D. W.; Errington, J. R.; Shen, V. K. Efficiency Comparison of Single- and Multiple-Macrostate Grand Canonical Ensemble Transition-Matrix Monte Carlo Simulations. *J. Phys. Chem. B* **2023**, *127* (13), 3041–3051.

(47) Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q. R. RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials. *Mol. Simul.* **2016**, *42* (2), 81–101.

(48) Fathieh, F.; Kalmutzki, M. J.; Kapustin, E. A.; Waller, P. J.; Yang, J.; Yaghi, O. M. Practical Water Production from Desert Air. *Sci. Adv.* **2018**, *4* (6), eaat3198.

(49) Mondloch, J. E.; Bury, W.; Fairen-Jimenez, D.; Kwon, S.; DeMarco, E. J.; Weston, M. H.; Sarjeant, A. A.; Nguyen, S. T.; Stair, P. C.; Snurr, R. Q.; et al. et al. Vapor-Phase Metalation by Atomic Layer Deposition in a Metal–Organic Framework. *J. Am. Chem. Soc.* **2013**, *135* (28), 10294–10297.

(50) Mahynski, N. A.; Blanco, M. A.; Errington, J. R.; Shen, V. K. Predicting Low-Temperature Free Energy Landscapes with Flat-Histogram Monte Carlo Methods. *J. Chem. Phys.* **2017**, *146*(7)

(51) Mahynski, N. A.; Errington, J. R.; Shen, V. K. Temperature Extrapolation of Multicomponent Grand Canonical Free Energy Landscapes. *J. Chem. Phys.* **2017**, *147* (5), 054105.

(52) Wilkinson, M. D.; Dumontier, M.; Aalbersberg, I. J.; Appleton, G.; Axton, M.; Baak, A.; Blomberg, N.; Boiten, J.-W.; da Silva Santos, L. B.; Bourne, P. E.; et al. et al. The FAIR Guiding Principles for Scientific Data Management and Stewardship. *Sci. Data* **2016**, *3* (1), 160018.

(53) Desgranges, C.; Delhommelle, J. Ensemble Learning of Partition Functions for the Prediction of Thermodynamic Properties of Adsorption in Metal–Organic and Covalent Organic Frameworks. *J. Phys. Chem. C* **2020**, *124* (3), 1907–1917.

(54) Desgranges, C.; Delhommelle, J. Towards a Machine Learned Thermodynamics: Exploration of Free Energy Landscapes in Molecular Fluids, Biological Systems and for Gas Storage and Separation in Metal–Organic Frameworks. *Mol. Syst. Des. Eng.* **2021**, *6* (1), 52–65.

(55) Evans, J. D.; Bon, V.; Senkovska, I.; Kaskel, S. A Universal Standard Archive File for Adsorption Data. *Langmuir* **2021**, *37* (14), 4222–4226.