Elucidating the particle size-dependent guest-induced structural transition of flexible metal–organic frameworks by exploring cooperative nature

Abstract

Flexible metal–organic frameworks (MOFs) exhibit a structural transition induced by adsorption of guest molecules. This guest-induced structural transition occurs at a certain gas pressure, resulting in an S-shaped adsorption isotherm. Consequently, these materials exhibit a high working capacity, making them highly competitive in energy-saving separation processes. However, the understanding of hysteresis loops between adsorption and desorption branches remains insufficient for industrial applications. Specifically, the particle size dependence of hysteresis behaviors is still actively being investigated. Generally, smaller particles of flexible MOFs show larger hysteresis loops. Herein, we constructed a simple multi-scale simulation model that couples molecular simulations for a unit cell with Ising lattice model-based simulations, in which solid–solid interactions for adjacent unit cells are considered, to address the cooperative nature within a flexible MOF particle. The solid–solid interactions strongly link unit cells in an identical state to form a domain, minimizing the heterointerface area. In transition states, the interfacial energy is independent of particle size, whereas the configurational entropy is significant for large particles, leading to a pronounced size dependence. This is applicable to real systems on the micron order, which is confirmed by the linear correlation between particle size and the free energy change of the unit cell over the hysteresis range. The correlation enables estimating particle size-dependent adsorption behavior, and consequently, tailoring the transition behaviors of flexible MOFs for target systems by controlling particle size. This study advances the understanding of hysteresis in guest-induced structural transitions and provides insights for designing adsorption-based separation processes.

---

1 Introduction

Recently, the disparity between the actions taken and those necessary to reduce greenhouse gas emissions has widened.¹ Consequently, the development of energy-efficient separation technologies, including adsorption and membranes, is an urgent priority.² Metal–organic frameworks (MOFs) have attracted significant attention for adsorption-based separation processes owing to their inherent potential arising from the infinite combinations of metal ions and coordinated organic linkers, which allows control over pore size and host–guest interactions.³ Several MOFs exhibit structural flexibility in response to external stimuli,^4–8 suddenly changing their structures from a closed-pore (cp) to an open-pore (op) form at a specific threshold pressure during adsorption. This guest-induced structural transition is called the “gate opening” behavior and adsorption isotherms on flexible MOFs have a characteristic S-shape. Other specific MOFs show a reversible transition behavior through a single adsorption process, occurring from a large-pore (lp) to a narrow-pore (np) form at a low pressure and subsequently from the np to lp form at a high pressure. This reversible transition, so-called the “breathing” behavior, exhibits the similar adsorption isotherms and occurs due to the same mechanism as the gate adsorption. A stepwise isotherm enables a larger working capacity with only a slight pressure change, compared to the conventional type I isotherm. Guest recognition^4,9 and intrinsic thermal management^10,11 capabilities also have positive effects on adsorption-based separation processes,¹² rendering flexible MOFs promising materials that should be urgently investigated.¹³

In general, adsorption behavior involving a phase transition is accompanied by hysteresis. A typical example is capillary condensation in nanopores, in which the evaporation (desorption) process starts at a lower pressure (P_des) than the condensation pressure during adsorption (P_ads).^14,15 The hysteresis behavior is not preferable for adsorption-based separation processes, because compared to an ideal adsorbent that exhibits an equilibrium transition at the same pressure for both adsorption and desorption branches, the adsorption and desorption processes require increasing and decreasing pressures to induce the transition, respectively, resulting in additional energy consumption.¹⁶ However, the capillary condensation/desorption pressures are determined only by the pore size and the surface properties of the pore walls,¹⁷ indicating that the width of the hysteresis loop cannot be controlled within the same material. Similarly, the guest-induced structural transition of flexible MOFs involves hysteresis and depends on the type of MOF;¹⁸ however, a crucial difference from the capillary condensation is that the hysteresis behavior depends on the particle size. Sakata et al. first reported that smaller particles of [Cu₂(bdc)₂(bpy)] (bdc = 1,4-benzenedicarboxylate, bpy = 4,4′-bipyridine) exhibit larger hysteresis.¹⁹ This tendency has also been observed in other flexible MOFs, such as breathing,²⁰ swelling,²¹ and linker-bending types.^22–26 Interestingly, an anisotropic size effect was observed for DUT-8(Ni) ([Ni₂(ndc)₂(dabco)]; ndc = 2,6-naphthalenedicarboxylate, dabco = 1,4-dizabicyclo[2.2.2]octane).^24,27 The crystal width, which corresponds to the planar direction of the sheet composed of ndc paddle-wheels, affected its hysteresis loop, whereas the crystal length, aligning with the direction in which dabco ligands link the ndc sheets, had no significant effect. These observations indicate that elucidating the size effect can uncover new possibilities for controlling the hysteresis behavior while maintaining the other desired properties of the target MOFs.

The theoretical understanding of hysteresis loops has been developed based on a thermodynamic framework for guest-induced structural transitions,²⁸ which explains the transition pressure as the point at which the free energy of the cp structure equals that of the op structure encapsulating the guest molecules. Based on molecular simulations using a toy model that mimics jungle gym-type MOF structures, we demonstrated the existence of an energy barrier between the cp and op phases,²⁹ and many subsequent corroborating studies have reported that the hysteresis loop is caused by the additional pressure variance required for the system to overcome the energy barrier.^30–33 However, the model cannot explain the size effect of the hysteresis behavior (which is discussed in detail below), mainly because it implicitly assumes a uniform transition within a bulk crystal, represented by a series of replicates of a nanoscale periodic unit cell. The current understanding of the energy barrier fails to incorporate aspects that are considered important when considering size dependence. These aspects include the effects of surface properties and the cooperative nature observed over long distances, which are only apparent in mesoscale systems wherein constituting unit cells are allowed to adopt different states. Several pioneering studies have been reported regarding this cooperative nature.^34–39 Vandenhaute et al.³⁷ developed a mesoscale simulation box of MIL-53(Al), showing that the simulation box exhibits a mixed-phase state during the transition induced by mechanical pressure, in which a unit cell within the simulation box tends to take on the same phase as its neighbors, thereby forming domain structures of cp and op phases. Mitsumoto and Takae³⁸ constructed a toy model that accounts for distortion energy upon structural deformation, revealing that the difference in the domain structures between adsorption and desorption processes causes hysteresis. Their results suggest that the transition state in actual systems does not lie in the intermediate structure between the cp and op structures, as suggested by earlier studies, but rather in the cooperative nature where the interactions between neighboring flexible motifs play a key role.^40,41 Therefore, a detailed investigation into the cooperative nature to redefine the “true” energy barrier would reveal the size-dependent hysteresis behavior.

In this study, we elucidated the guest-induced structural transition including the particle size dependence of the hysteresis loop. We developed a simple multi-scale simulation model that couples molecular simulations for a nanoscale unit cell with Ising lattice model-based simulations for a mesoscale system, thereby addressing the cooperative nature between motifs within frameworks of a MOF particle. As a model case, we focused on ELM-11 ([Cu(BF₄)₂(bpy)₂]),^42,43 which possesses a stack-layered structure and exhibits a guest-induced structural transition arising from the expansion of its interlayer widths. The material displays a high degree of selectivity for CO₂/CH₄ separation, thereby realizing an energy-efficient process for the treatment of industrial exhaust gas.¹²

2 Theoretical procedures

2.1 Thermodynamics of guest-induced structural transition

The osmotic free energy of a unit cell, as part of a bulk MOF crystal having a stack-layered structure, denoted as ω_os, is represented by

ω_os(h,P) = f_host(h) + ω_guest(h,P) (1)

where h is the interlayer width, P is the pressure, and f_host and ω_guest are the free energies associated with host-framework deformation and guest adsorption, respectively. We define h_cp as the interlayer width at which ω_os reaches its minimum value when P = 0, corresponding to the cp phase. The free energy change accompanying layer expansion from h_cp to h, Δω_os, is expressed aswhere Δf_host and Δω_guest are changes in the f_host and ω_guest values from h = h_cp, respectively. The equilibrium phase transition between the cp (h = h_cp) and op (h = h_op) states can be explained as occurring at P = P_eq that satisfies Δω_os(h_op,P_eq) = 0. Most theoretical studies on the guest-induced structural transition have relied on this thermodynamic framework²⁸ (hereafter designated as the primitive theory). However, discussions based on the periodic unit cell implicitly assume a bulk crystal in which all unit cells are in the same state. Thus, when explicitly considering a MOF particle consisting of M unit cells, the primitive theory treats the entire free energy change of the particle, ΔΩ_os, as:

ΔΩ_os = MΔω_os(h,P) (3)

However, it is better to express ΔΩ_os as a function of interlayer widths composed of unit cells, h = (h₁, h₂, …, h_M). ΔΩ_os(h,P) should include the free energy term for each unit cell, Δω_os(h_i,P), and interactions between unit cells representing a penalty for adjacent cells in different states. In this context, the primitive theory corresponds to the assumption of an infinite inter-unit-cell penalty, which enforces all unit cells to have an identical interlayer width. Note that in the derivation of the equation for the structural transition-type adsorption (STA equation³⁹) reported by our group, ΔΩ_os(h,P) is modeled as follows:

where s represents the average number of unit cells that simultaneously deform as one domain. Namely, eqn (4) is a simplified expression of the cooperative nature; the essence of inter-unit-cell interaction is abstracted as an effective domain size within which the unit cells undergo cooperative deformation.

In contrast, the present study directly addresses the inter-unit-cell interaction. Specifically, ΔΩ_os(h,P) was modeled as

where F_IUC is the free energy associated with the inter-unit-cell interaction. As an initial approach, this study introduced several assumptions. First, the transition states of the unit cells were disregarded, allowing the unit cells to be categorized exclusively into either cp or op states. This assumption converts the first term on the right-hand side of eqn (5) intowhere m is the number of op cells in the particle, and Δω_os(h_op,P) is simply denoted as Δω^op_os(P). Second, inter-unit-cell interactions were limited to only between adjacent cells. Under these two assumptions, the inter-unit-cell interaction energy, E_IUC, for a given configuration can be written aswhere r^(m) is the index vector that indicates the positions of op cells when the number of op cells is m. n_i is the number of interfaces at which different states are adjacent along the i-direction, and J_i (>0) is the interfacial energy, serving as a penalty for adjacent different states. Thus, F_IUC for a given m is expressed as

Through the last transformation in eqn (8), the meaning of r changes from the index vector to the integer coordinates for m distinguishable op “particles.” Consequently, the summation is divided by m! to eliminate duplication. Furthermore, E_IUC is also extended to return ∞ if any two of the m particles have identical coordinates, ensuring that exp(−E_IUC/k_BT) = 0. Taken together, ΔΩ_os can be modeled as a function of m and P as

where r^m = (r₁, r₂, …, r_m) for simplicity. Eqn (9) presents the standard form of the equation that expresses the relationship between energy and free energy, implying that if E(r^m,P) = mΔω^op_os(P) + E_IUC(r^m) can be evaluated, Monte Carlo (MC) techniques are applicable for analysis. Specifically, by using canonical MC simulations, in which the transition to a proposed configuration r_j from r_i is accepted with a probability:

we can obtain the ensemble average of m, 〈m〉, at a given P. Consequently, the average amount adsorbed, N_guest, can be determined using the adsorbed amount for a unit cell in the op phase, n^op_guest, such that:

N_guest(P) = 〈m〉n^op_guest(P) (11)

Furthermore, the test particle method⁴⁴ enables us to evaluate the profile of ΔΩ_os(m,P) as

where E⁺ is the energy received when an additional op cell is inserted at random into a system with a configuration of r^m (see Section S1, ESI† for details on the derivation of eqn (12)).

2.2 Simulation model

To evaluate E(r^m,P), we performed two types of simulations. The first, named “nanoscale simulation,” involves conducting molecular simulations for a toy model that mimics the framework structure of ELM-11 to obtain variables regarding the unit cell, such as n^op_guest(P) and Δω^op_os(P) (Fig. 1a). We used the same toy model as in our previous studies,^31–33 where the stack-layered structure of ELM-11 is represented by a pair of planes with a pillar atom and the cell size of 10σ_gg × 10σ_gg × h (σ_gg: the size parameter of the Lennard-Jones (LJ) potential for a guest molecule). The second, named “mesoscale simulation,” represents the cooperative behavior within a particle. We constructed a particle model consisting of M_x × M_y × M_z (= M) unit cells that adopt either the cp or op phase, with interfacial energy J_i imposed if the adjacent cell in the i-direction is in the different state. MC simulations were conducted for this particle model to obtain N_guest(P) and ΔΩ_os(m,P) according to eqn (11) and (12) (Fig. 1b). As the stack-layered structure should have a much stronger connection in the x- and y-directions due to coordination bonds than in the z-direction which stems from van der Waals interactions, we assumed J_x = J_y (≡ J_#), J_z = 0, and M_x = M_y (≡ M_#) and varied J_#, M_#, and M_z as a first step. For more details, refer to Section 5.

Fig. 1 Schematics of (a) the toy model that mimics the framework structure of ELM-11 and (b) the particle model where the unit cells can undergo transitions while interacting with adjacent cells.

Note that all variables are discussed in their dimensionless forms using the LJ parameters of the guest molecule, σ_gg and ε_gg, where ε_gg is the depth of the LJ potential. Throughout all simulations, the temperature was fixed at k_BT/ε_gg = 1.

3 Results and discussion

3.1 Particle size-dependent behaviors of ELM-11

First, we experimentally examined the size-dependent behaviors of ELM-11 (see Section 5 for synthesis and analysis procedures). Fig. 2a–e show scanning electron microscopy (SEM) images of ELM-11 particles synthesized under various concentration conditions. X-ray diffraction (XRD) measurements confirmed that these particles are the hydrated form of ELM-11 (ref. 45) (Fig. 2f), and the particles sizes on the longer side measured 22.6 ± 10.0, 10.2 ± 4.1, 7.8 ± 2.6, 7.6 ± 2.7, and 6.7 ± 2.6 μm, respectively (particle size distributions are provided in Section S2, ESI†). Higher concentrations resulted in smaller particles, which agreed with the classical nucleation theory and the results of ZIF-8 synthesis.^46,47Fig. 2g shows CO₂ adsorption isotherms on the ELM-11 samples at 273 K, indicating a clear trend wherein smaller particles exhibit larger hysteresis and less distinct S-shaped curves. This aligns with the tendencies exhibited by other flexible MOFs.^{19,20,22–26,48}

Fig. 2 (a)–(e) SEM images, (f) XRD patterns, and (g) CO₂ adsorption isotherms (273 K) of ELM-11 particles synthesized under various concentration conditions. The closed and open symbols in (g) represent adsorption and desorption branches, respectively.

3.2 Nanoscale simulation of the unit cell

Fig. 3a depicts the free energy landscape of the unit cell along the interlayer width. At P = 0, the free energy increases as h increases, which represents the free energy change of host-framework deformation, Δf_host, as no guest molecules are adsorbed. With increasing pressure, guest adsorption stabilizes the system, leading to a decrease in Δω_os for a larger h. The bistable state is established at P = P₃, indicating that P₃ is the equilibrium transition pressure for a unit cell (P_eq). As the stable state switches across P_eq, the adsorption isotherm exhibits a stepwise shape, as depicted by the dashed line in Fig. 3b. From these results, we determined h_op = 2.00σ_gg and obtained n^op_guest(P) and Δω^op_os(P) profiles from the data, as shown in Fig. 3c and d, which were used in the subsequent mesoscale simulations.

Fig. 3 (a) Free energy profiles of the unit cell at various pressures. Filled star symbols denote the stable state at each pressure, while the open star symbol denotes a metastable state. (b) Adsorption isotherms of the system. The equilibrium transition from cp to op states at P = 0.00039 (P₃) illustrates a step-wise adsorption behavior (dashed line), and the kinetic transition occurs with a hysteresis loop when the height of energy barriers is lower than 2k_BT. (c) The resulting n^op_guest(P) and (d) Δω^op_os(P) profiles, which are used in the mesoscale simulations.

Note that the primitive theory attributes the cause of the hysteresis loop to the energy barriers found in the free energy profiles shown in Fig. 3a.³¹ For example, assuming that the system can overcome an energy barrier lower than 2k_BT through thermal fluctuation, the transition during adsorption would occur at P = P₅ (> P_eq), while during desorption, it would occur at P = P₂ (< P_eq), as depicted by the solid lines in Fig. 3b. However, this approach has resulted in an incomplete understanding in two aspects: first, the simulated adsorption isotherms exhibit a stepwise shape at the transition pressure, whereas experimental results show gradual changes within narrow pressure ranges (Fig. 2g). Second, according to eqn (3), ΔΩ_os is the product of Δω_os and the number of unit cells. This suggests that a larger particle has a higher energy barrier, and thus, a larger hysteresis loop (details are present in Section S3†), which contradicts the experimental observations (Fig. 2g). To explain the experimental results, the energy barrier must be lower for larger particles.

3.3 Mesoscale simulation for the particle model

Fig. 4a depicts adsorption isotherms for the particle model (M_# = M_z = 10) with various J_#. These adsorption behaviors showed three features: (1) weak interaction between unit cells, as represented by J_# = 0.6k_BT, resulted in no hysteresis loops. (2) Moderate interaction caused a step-by-step change in the amount adsorbed, accompanied by hysteresis. (3) Strong interaction expanded the hysteresis loop. All these points can be explained by the free energy landscapes at P = P_eq, i.e., ΔΩ_os(m,P_eq). Fig. 4b shows free energy profiles with various J_# values at P = P_eq. A small J_# exhibited a downward-convex profile with no energy barriers. In contrast, the profiles with a larger J_# exhibited repeated upward convexes; these energy barriers inhibited the equilibrium transition and caused hysteresis. Further increasing J_# increased the energy barriers and consequently widened the hysteresis loop.

Fig. 4 (a) Adsorption isotherms for the particle model and (b) free energy profiles with various J_#. (c) and (d) Display snapshots of the particle model at the points marked with star symbols in (b), where pink and blue cubes represent op and cp cells, respectively. The layers are displayed with a slight offset.

In principle, free energy is determined by the balance between internal energy and entropy. Because Δω^op_os equals zero at P = P_eq, the free energy profiles illustrated in Fig. 4b consist solely of F_IUC, which is determined by the balance between how comfortable a cell feels towards its neighbors and how freely op cells are placed. The downward-convex profiles with smaller J_# demonstrate that, under these conditions, the entropic contribution controls the system as the interfacial energy is too small. A snapshot at the stable point, highlighted with a star symbol, illustrates the seemingly random coexistence of cp and op cells, supporting this explanation (Fig. 4c). The free energy profile remained downward convex with increasing pressure while the value of m at which ΔΩ_os reaches its minimum varied, resulting in a gradual increase in adsorption without hysteresis (details are provided in Section S4†). In contrast, as shown in Fig. 4d, increasing J_# tends to make all unit cells within the same layer share an identical state so that the interfacial contribution is minimized, leading to M_z + 1 local minima when m is a multiple of M_#². Therefore, the system needs to overcome an energy barrier M_z times, resulting in a step-by-step adsorption increase with a hysteresis loop. These findings indicate that an experimentally obtained S-shaped isotherm with a large increase in uptake within a narrow pressure range originates from the layer-by-layer transition in a MOF particle.

3.4 Size dependence of a guest-induced structural transition

Fig. 5a shows adsorption isotherms for a particle with various M_# under the conditions of M_z = 10 and J_# = 2.2k_BT, demonstrating that decreasing M_# resulted in a wider hysteresis loop. Because M_# can be regarded as being proportional to particle size, the results qualitatively coincide with the experimental observation (Fig. 2g). Fig. 5b shows that the energy barrier at a pressure higher than P_eq increases as M_# decreases, which should be the cause of the widening hysteresis loop. To clarify the mechanism behind the obtained results, we attempt to simplify ΔΩ_os(m,P) (eqn (9)) under appropriate assumptions. As discussed in the previous section, the process that causes an energy barrier fundamentally involves a transition within a layer. Because J_z was set to 0 in the simulations, each layer can be regarded as independent from the other layers. This means that one simulation with M_z layers is statistically identical to the average of M_z individual simulations with a single layer (details are provided in Section S5†). Therefore, the key is to understand the free energy profile of a single layer, as shown in Fig. 5c where M_# = 1 and J_# = 2.2k_BT. The snapshots in Fig. 5c show that, with such a large J_#, op cells exist adjacent to each other to form a domain, which allows categorizing the transition into three stages: both the initial and final stages feature square-shaped domains, albeit with different phase compositions, while the intermediate stage presents a rectangular shape with the long side of M_#. Both square and rectangular shapes have the minimum interfacial areas under periodic boundary conditions. Given that the intermediate stage presents a plateau region in ΔΩ_os(m,P_eq) and that ΔΩ_os(m,P) can be expressed as mΔω^op_os(P) + ΔΩ_os(m,P_eq), the transition state having a local maximum energy emerges in the initial stage during the adsorption process (Δω^op_os < 0), as demonstrated in Fig. 5b. Here, ΔΩ_os(m,P) can be rewritten in a different form from eqn (9) using the energy levels of the interfacial interaction, ε_IUC,c (c: the index for the energy level), and the number of configurations satisfying E_IUC = ε_IUC,c, denoted as w_c, as follows:

Fig. 5 (a) Adsorption isotherms for the particle model with various M_#. (b) Free energy profiles at a higher pressure than P_eq, indicating that a smaller M_# exhibits a higher energy barrier. (c) The free energy profile of a single layer (M_z = 1) with M_# = 10 and J_# = 2.2k_BT at P = P_eq, showing three stages with different domain shapes. The red line represents an approximated curve of the free energy profile during the initial stage (eqn (14)). (d) Schematic representations of the states at m = 16.

Note that in eqn (13), the op cells are regarded as indistinguishable. Assuming a square number for m, the energy level in the initial stage should be in ascending order, as illustrated in Fig. 5d. Under the conditions where the interfacial energy is governed, considering only the lowest energy level state should provide a good approximation. Therefore, eqn (13) can be transformed as

The red line in Fig. 5c represents eqn (14), demonstrating a reasonable approximation for the profile in the initial stage, including non-square numbers of m. Differentiating eqn (14) with respect to m provides the transition state as

Therefore, the height of the energy barrier, ΔΩ_tr, can be approximated as:

While the first energetic term remains constant regardless of layer size, the second entropic term reduces the energy barrier more significantly for larger layers. Therefore, we conclude that particle size-dependent hysteresis primarily arises from the entropic contribution.

Eqn (16) also indicates an important relationship between particle size and transition pressure. Assuming that the system can overcome the energy barrier with the height of κ through its thermal fluctuations, eqn (16) can be rewritten as

where we explicitly express the size dependence of the transition pressure. As κ seems to be independent of M_#, eqn (17) represents the linear relationship between ln M_# and 1/Δω^op_os. Fig. 6a shows the plot of ln M_#vs. 1/Δω^op_os obtained from simulations for single layers with M_# = 10, 15, 30, 50, and 100, where Δω^op_os was evaluated at the pressure at which the adsorption branch reaches its transition ratio m/M_#² of 0.5 (see Fig. S6† for the obtained isotherms). Fig. 6a demonstrates clear linearity between ln M_# and 1/Δω^op_os, which agrees well with eqn (17). However, the slope and intercept of the regression line indicate that J_# = 1.49k_BT and κ = 10.4k_BT, which differ slightly from the exact values: J_# = 2.2k_BT (a set value) and κ = 23k_BT (the average energy barrier heights in the actual free energy profiles for layers with various M_#). These deviations probably arise from the assumption that only the lowest energy state is considered. Nonetheless, this linearity is useful to validate the applicability of our simulation to experimental systems, as discussed below.

Fig. 6 Plot of ln M_#vs. 1/Δω^op_os for (a) the simulation results and (b) the experimental results (Fig. 2g). The dashed lines represent the regression lines.

3.5 Validation of the simulation model through comparison with experimental results

Due to computational costs, our simulation was limited to a maximum size of M_# = 100. Considering that the unit cell size of ELM-11 is approximately 1.1 nm (the length of Cu–bpy–Cu), our results only cover less than 100 nm. Therefore, we examined the linear relationship between ln M_# and 1/Δω^op_os in experimental results (Fig. 2g) to indirectly demonstrate the applicability of our simulation model to real-system conditions, specifically at the micrometer scale. According to the primitive theory, Δω^op_os(P_ads) can be obtained by²⁸where V_m is the molar volume of the guest. P_eq is impossible to observe in real systems due to the hysteresis loop. Our previous studies proposed that P_eq is closer to P_des than P_ads because pressure responsiveness to reducing the energy barrier for a nanoscale transition is different (Fig. 3a and b); the barrier height decreases more drastically around P_des than P_ads.^29,31 However, the present study indicates that the actual energy barrier controlling the kinetic behavior of the system is based on the cooperative nature at the mesoscale, in which the free energy profile at P = P_eq displays symmetry (Fig. 5c). The symmetry suggests that the transition state during the desorption process can be identified in the final stage in Fig. 5c, and that the height of the energy barrier during the desorption process can be modeled in a manner analogous to that described in eqn (14)–(16). Considering that the height of the energy barrier at P_ads and P_des should be identical, the following must be satisfied:

Δω^op_os(P_ads) = −Δω^op_os(P_des) (19)

Therefore, eqn (18) can be transformed without P_eq as

In the second transformation, n^op_guest was modeled by the Langmuir equation, characterized by the parameters n_∞ and K. Note that eqn (19) also provides the important insight that P_eq can be numerically determined from an experimental isotherm (see Section S7 for details†). As n_∞ and K of the ELM-11–CO₂ system were determined to be 3.910 mmol g⁻¹ and 0.201 kPa⁻¹ at 273 K, respectively,³⁹ Δω^op_os can be evaluated from the experimental adsorption and desorption branches. Fig. 6b shows the plot of ln M_#vs. 1/Δω^op_os for ELM-11 derived from experimental results, where the M_# values were calculated by dividing the mean particle size by 1.1 nm. The experiments also confirmed clear linearity, suggesting our transition model is valid in actual MOFs. Substituting eqn (20) into eqn (17) yields:

By extrapolating the regression line in Fig. 6b, we can estimate how the width of the hysteresis loop changes with particle sizes. For example, if we could synthesize one order of magnitude larger ELM-11 particles (100 μm), the width of the hysteresis loop would be 4.8 kPa, which is 1.7 times narrower than that of 22.6 μm ELM-11 particles (Fig. 2g). This would provide guidelines for synthesizing a flexible MOF with a narrow hysteresis loop, thereby reducing the additional pumping cost in adsorption-based separation processes. However, eqn (21) also indicates that the hysteresis loop never disappears even if millimeter-sized single crystals are synthesized (e.g., 3.2 kPa for 1 mm and 2.4 kPa for 10 mm), because the size responsiveness becomes milder for larger particles, which follows an experimental trend.²⁴

Finally, we discuss the validity of the approach taken in this study, which posits that the transition state controlling the system's behavior stems from the interfacial energy rather than the intermediate structure during the layer expansion of the unit cell. At P = P_eq where Δω^op_os = 0, the transition state of the entire particle is at the intermediate stage of Fig. 5c. Considering only the lowest energy state, the energy barrier at the intermediate stage can be approximated as

ΔΩ_tr(P_eq,M_#) = 2M_#J_# − k_BT ln 2M_# (22)

To disregard the effect of the intermediate structure during the layer expansion, the following points should be verified: first, Δω_tr is much smaller than ΔΩ_tr, thereby allowing the height of the energy barrier to be regarded as ΔΩ_tr. Second, the total energy barrier when all the unit cells undergo transition simultaneously, i.e., M_#²Δω_tr, is much greater than the ΔΩ_tr value, thereby making the transition state stemming from the interfacial energy an actual saddle point. Thus,

Δω_tr ≪ 2M_#J_# − k_BT ln 2M_# ≪ M_#²Δω_tr (23)

is the required condition. By assuming Δω_tr ≃ 4.6k_BT (Fig. 3c) and 10 μm particles (M_# ∼ 10⁴), the order of J_# should be much larger than 10⁻³k_BT and smaller than 10⁴k_BT. It is plausible to assume that the energy arising from the interfaces is larger than the thermal motion (∼k_BT) but smaller than the reaction heat (∼10²k_BT) involved in cleaving and connecting chemical bonds,⁴⁹ which satisfies the requirement condition. Therefore, we can conclude that the intermediate structure during the layer expansion does not severely affect the system's transition behavior nor hysteresis behaviors.

4 Conclusion

This study explored the particle size dependence of guest-induced structural transitions from experimental and theoretical perspectives. We conducted mesoscale simulations of multiple unit cells, considering the interfacial interactions between adjacent cells in different states, to understand the size effect. These interactions strongly bind cells in the same state, creating a domain that reduces the heterointerface area. While the energetic term contributing to the free energy profile is independent of particle size, the configurational entropic term increases with particle size, resulting in a pronounced size dependence of the hysteresis loop.

Finally, we note three limitations of this study and one perspective. The limitations are as follows:

(1) In our model, we consider interactions only between adjacent cells, which might be insufficient for modeling guest-induced structural transitions that involve drastic deformation. To accurately model the distortions, considering interactions with non-adjacent or diagonally adjacent cells would be effective, as some studies have demonstrated in exploring solid–solid phase transitions in complex magnetic materials.^50,51

(2) Eqn (17) has room for improvement. Although the experimental results showed a linear correlation between ln M_# and 1/Δω^op_os (Fig. 6b), comparing the regression line to eqn (17) indicates that κ is less than zero, i.e., a negative thermal fluctuation. This unnatural outcome may stem from the accuracy of eqn (14) (as discussed above), the discrepancy between an actual ELM-11 and the ideal condition of J_z = 0, and the way that M_# was derived from experiments (where particle size was considered in Fig. 6b, though crystallite size might be more appropriate). An improved linear equation that addresses these concerns is highly promising, offering the potential to gain quantitative insights into interfacial energy and thermal fluctuation in real systems.

(3) Although our model demonstrates the typical size dependence shown by specific flexible MOFs, such as ELM-11, Cu₂(bdc)₂(bpy),¹⁹ MIL-53,²⁰ and DUT-8,^22–26 some other unique size dependence cannot be explained. Specifically, ZIF-8 exhibits a shift to a higher transition pressure in both the adsorption and desorption branches as the particle size decreases. A previous study suggested that this phenomenon is caused by the weaker adsorption potential near the surface; therefore, a smaller particle with a larger specific surface area requires additional pressurization to undergo a transition in both the adsorption and desorption branches.⁵² This exception highlights the necessity to expand our model to incorporate surface contributions for a comprehensive understanding of the size-dependent hysteresis behavior.

Despite these limitations, our simulation model is useful not only for understanding the fundamental mechanism behind size-dependent hysteresis behavior, but also for exploring dependences on framework type. In this study, we focused on stack-layered MOFs having stronger connections within their two-dimensional layers. However, some MOFs display stronger connections along one-dimensional paths, such as [Cu(BF₄)₂(bpp)₂] (bpp = 1,3-bis(4-pyridyl)propane).^53,54 The behavior of these one-dimensional chain MOFs could be simulated using our model with J_z > 0 and J_# = 0, preliminary results of which are shown in Section S8.† Although a trend similar to that observed in the stack-layered MOF was obtained with an increase in interfacial energy, the S-shaped uptake upon structural transition was less distinct compared to those in Fig. 4a. This observation aligns with the experimental results: [Cu(BF₄)₂(bpp)₂] exhibits a more gradual adsorption increase compared to ELM-11.^53,54 Furthermore, when J_# and J_z have specific non-zero values, a decrease in M_# resulted in a widening of the hysteresis loop, whereas M_z had no impact on the behavior. This tendency agrees with the anisotropic size effect observed in DUT-8(Ni)²⁴ (see Section S8† for prospective results). Additionally, our simulation model also has the potential to elucidate a more generalized system with three-dimensional interfacial interactions (Fig. S9†). A systematic investigation using our model, considering the type of connectivity in the framework, would provide a comprehensive understanding of the size-dependent hysteresis behavior of flexible MOFs; such a study is currently underway and will be reported elsewhere.

5 Methods

5.1 Synthesis of ELM-11 samples

ELM-11 was synthesized by the following procedure, which is based on our previous report.^45,55 First, 0.295 g 4,4′-bipyridine (bpy; 98.0%, Tokyo Chemical Industry Co., Ltd) was dissolved in 1.181 mL methanol (99.%, Kishida Chemical Co., Ltd). Second, 0.283 mL of a Cu(BF₄)₂ 45% aqueous solution (Kanto Chemical Co., Inc.) was mixed with 0.738 mL deionized water prepared using an Arium^®mini plus (Sartorius) water purification system. Subsequently, the bpy methanol solution was added dropwise into a Cu(BF₄)₂ aqueous solution in a vial while stirring with a magnetic stirrer at room temperature (∼297 K). The dropwise addition was completed slowly over 2 h, followed by mixing for several days. The synthesized particles were collected by vacuum filtration, washed several times with deionized water, and then dried under vacuum overnight. The particle size of ELM-11 becomes smaller as the concentration of bpy and Cu(BF₄)₂ is higher. Therefore, we varied the concentrations with a fixed ratio to control the size range of the resulting particles ([bpy], [Cu(BF₄)₂]) = (0.8 M, 1.6 M), (1.0 M, 2.0 M), (1.2 M, 2.4 M), (1.4 M, 2.8 M), and (1.6 M, 3.2 M); the particles obtained under these conditions correspond to Fig. 2a–e, respectively.

5.2 Characterization

XRD patterns were obtained using an UltimaIV/285/DX (Rigaku Corp., Japan) with Cu Kα radiation (40 kV and 20 mA) at a scan rate of 10° min⁻¹ and step size of 0.01°. SEM images were obtained using an SU8220 field emission SEM (Hitachi High-Tech Corp., Japan). The mean particle size was calculated by measuring the length of at least 70 particles using the SEM images. CO₂ adsorption isotherms were measured at 273 K using a BELSORP-mini (MicrotracBEL Corp., Japan), in which the equilibrium setting was pressure variation within 30 Pa for 600 s and the sample amount was approximately 50 mg. The samples were dehydrated and activated by degassing at 393 K for 10 h under vacuum below 10 Pa before adsorption measurements.

5.3 Nanoscale simulations for the unit cell

Grand canonical MC simulations were conducted on 36 models with slightly different interlayer widths from h = 1.70σ_gg to h = 2.05σ_gg at various pressures (Fig. 1a), obtaining the amount adsorbed on the unit cell, n_guest(h,P). By computing relative free energy of each structure without guest molecules and the grand potential of the guest molecules, ω_os(h,P) was obtained. Details of the simulation are present elsewhere³¹ or in Section S9.†

5.4 Mesoscale simulations for the particle model

We performed two types of mesoscale simulations using the particle model, in which periodic boundary conditions were imposed in all directions.

The first was to obtain the adsorption isotherms, which resembles a standard simulation method for the Ising lattice model. The simulation procedures were as follows:

(1) Set a pressure and obtain the Δω^op_os(P) value from the nanoscale simulation.

(2) Choose one cell at random and convert its state into the other one.

(3) Accept the conversion with the probability defined by eqn (10).

(4) Perform Steps 2 and 3 M × 10, 000 times and calculate the ensemble average number of op cells, 〈m(P)〉, in which the first half of the trials are used for equilibration, and the sampling is conducted in the latter half.

(5) Repeat Steps 1–4 with increasing/decreasing pressure to obtain adsorption/desorption isotherms.

Except for the results shown in Fig. 4a, the adsorption and desorption isotherms were averaged over 10 separate simulations to obtain smooth profiles.

The second method mimicked Widom's test particle technique to calculate free energy profiles. The simulation procedures are outlined as follows:

(1) Set the number of op cells, m, in the particle model and the pressure to P = P_eq, where Δω^op_os(P_eq) = 0.

(2) Choose one of op cells at random and convert it into the cp state. Then, choose another cell and convert it into the op state, which corresponds to a movement trial.

(3) Accept the movement trial based on the probability defined by eqn (10):

(4) Perform Steps 2 and 3 1 000 000 000 times.

(5) Every 1000 MC steps after 5 000 000 MC steps during Step 4, calculate E⁺(r^m′,P), which is the energy received by a cell when converting it into the op phase. This calculation is conducted sequentially for all cells, and the resultant values are averaged.

(6) Repeat Steps 1–4 with increasing m to obtain free energy profiles according to eqn (12).

Data availability

The data supporting this article have been included as part of the ESI.†

Author contributions

Homare Arima: conceptualization, formal analysis, investigation, writing – original draft, visualization Shotaro Hiraide: conceptualization, methodology, writing – review & editing, visualization, supervision Satoshi Watanabe: writing – review & editing, supervision.

Conflicts of interest

The authors declare no conflict of interest.

Acknowledgements

We are grateful to H. Nagano, K. Nishimoto, and K. Shota for their assistance with this study. This study was financially supported by a Grant-in-Aid for Scientific Research (B) (No. 23K23116 and 24K01235), JST, the establishment of university fellowships towards the creation of science technology innovation (No. JPMJFS2123), the Grant-in-Aid for JSPS Fellows (No. 23KJ1198), and the Hosokawa Powder Technology Foundation (Grant Number HPTF22501).

References

1. N. Höhne, M. d. Elzen, J. Rogelj, B. Metz, T. Fransen, T. Kuramochi, A. Olhoff, J. Alcamo, H. Winkler, S. Fu, M. Schaeffer, R. Schaeffer, G. P. Peters, S. Maxwell and N. Dubash, Nature, 2020, 579, 25–28.
2. D. S. Sholl and R. P. Lively, Nature, 2016, 532, 435–437.
3. J. Y. S. Lin, Science, 2016, 353, 6295.
4. S. Horike, S. Shimomura and S. Kitagawa, Nat. Chem., 2009, 1, 695–704.
5. A. Schneemann, V. Bon, I. Schwedler, I. Senkovska, S. Kaskel and R. A. Fischer, Chem. Soc. Rev., 2014, 43, 6062.
6. S. Krause, N. Hosono and S. Kitagawa, Angew. Chem., Int. Ed., 2020, 59, 15325–15341.
7. Y. Li, Y. Wang, W. Fan and D. Sun, Dalton Trans., 2022, 51, 4608.
8. I. Senkovska, V. Bon, L. Abylgazina, M. Mendt, J. Berger, G. Kieslich, P. Petkov, J. L. Fiorio, J.-O. Joswig, T. Heine, L. Schaper, C. Bachetzky, R. Schmid, R. A. Fischer, A. Pöppl, E. Brunner and S. Kaskel, Angew. Chem., Int. Ed., 2023, 62, e202218076.
9. T. K. Maji, R. Matsuda and S. Kitagawa, Nat. Mater., 2007, 6, 142.
10. J. A. Mason, J. Oktawiec, M. K. Taylor, M. R. Hudson, J. Rodriguez, J. E. Bachman, M. I. Gonzalez, A. Cervellino, A. Guagliardi, C. M. Brown, P. L. Llewellyn, N. Masciocchi and J. R. Long, Nature, 2015, 527, 357.
11. S. Hiraide, H. Tanaka, N. Ishikawa and M. T. Miyahara, ACS Appl. Mater. Interfaces, 2017, 46, 41066–41077.
12. S. Hiraide, Y. Sakanaka, S. Kajiro, S. Kawaguchi, M. T. Miyahara and H. Tanaka, Nat. Commun., 2020, 11, 3867.
13. J. Fujiki, H. Kajiro, Y. Takakura, T. Yajima and Y. Kawajiri, Chem. Eng. J., 2023, 460, 141781.
14. A. V. Neimark, P. I. Ravikovitch and A. Vishnyakov, Phys. Rev. E, 2000, 62, R1493–R1496.
15. A. Vishnyakov and A. V. Neimark, J. Phys. Chem. B, 2001, 105, 7009–7020.
16. R. Maity, M. Gholami, S. A. Peter, M. Schoukens, S. Tiriana, G. V. Baron and J. F. M. Denayer, Small, 2023, 19, 2302893.
17. T. Hiratsuka, H. Tanaka and M. T. Miyahara, ACS Nano, 2017, 11, 269–276.
18. S. Rahman, A. Arami-Niya, X. Yang, G. Xiao, G. K. Li and E. F. May, Commun. Chem., 2020, 3, 186.
19. Y. Sakata, S. Furukawa, M. Kondo, N. Horike, Y. Takashima, H. Uehara, N. Louvain, M. Meilikhov, T. Tsuruoka, S. Isoda, W. Kosaka, O. Sakata and S. Kitagawa, Science, 2013, 339, 193–196.
20. V. Bon, N. Busov, I. Senkovska, N. Bönisch, L. Abylganzia, A. Khadiev, D. Novikov and S. Kaskel, Chem. Commun., 2022, 58, 1492–1495.
21. R. Bose, V. Bon, N. Bönisch, P. Selvam, N. S. Kaisare and S. Kaskel, Chem. Mater., 2023, 35, 7825–7838.
22. S. Ehrling, I. Senkovska, V. Bon, J. D. Evans, P. Petikov, Y. Krupskaya, V. Katev, T. Wulf, A. Krylov, A. Vtyurin, S. Krylova, S. Adichtchev, E. Slyusareva, M. S. Weiss, B. Büchner, T. Heine and S. Kaskel, Front. Mater., 2019, 7, 21459.
23. S. Krause, F. S. Reuter, S. Ehrling, V. Bon, I. Senkovska, S. Kaskel and E. Brunner, Chem. Mater., 2020, 32, 4641–4650.
24. L. Abylganzia, I. Senkovska, R. Engemann, S. Ehrling, T. E. Gorelik, N. Kavoosi, U. Kaiser and S. Kaskel, Front. Chem., 2021, 9, 674566.
25. L. Abylgazina, I. Senkovska, S. Ehrling, V. Bon, P. S. Pektov, J. D. Evans, S. Krylova, A. Krylov and S. Kaskel, CrystEngComm, 2021, 23, 538–549.
26. M. Maliuta, I. Senkovska, R. Thümtmter, S. Ehrling, S. Becker, V. Romaka, J. D. Evans and S. Kaskel, Dalton Trans., 2023, 52, 2816.
27. N. Klein, C. Herzog, M. Sabo, I. Senkovska, J. Getzscgmann, S. Paasch, M. R. Lohe, E. Brunner and S. Kaskel, Phys. Chem. Chem. Phys., 2010, 12, 11778–11784.
28. F.-X. Coudert, M. Jeffroy, A. H. Fuchs, A. Boutin and C. Mellot-Draznieks, J. Am. Chem. Soc., 2008, 130, 14294–14302.
29. S. Watanabe, H. Sugiyama, H. Adachi, H. Tanaka and M. T. Minoru, J. Chem. Phys., 2009, 130, 164707.
30. H. Sugiyama, S. Watanabe, H. Tanaka and M. T. Miyahara, Langmuir, 2012, 28, 5093–5100.
31. R. Numaguchi, H. Tanaka, S. Watanabe and M. T. Miyahara, J. Chem. Phys., 2013, 138, 054708.
32. S. Hiraide, H. Arima, H. Tanaka and M. T. Miyahara, ACS Appl. Mater. Interfaces, 2021, 13, 30213–30223.
33. H. Arima, S. Hiraide, M. T. Miyahara and S. Watanabe, ACS Appl. Mater. Interfaces, 2023, 15, 36975–36987.
34. C. Triguero, F.-X. Coudert, A. Boutin, A. H. Fuchs and A. V. Neimark, J. Chem. Phys., 2012, 137, 184702.
35. F.-X. Coudert, A. Boutin, A. H. Fuchs and A. V. Neimark, J. Phys. Chem. Lett., 2013, 4, 3198–3205.
36. S. M. J. Rogge, Faraday Discuss., 2021, 225, 271–285.
37. S. Vandenhaute, S. M. J. Rogge and V. V. Speybroeck, Front. Chem., 2021, 9, 718920.
38. K. Mitsumoto and K. Takae, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, 30.
39. S. Hiraide, Y. Sakanaka, Y. Iida, H. Arima, M. T. Miyahara and S. Watanabe, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, 31.
40. L. Onsager, Phys. Rev., 1944, 65, 117.
41. B. Fultz, Phase Transition in Materials, Cambridge university press, 2nd edn, 2020, pp. 10–12.
42. D. Li and K. Kaneko, Chem. Phys. Lett., 2001, 335, 50–56.
43. A. Kondo, H. Noguchi, S. Ohnishi, H. Kajiro, A. Tohdoh, Y. Hattori, W.-C. Xu, H. Tanaka, H. Kanoh and K. Kaneko, Nano Lett., 2006, 6, 2581–2584.
44. B. Widom, J. Phys. Chem., 1982, 86, 869–872.
45. S. Watanabe, S. Hiraide, H. Kunimitsu, A. Fujiwara and M. T. Miyahara, Front. Mater., 2022, 9, 825592.
46. X.-C. Huang, Y.-Y. Lin, J.-P. Zhang and X.-M. Chen, Angew. Chem., Int. Ed., 2006, 45, 1557–1559.
47. S. Watanabe, S. Ohsaki, T. Hanafusa, K. Takada, K. Mae and M. T. Miyahara, Chem. Eng. J., 2017, 313, 724–733.
48. S. Ehrling, H. Miura, I. Senkovska and S. Kaskel, Trends Chem., 2021, 3, 291–304.
49. D. D. Wagman, W. H. Evans, V. B. Parker, R. H. Schumm, I. Halow, S. M. Bailey, K. L. Churney and R. L. Nuttall, J. Phys. Chem. Ref. Data, 1982, 11, 2.
50. R. A. d. Anjos, J. R. Vianna, J. R. d. Sousa and J. A. Plascak, Phys. Rev. E, 2007, 76, 022103.
51. D. P. Landau, Phys. Rev. B: Condens. Matter Mater. Phys., 1980, 21, 1285–1297.
52. S. Ohsaki, S. Watanabe, H. Tanaka and M. T. Miyahara, J. Phys. Chem. C, 2017, 121, 20366–20374.
53. R. Kotani, A. Kondo and K. Maeda, Chem. Commun., 2012, 48, 11316–11318.
54. T. Suzuki, R. Kotani, A. Kondo and K. Maeda, J. Phys. Chem. C, 2016, 120, 21571–21579.
55. S. Hiraide, K. Nishimoto and S. Watanabe, J. Mater. Chem. A, 2024, 12, 18193.

---

Footnote

† Electronic supplementary information (ESI) available: Derivation of eqn (12). Particle size distribution for the ELM-11 samples. Explanation of particle-size dependence based on the primitive theory and its shortcomings. Downward-convex profile inducing a gradual increase in the amount adsorbed. Particle-thickness dependence. Simulated adsorption isotherms with various M_#. Determination of the equilibrium transition pressure. Cooperative transition behavior under various conditions. Details of the nanoscale simulation. See DOI: https://doi.org/10.1039/d4ta04222k

---