Thermodynamic insight into stimuli-responsive behaviour of soft porous crystals

Knowledge of the thermodynamic potential in terms of the independent variables allows to characterize the macroscopic state of the system. However, in practice, it is difficult to access this potential experimentally due to irreversible transitions that occur between equilibrium states. A showcase example of sudden transitions between (meta)stable equilibrium states is observed for soft porous crystals possessing a network with long-range structural order, which can transform between various states upon external stimuli such as pressure, temperature and guest adsorption. Such phase transformations are typically characterized by large volume changes and may be followed experimentally by monitoring the volume change in terms of certain external triggers. Herein, we present a generalized thermodynamic approach to construct the underlying Helmholtz free energy as a function of the state variables that governs the observed behaviour based on microscopic simulations. This concept allows a unique identification of the conditions under which a material becomes flexible.

1. The equation in page 2 (last part) should have an explanation for the term of Sigma(uidNi). 2. Figure 2 has many MOF names but many of them are not dealt with in this paper. NOTT-300, Zn(CN)2, STA-12(Co) etc. The information confuses the readers and I recommend them to be moved to the supporting information or delete. Also there are many abbreviation of MOFs and it is better to describe the precise formula of each MOF in the main manuscript, otherwise it is not friendly to the non-MOF scientists. Figure 1 is not sufficient to let people understand the metalligand connectivity and assembled open structure, DMOF-1 has two ligands but the current Figure  1 does not have such information. 3. In principle they employed simple MOFs which means the structures have high symmetry, simple composition, not large cell. On the other hand, there are many important MOFs having low symmetry, large unit cell, and often various organic ligands are involved in one open structures (e.g. Angew. Chem. Int. Ed. 2010, 49, 4820-4824.). How about the (potential) solution to understand the flexible nature by the proposed method? In recent, mixed ligand approach is specifically important to tune the flexibility for application. Please provide some additional data or at least discussion. 4. References are required for the parts "...contracted phase occurs at high transition pressure of about 120 MPa and returns back to the open phase upon releasing the pressure below 20 MPa." (page 6) or "...which was found to be negative for MIL-53(Al)-F in contrast to its OH analogue." (page 8).
If the authors response these points and provide appropriate revision or additional results, I consider the manuscript is acceptable for this journal.
Reviewer #3 (Remarks to the Author): The manuscript reports on a systematic analysis of the thermodynamic behavior of flexible porous materials, in particular metal-organic frameworks, obtained through theoretical calculations based on molecular dynamics simulations.
The manuscript is well written and well organized. The work that has been carried out leads to significant results of general and wide interest to the scientific community. I believe the paper is worthy to be published on Nature Communications but some revision is needed.
I have a few concerns about the theoretical calculations. 1) Authors used a force-field (FF) that has been developed specifically for treating flexible MOFs. As demonstrated by published articles, the FF has been carefully tested, but I wonder how results are affected by the functional form of the FF. Could they try to repeat some of the calculations with other "standard" FFs and see how they work with respect to the QuickFF. 2) Prediction of thermodynamic properties through molecular dynamics requires rather long simulations to be sure to have a reliable sampling of the PES. According to the supplementary information, authors run simulations for a given time. How can they be sure that the simulation time was long enough? Could they repeat the simulations for different time scales at least for one system? 3) Authors adopts a molecular dynamics (MD) approach to predict the thermodynamic behavior of the flexible systems, but they do not mention that the same information could be obtained through lattice dynamics (LD) in a quasi-harmonic approximation, for instance. To confirm their findings and their wider reliability, authors should also try to use LD. In that respect, LD is more suitable to ab initio quantum mechanical methods and not only FF. So, this could be an even stronger evidence of the results obtained with classical simulations. 4) In the manuscript, the first time the authors state that the work is based on molecular dynamics simulations is at page 5. Before that, they generically use "microscopic approach". That's true, but in my opinion, the level of theory adopted in the "microscopic approach" should be clearly specified at the beginning of the paper.
We thank the reviewer for his/her constructive comments. We were very pleased with the assessment of the reviewer about the topic of the work which was judged to be highly suitable for Nature Communications. The reviewer encouraged us to show even more the predictive power of the proposed thermodynamic model, by studying some typical and important MOFs for which no experimental data are available. We were highly triggered by this comment and decided to conduct an in-depth study on some of the most challenging highly flexible materials (Co(bdp) and DUT-49(Cu) -as outlined in detail further in the response letter) for which no complete experimental insight on their stimuli responsive behavior is available. We decided to focus on the two proposed MOFs to make the exploitation of our thermodynamic model even more decisive. Indeed, we used the thermodynamic model to predict the pressure-induced behavior of Co(bdp) and to rationalize one of the most spectacular and unexpected guest-assisted breathing behaviors of MOFs, i.e. the negative gas adsorption observed in DUT-49(Cu). Apart from this major point of the reviewer, we also improved the manuscript by taking into account all other raised comments. Overall, we are convinced that with the input of the reviewer, we were able to bring the manuscript to a substantially higher level. We appreciate his/her input and the constructive final assessment very much. We hope that the study is now suitable to be published in Nature Communications. Below the comments of the reviewer are copied (in italics) with an answer from our side.
I was asked to review the manuscript by Vanduyfhuys and co-workers who developed a thermodynamic model that allows the computational study of the underlying thermodynamics in flexible MOFs. The target of the work, fundamental understanding of flexible MOFs, is highly topical and presents a fascinating research field of great current multidisciplinary (!) interest. Therefore, I regard the general topic of this field as highly suitable for Nat. Commun.

General considerations:
The target of the work itself, development of a model describing the underlying thermodynamics of the structural flexibility in MOFs is interesting. The study is well-organised, nicely written and the reader carefully guided through the development of the thermodynamic model of which the necessary calculations have been performed by force field based molecular simulations. Importantly, the model not only allows for quantifying the thermodynamics of an observed phase transition in a material, but further -and this the strength of the developed model-the identification of flexibility in new frameworks computationally.
My current concern directly follows up on the strength of the model, the prediction of flexibility upon a certain stimulus, which, however, has not been verified convincingly by experiment (yet). I agree with the authors that the developed model has a large potential, however, this potential need to be quantified for one or two typical and important flexible MOFs. In the current state, I feel that the study is better suited in a more specialised journal and therefore I am recommending rejection of the manuscript, however, with highly motivating the authors of resubmission. For (pressure). In this context, it will be also exciting to see if the negative gas adsorption behaviour of DUT-49(Cu) falls within one of the given categories.
Response on the predictive power of the thermodynamic model and its application to some more typical and important flexible MOFs: We studied both the stimuli responsive behavior of the Co(bdp) in absence and presence of methane molecules. For the methane-induced responsive behavior, experimental data are available, however, the pressure-induced flexibility of the empty frameworks has never been explored so far. Our thermodynamic model predicted that the Co(bdp) behaves as a triggered disperser when exposed to mechanical pressure. Furthermore, we investigated in how far our thermodynamic model enables to rationalize the observed negative gas adsorption with methane in DUT-49(Cu).
The Co(bdp) (BDP 2-=1.4-benzenedipyrazolate) was originally synthesized by the group of Jeffrey R. Long in 2008. 1 The structure consists of pairs of Co 2+ ions along a one-dimensional chain, which are tetrahedrally coordinated by N atoms from four independent BDP 2ligands. Two structures were resolved in the original paper: the Co(BDP).2DEF.H2O (as-synthesized) and the fully desolvated Co(BDP). Some cell parameters and characteristic volumes, as determined experimentally are taken up in Table  R1.1 of this review letter. The as-synthesized material was determined to have a tetragonal crystal structure with a cell volume of 2458 Å 3 . After a full desolvation, the X-ray powder diffraction pattern indicated a complete conversion of the original tetragonal structure to a contracted monoclinic phase with a substantial volume change towards 1183 Å 3 . In 2010 a spectacular structural response of the material was experimentally evidenced by Salles and co-workers upon adsorption of nonpolar molecules such as N2 and H2. 2 N2 gas sorption measurements at 77 K revealed a five-step adsorption process, as schematically shown in Figure R1.1. Later in 2015 the response of the material was studied upon methane adsorption by Mason et al. 3 It was found that at low methane pressures a contracted phase 1 with a volume of 1183 Å 3 was observed, whereas at higher methane pressures a large pore phase was obtained with a volume of 2294 Å 3 . The unit cell properties of all resolved structures are summarized in Table R1.1 of this response letter together with structural data of the simulated structures performed for this paper.  We applied the thermodynamic model on both the empty framework and structures loaded with a varying amount of methane molecules. To this end, we first constructed a new force field, with our QuickFF routine following the computational details described in the main manuscript. A series of periodic Density Functional Theory calculations were performed using VASP, on initial structures resembling the contracted pore and large pore phases as reported in reference. 2 We constructed an equation-of-state around the large pore phase and used these input data for the construction of the force field. The equation-of-state resulting from the DFT calculations at 0 K are shown in Figure R1  Structural transitions induced by a mechanical pressure were studied by constructing the thermodynamic potential in terms of the mechanical equations of state and the volume versus pressure curves. The results are shown in Figure R1.3. Interestingly two empty phases are found: a large pore phase with a volume of 2470 Å 3 and a contracted pore phase with a volume of 1060 Å 3 . The contracted pore phase is substantially more stable by about 45 kJ/mol at 300 K. This structure resembles the experimentally observed desolvated phase, which is stabilized by  interactions between the aryl rings of the neighboring BDP 2ligands. Following the nomenclature introduced in the paper, Co(bdp) belongs to the class of Type II materials and behaves as a triggered disperser (Type IIb). Experimentally this pressure induced behavior of the material was not observed yet, as the current synthesis procedures only obtained the contracted pore phase when fully desolvating the as-synthesized structure. These predictions are expected to boost the experimentalists, to develop another synthesis/activation procedure in order to prepare the material in its large pore phase and to further control its performance in a compression/decompression cycle.
We have included the results of the Co(bdp) in Figure 4(a) and discussed the results on page 6 of the manuscript. The following paragraph was added: "Similar pressure-induced flexibility is predicted from our model for the Co(bdp) material, not belonging to the MIL-53 series. Current experimental procedures only allowed to synthesize a fully desolvated structure in the contracted pore phase, whereas our model predicts the existence of a large and contracted empty pore phase which is only slightly less stable (45 kJ/mol at 300 K). This material might be interesting for further testing in a compression/decompression cycle" In a next step, we used the proposed thermodynamic model to study the methane induced flexibility of the Co(bdp) material which was experimentally observed. We adopted the procedures explained in the paper for the case where the trigger X≠P. As a result, a series of profiles were constructed with the newly developed force field. The thermodynamic potential in terms of the volume, the mechanical equations of state and () VN are shown in Figure R1.3 for a varying number of molecules per unit cell. Our model predicts that the material undergoes a phase transition from a contracted pore to a large pore phase, when initially starting from a contracted pore phase, for a loading of methane per unit cell varying between 2 and 4 molecules. The characteristics of the unit cell and volume of the various phases are indicated in Table R1.1. For phases derived from the MD simulations at 300 K, the cell parameters were obtained by averaging over the MD (NV = 0T) simulations with a volume taken from the grid point which is the closest to the minimum (i.e. with the pressure ≈1 bar). Experimentally, a transition from a contracted pore phase to a large pore phase was found for a methane pressure between 16 bar corresponding to 1 CH4 adsorbed molecule per unit cell volume of 1183 A 3 and 26 bar corresponding to 7 CH4 adsorbed molecules per unit cell volume of 2305 A 3 . 3 Our thermodynamic model agrees with the experimentally-observed methane induced flexibility of the Co(bdp).
We incorporated our results of methane adsorption in Co(bdp) in the manuscript. The results are shown in Figure 4 (c), which is repeated here in Figure R1.5. The original results on xenon adsorbed in MIL-53(Al) are still discussed in the main paper but the corresponding figures with the mechanical equations of state, pressure versus external trigger and volume versus number of particles are now shown in the Supplementary Information. As such we meet the comment of the reviewer to incorporate results on more exciting materials of the recent literature.
The text has also slightly been adapted to incorporate the results of Co(bdp) in the main paper as shown below: "Guest-induced phase transformations are illustrated for xenon adsorption in MIL-53(Al) ( Figure S2 of the SI), methane adsorption in Co(bdp) (Figure 4c) and benzene adsorption in DMOF-1(Zn) (Figure 4d). We investigated the response of all these materials upon guest loading by constructing a series of ( , ) Legendre transformation of the corresponding Helmholtz free energy to the grand canonical potential. 4 However, since we here investigate the number of particles as a trigger, we proceed by localizing the stable points at the intersection of the pressure profiles with a fixed mechanical pressure of 1 bar. Hence, the () VN response curves upon guest loading and at fixed mechanical pressure, corresponding to atmospheric pressure, are constructed (Figure 4c for methane@Co(bdp) and Figure S2  Finally we also studied the pressure induced behavior of the DUT-49(Cu) as requested by the reviewer. DUT-49(Cu) is composed of copper paddlewheels connected through 9.9'-([1,1'-biphenyl]-4,4'diyl)bis(9H-carbazole-3,6-dicarboxylate)(BBCDC) ligands and it was originally synthesized by the group of Kaskel in 2012. 7 In 2016 the same group showed a spectacular negative gas adsorption for methane and n-butane. 8 This negative gas adsorption was associated with a phase transition from the open pore phase to a contracted pore phase, characterized on the microscopic scale by the buckling of the BBCDC ligand. Very recently, Evans and co-workers performed a theoretical study on the DUT-49(Cu) material to rationalize the negative gas adsorption. They investigated the influence of the external pressure and methane adsorption on the flexibility of this material. 9 To this end, free energy versus volume curves were constructed using the procedure outlined in reference. 10 These simulations were performed in the NVh0T ensemble where the shape fluctuations are not taken into account.
Starting from the force field reported by Evans et al. we studied the mechanical response of the empty DUT-49(Cu) material using our thermodynamic model. Within the time frame of this review, it was impossible to use our own QuickFF routine to derive an in-house developed force field, given the large amount of atoms in the unit cell for which DFT calculations would have to be performed. The empty DUT-49(Cu) framework contains 1728 atoms. Because of symmetry, this number can be reduced by at least factor 4. Therefore, we started from the force field proposed by Evans to study the mechanical induced flexibility which was not yet experimentally measured. The results are shown in Figure R1.4 We find similar results as reported by Evans et al. and predict that DUT-49(Cu) can be classified as a Type IIa material and is hence a potential candidate as shock absorber. If the material is originally in its large pore phase, it will undergo a phase transition to the contracted pore phase upon external pressures higher than 50 MPa. When releasing the pressure the material would remain in the metastable contracted phase. We added the pressure-induced-flexibility behavior of DUT-49(Cu) in the main manuscript, in the section where Type IIa materials are discussed on page 6.
"DUT-49(Cu) is another example of a potential shock absorber ( Figure 4a). 7 According to our simulations, DUT-49(Cu) is predicted to have the ability to store more mechanical energy as MIL-53(Al)-F (57.8 J/g vs 7.9 J/g). This might seem contradictory, as the transition pressure of DUT-49(Cu) is lower, however, as explained in the SI, the larger volume variation (per gram material) is the main cause for the higher amount of absorbed mechanical energy. " The following section has been added to the Supplementary Information about the mechanical energy stored in MIL-53(Al)-F and DUT-49(Cu):

Approximate calculation of mechanical energy stored in MIL-53(Al)-F and DUT-49(Cu):
An approximate method to calculate the mechanical energy stored in the material during one compression/decompression cycle consists of considering the transition pressure times the volume change associated with the phase transition. This is only an approximate calculation since it assumes that the transition occurs reversibly at a fixed pressure given by the transition pressure. However, during the transition itself, the system is not in mechanical equilibrium with its environment. As a result, the transition is not a quasi-static process, which also makes it irreversible and this method is hence only approximate. Using this method, we find values of 25. Response in how far the negative gas adsorption could be detected by the thermodynamic model. This is a very interesting point, which warrants a more fundamental thermodynamic discussion. We shortly mentioned the potential extension of the model to explain negative gas adsorption as observed by Kaskel et al. 8 in the manuscript on page 10 (vide infra).
The starting point of our thermodynamic model is the Helmholtz free energy ( , , ) F N V T with state variables temperature, volume and number of particles. However, one can transform one thermodynamic potential to another by means of a Legendre transformation, as noted in the paper on page 3. 4 In an experimental set-up, one typically controls the mechanical pressure on the MOF, the chemical potential of the adsorbent species in the environment and the temperature. In many cases, those discussed below, the mechanical pressure P is equal to the vapor pressure exerted by the gas. The chemical potential and the vapor pressure are coupled by means of the equation of state. To detect the phenomenon of negative gas adsorption, one needs to construct the adsorption isotherm ( ) which yields the number of adsorbed particles in terms of the vapor pressure. A negative gas adsorption corresponds to a region where the number of particles decreases when increasing the vapor pressure. This phenomenon was experimentally observed for DUT-49(Cu) as shown in Figure R1 Such a procedure was already applied by Vanduyfhuys et al. 4 In the latter paper, the minimization was performed in two steps: first a minimization is performed in terms of the number of particles for fixed values of  and P : is an intermediate potential in which N and V are treated as independent variables and  and P are fixed parameters. As a result the thermodynamic potential is obtained ̃( ; , ), in which the chemical potential indeed acts as a trigger, but also the amount of adsorbed guest particles ( ; , ). As was mentioned earlier, the mechanical pressure P equals the vapor pressure . Therefore, from now on, we will drop either  or P from the argument list of ̃ and and assume  to be a function of through the equation of state. Furthermore, the pressure term PV in ̃ can be neglected for vapor pressures in the order of 1 bar. In the following, we ignore this term and as a result, the intermediate potential ̃( ; , ) actually becomes the grand canonical potential ( ; ) in the ( , , ) ensemble.
We illustrate the procedure for xenon@MIL-53(Al) at 300 K. Figure R1.7(a) shows the ̃( ; ), in terms of the volume whereas Figure R1.7(b) shows the number of adsorbed guest molecules in terms of the volume for various values of the vapor pressure. In principle one needs to perform the second minimization step in terms of the volume on this thermodynamic potential. By applying a same procedure as explained in Figure    We now rationalize how the profiles of ̃( ; ) and ( ; ) look in case of positive and negative gas adsorption. Suppose we start for the xenon@MIL-53(Al) initially in the large pore phase. This corresponds to the empty phase most encountered experimentally, since the material is activated at high temperature. When increasing the vapor pressure, a new minimum appears corresponding to an intermediate phase. The transition from the large pore phase to the intermediate phase only occurs when the large pore minimum in the osmotic potential disappears. Further increasing the vapor pressure finally yields back the large pore phase. In this case positive gas adsorption is observed, as can be detected by inspecting the ( ; ) profiles. As the vapor pressure increases, the structure transforms at a vapor pressure of 0.1 bar (see arrow for the profile at 0.16 bar in Figure R1.7). At this transition, the number of particles increases in case of the xenon@MIL-53(Al). To truly compare with experiment one can perform the second minimization step of the Legendre transformation and construct the adsorption isotherm. The results of this second minimization step for xenon@MIL-53(Al) are shown in Figure R1.8(b).
To test the hypothesis further we started from the free energy profiles at 120 K published by Evans et al. 9 (see Figure R1.9). Following the procedure with the Legendre transform sketched above, Figure R1.10 shows the results of applying the first Legendre transform to these profiles, i.e. the thermodynamic potential ̃( ; ) as well as the number of particles ( ; ) for various values of the vapor pressure of methane. From these results, one would also conclude that no structural transition is observed when collective behavior is assumed and the system always remains in the large pore phase, because this large pore minimum never disappears. However, it could be that the transition is not observed due to missing data  To further illustrate this, Figure R1.11(a) also shows the response of the cell volume of DUT-49(Cu) to the vapor pressure of methane and Figure R1.11(b) shows the adsorption isotherm of methane. As we can see in the adsorption isotherm, a strong increase in the amount of adsorbed methane for the large pore (blue squares) occurs at around 6 bar, and as a result the amount of adsorbed methane in the large pore (blue squares) becomes larger than in the narrow pore (red circles). Hence, if the transition occurs at pressures higher than 6 bar, negative gas adsorption would indeed occur. "In principle our thermodynamic model yields all necessary components to also explain unexpected effects such as negative gas adsorption, which was observed for DUT-49(Cu) upon exposure of methane, as explained in section 5 of the SI."

Minor points:
Is it possible to investigate the role of entropy more explicitly? Particularly, the temperature driven flexibility is determined by the balance of dispersion interactions and lattice entropic gain, hence a study of entropy would be very interesting at least for X = T.
The temperature driven flexibility is indeed believed to be governed by a subtle balance between longrange dispersion interactions in the framework and entropy contributions at finite temperature. 11 A proper treatment of entropy for flexible materials belongs still to one of the most challenging problems of materials modelling. It requires accurate sampling of the potential energy surface and methods which go beyond the harmonic oscillator approximation. [12][13][14] From the free energy estimation methods used in this paper, we can deduce the entropy contribution at a given temperature T in terms of the volume. Although it is indeed very challenging to compute the accurate entropy differences in terms of different temperatures ∆ = ( 2 ) − ( 1 ), we can, however, extract entropy differences in terms of different volumes ∆ = ( 2 ) − ( 1 ) by computing the free energy difference as well as the difference in internal energy and applying the formula ∆ = ∆ −∆ . Furthemore, we can compute ∆ for various temperatures to investigate the influence of temperature on entropy differences. Note that this procedure only allows us to compute ∆ ∆ and not ∆ .
We studied more in depth how the entropy varies in terms of the volume and for various temperatures for MIL-53(Al)-F and for xenon@Mil-53(Al). Such an analysis should allow to give more insight into the larger contribution of entropy to the large pore phase, which is thought to be governing the phase transformation. Figure R1.12 shows the decomposition of the free energy of the MIL-53(Al)-F into the internal energy and the entropy contribution for the force-field based MD simulations at various temperatures. Note that the reference level for both the internal energy and entropy is arbitrarily chosen to enable a visualization of all curves on the same plot. Hence only energy/entropy differences between different volumes ∆ are meaningful. As mentioned before we are unable to deduce the entropy explicit in terms of temperature. 14 The internal energy is almost independent of the temperature (around 2 kJ/mol variation on the large pore-contracted pore stability). Our simulations using classical molecular dynamics indicate that the entropy variation in terms of the volume is nearly independent on the temperature. Such effects were already previously observed in other systems. 15 However, further investigation is warranted on this topic, as for instance nuclear quantum effects may play an important role to determine the entropy more accurately. Current results show that the entropy for the large pore phase is substantially higher than for the contracted pore phase ( ( , ) − ( , ) ≈ 33 J/(mol.K) for a large pore volume and a contracted pore volume taken as the 300 K equilibria). At higher temperatures, the large pore phase is indeed systematically more stable compared to the contracted pore phase, since the entropy contribution () lp cp T S S  to the free energy becomes more dominant due to its linear dependency on the temperature and yields a systematically larger contribution to the overall free energy. Further investigation would be required to study in detail the temperature independence in terms of the volume but this is beyond the scope of the current paper.
A similar analysis was performed for the xenon@MIl-53(Al). A decomposition of the free energy into the internal energy and entropy is performed for various xenon loadings. In this case the internal energy does depend on the xenon concentration, which is related to the contributions from the xenon-host and xenon-xenon interactions. Also the entropy is largely dependent on the number of particles, as expected since higher loadings mean less free space and hence lower entropy.
We added in the methods section the following sentence and reference to the Supplementary Information: From the output of such molecular dynamics simulations, the time average of the instantaneous hydrostatic pressure Pi is computed, which is defined at each time step as Pi=Tr(σi)/3 with σi the instantaneous, internal stress tensor resulting in () PV and () FV profiles. More information about the procedure can be found in the work of Rogge et al. 10 For a selected set of materials, the separate contribution of internal energy and entropy was deduced from the simulations to obtain more insight into the effects contributing to the stabilization of one or the other phase (Section 6 of the SI). Figure 1 is somewhat confusing. I suggest updating the figures, in particularly the metal nodes and include the coordination environment of the metals.

Figure1: for a non-expert in the field, the inorganic building unit shown in
We updated the figure as requested by the reviewer 1. Furthermore, the figure was extended to include also Co(bdp) and DUT-49(Cu) as these material is now discussed in the revised version. The new figure is shown below. The sentence has been adapted as requested by the reviewer.
"Experimentally, framework flexibility can be followed by monitoring the response of the material to the external stimulus (schematically shown in Figure 2) with techniques such as mercury-intrusion porosimetry 16 , high-pressure X-ray diffraction 17 , differential scanning calorimetry 18 or spectroscopic techniques. 19 " Typo at page 6, very top "remain remain".

All minor corrections have been adapted.
We thank the reviewer for his/her constructive comments, which helped to improve the manuscript to a great extent. We have taken all comments of the reviewer into account and adapted the manuscript accordingly. Below the comments of the reviewer are copied (in italics) with an answer from our side. Figure 3 and the section of conclusion.

I read the manuscript carefully, and found the proposed results and discussion is of significant to understand the unique nature of soft MOF and realized that the approach is powerful to illustrate the different types of flexible behaviors in various MOFs. I found several new information which are not evaluated by the experiments, and the results herein support the experiments and prediction of the flexiblity-based functionality of MOFs. The paper mostly focuses on the three-dimensional structured MOFs but there are many flexible MOFs with two-dimensional or even one-dimensional coordination polymers (e.g. Nano Lett 2006, 6, 2581-2584.). I am not sure how the present approach could be available for the other flexible MOFs/CPs in terms of creation of force field, and the authors should add explanation about this point. I also wondered how the behavior of crystal-to-amorphous (reversible)
transformation which is recently highlighted will be explained based on the proposed approach (e.g. Acc. Chem. Res. 2014, 47, 1555-1562. or APL Materials 2014. The authors only studied on the crystalline nature but the expansion of the discussion also for glass or amorphous MOFs would become more strong to generalize the work in here. At least the authors should mention/discuss the potential or perspective to access this issue in this manuscript. The comments raised by the reviewer towards the extension to two-dimensional or one-dimensional coordination polymers, crystal-to-amorphous transformation or amorphous MOFs is particularly interesting. The proposed thermodynamic model is in principle broadly applicable, only a few assumptions are made. First the volume needs to represent a representative collective variable which enables to detect the breathing flexibility investigated here (as mentioned on Page 1 of the paper). Second, collective behavior is assumed as mentioned on Page 8 of the paper [A transition to the contracted pore phase is predicted to be only possible if the open pore minimum in the free energy profile would disappear, according to the assumption of collective behavior. 40 ]

Application towards 1D and 2D coordination polymers
While our model is conceptually applicable for 1D or 2D coordination polymers, deriving reliable force fields for these type of materials is not a straightforward task. For the mentioned materials, the structure is very much dominated by non-bonding interactions between 1D chains or between the layers. For the molecular framework Ag(tcm) (tcm= tricyanomethanide) 20 with a layer-like topology, we derived a force field based on Density Functional Theory based input. However after optimization of the structure with the new force field, the layers drifted apart, which is an indication that the non-bonding interactions from the MM3 force field, which are used most commonly complementary to the covalent force field derived by QuickFF, are not fully suited. Such a study on the non-bonding interactions in lower dimensional coordination polymers is beyond the scope of the current paper. Indeed, as mentioned in literature, non-bonding interactions can only be described well by expensive high-level calculations in both 1D 21 and layered 2D systems. 22

Application towards amorphous materials
Soft porous materials lacking long-range order or soft porous materials for which the phase transition is accompanied by a loss of long-range order cannot be directly simulated with the protocol suggested here, as our procedure relies on periodic boundary conditions. Using these periodic boundary conditions, the length scale on which the loss of order can be described is limited to the size of the unit cell, such that only short-range loss of order can be described within our model. In contrast, amorphous metal−organic frameworks (aMOFs) retain the basic building blocks and connectivity of their crystalline counterparts, and are hence highly ordered on these short length scales, but lack any long-range periodic order, as mentioned in the review of Bennett et al. 23 To date, simulation of such systems hence remains very challenging with atomistic based models, as constructing sufficiently large unit cells of these materials is computationally restrictively expensive. Hence, to describe these long-range effects, one needs to rely on an approximate description of the system with fewer degrees of freedom, such as the simplified Hamiltonian to describe the adsorption-induced layer-by-layer LP-to-CP transition in MIL-53 24 , or coarse-grained models, which have already proven to be applicable for the description of biomolecules, 25 but are only slowly being introduced for framework materials. 26 Following paragraph has been added to the main manuscript at page 4 and 10: The approach is generic and may be applied directly to any soft porous crystal exhibiting long-range order. Materials for which the phase transition is accompanied by loss of long-range order, such as amorphous MOFs, 23 cannot be straightforwardly simulated as the present approach relies on periodic boundary conditions. Effects such as amorphization under elevated pressures may however be detected both experimentally and theoretically by a broadening of the peaks in the radial distribution function. 27 Our thermodynamic model is in principle generic but relies on reliable force fields. As a typical illustration, we compared the results obtained for DMOF-1(Zn) with various force fields (section 9.1 of the SI). For some materials such as two-dimensional polymers, force field development might be very challenging due to critical dependence on the non-bonding interactions between the layers or the chains.. 21,22 This work has enough generality and novelty to explain/predict the flexible nature of various MOFs and it would contribute to explore the new functions and application in the future which would be worth to be published in this journal, but I address some points to be revised to consider the final decision as following.

The equation in page 2 (last part) should have an explanation for the term of Sigma(uidNi).
If we start from the well-known fundamental relation of thermodynamics for closed systems = − , we can extend this equation to account for adsorption of guest molecules in open systems. For each species i with Ni particles and chemical potential , an extra term can be included in the fundamental relation, which accounts for the work needed to adsorb particles in the framework. As such, we arrive at Finally, by using the definition of the Helmholtz free energy = − , we arrive at: Figure 2 has many MOF names but many of them are not dealt with in this paper. NOTT-300, Zn(CN)2, STA-12(Co) etc. The information confuses the readers and I recommend them to be moved to the supporting information or delete. Also there are many abbreviation of MOFs and it is better to describe the precise formula of each MOF in the main manuscript, otherwise it is not friendly to the non-MOF scientists. Figure 1 is not sufficient to let people understand the metal-ligand connectivity and assembled open structure, DMOF-1(Zn) has two ligands but the current Figure 1 does not have such information.
As this is a paper intended for a general audience, we agree with the reviewer that MOF-specific terminology should be introduced appropriately. However, we opted not to remove all materials from Figure 2, as they give an indication that the represented phenomena were found in various materials. Instead, we have now included a table in the Supplementary Information (as shown below) yielding the precise formula for each material mentioned in the paper.   While MOFs with a larger unit cell and a lower symmetry require more CPU time to reach convergence, our protocol does not directly depend on the size of the unit cell nor the symmetry of the material. As a proof of this statement and as requested by Reviewer 1, we included DUT-49(Cu) as a more challenging material with a larger unit cell (1728 atoms in the empty unit cell). However, our protocol requires the volume to be a good collective variable, distinguishing between the different (meta)stable states such as the large pore and contracted pore phases. Hence, one has to be careful to identify all metastable states at the given thermodynamic conditions. Whereas high symmetric materials often have only one or two (meta)stable states as a function of the volume, less symmetric materials may exhibit multiple (meta)stable states which are degenerate, which makes it more difficult to identify all (meta)stable states. In that case, one has to ensure that all (meta)stable states are sampled to obtain reliable pressure and free energy profiles. The materials from the paper proposed by the reviewer, are indeed very challenging to model, both from the perspective of force-field development (van der Waals interactions will be of crucial importance to accurately describe the interdigitated nature of the 2D layers) as well as regarding to the choice of the collective variable (the volume is not a good collective variable for gateopening). We have inserted the appropriate references in the text.

If the authors response these points and provide appropriate revision or additional results, I consider the manuscript is acceptable for this journal.
We thank the reviewer for his/her constructive comments, which helped to improve the manuscript to a great extent. We have taken all comments of the reviewer into account and adapted the manuscript accordingly. Below the comments of the reviewer are copied (in italics) with an answer from our side.
The manuscript reports on a systematic analysis of the thermodynamic behavior of flexible porous materials, in particular metal-organic frameworks, obtained through theoretical calculations based on molecular dynamics simulations.
The manuscript is well written and well organized. The work that has been carried out leads to significant results of general and wide interest to the scientific community. I believe the paper is worthy to be published on Nature Communications but some revision is needed.

I have a few concerns about the theoretical calculations. 1) Authors used a force-field (FF) that has been developed specifically for treating flexible MOFs. As demonstrated by published articles, the FF has been carefully tested, but I wonder how results are affected by the functional form of the FF. Could they try to repeat some of the calculations with other "standard" FFs and see how they work with respect to the QuickFF.
The sensitivity of the results with respect to the applied force field warrants indeed some attention. We choose to investigate the flexibility behavior of DMOF-1(Zn) more in depth to respond to the question of the reviewer. For this material, we used in the submitted paper the force field of Grosch et al. 28 In this force field, the covalent terms that describe interactions within the BDC and DABCO ligands are taken from the Generalized Amber Force Field (GAFF). 29 Covalent interactions involving the metal unit are parametrized using DFT calculations with the M06-2X functional on framework fragments. 30 Electrostatic interactions are described using point charges, which are obtained by fitting to the ab initio electrostatic potential (from the DFT framework-fragment calculations) using the CHELPG method. 31 The van der Waals interactions are modeled using a Lennard-Jones potential using the corresponding GAFF parameters.
We now also generated a QuickFF force field for the DMOF-1(Zn) framework, starting from an ab initio determined Hessian using periodic plane waves DFT calculations. Similar to the other materials in this study, we employed the PBE exchange-correlation functional together with DFT-D3(BJ) dispersion corrections. More details can be found in the section 1.1 of the SI. Furthermore we also applied a general purpose force field namely UFF4MOF. 32 UFF4MOF was developed as an extension of UFF towards MOFs, with new parameters specifically valid for MOF materials. UFF4MOF has been shown very accurate in reproducing unit cell dimensions of a series of MOFs, however it needs to be tested in how far the parameter set is accurate enough to simulate physical phenomena that are more sensitive to the specific shape of the potential energy surface and temperature corrections such as breathing. A recent assessment of various force fields to predict bulk material's properties like the bulk modulus and linear thermal expansion for a series of MOFs was performed by Boyd et al. 33 Their study only included rigid materials like IRMOF-1, IRMOF-10, HKUST-1 and UiO-66. Whereas the bulk properties were recently predicted, properties which are sensitive to the framework vibrational modes are prone to larger deviations.
For the three tested force fields the results of the thermodynamic model are given in Figure R3.1. All force fields predict that the volume decreases as the benzene loading increases. The results of the QuickFF and the Grosch force field are the closest together. The UFF4MOF force field predicts substantially different results than the other two force fields. UFF4MOF predicts a second branch in the V(N) profile, which is the most stable as soon as two or more benzene molecules are present per unit cell. This branch shows volumes around 1900 Å 3 which is not in agreement with experiment.
The results of QuickFF and the Grosch force field agree within acceptable limits. The QuickFF force field predicts in general larger volumes than the other force fields as well as experiment. Note however that for zero loading it is close to the ab initio volume to which it was fitted. Grosch et al. shows a more pronounced decrease in volume upon adsorption of 2 benzene molecules per unit cell. This is accompanied by a transition from a square to a rectangular cell. For QuickFF, the cell remains (more or less) square.
A more thorough investigation on the sensitivity of the results with respect to various force fields would have to be done to generalize the conclusions taken here. However, from the current analysis, it is not advised to use general purpose force fields like UFF4MOF to describe the flexibility induced behavior of soft porous materials. In literature, also MOF-FF -also a first-principles derived force field -was used to generate free energy profiles of the DUT-49(Cu) material in the work of Evans et al. 9 The results obtained in that work gave reasonable agreement with experiments. MOF-FF was developed by the group of Schmid et al. and was designed specifically to describe Metal-organic frameworks. 34,35 The force field energy expression is very similar to the MM3 expression, the electrostatic interactions are described by Gaussian charges and the covalent parameters are fitted to ab initio cluster data using a genetic algorithm. The MOF-FF concept assumes also a building block approach and is in this sense similar to the first generation force fields developed with QuickFF. 36 However current QuickFF routines enable to generate the input data from periodic Density Functional Theory calculations as it was performed in this paper.
We added the part on the force field testing in the Supplementary Information (section 9.1) of the paper and made a critical note on the force field sensitivity in the main paper in the Methods section.
"The thermodynamic model is in principle generic but relies on reliable force fields. For DMOF-1(Zn) we compared the results obtained with various force fields (section 9.1 of the SI). For some materials such as two-dimensional polymers, force field development might be very challenging due to critical dependence on the non-bonding interactions between the layers or the chains. 21 2) Prediction of thermodynamic properties through molecular dynamics requires rather long simulations to be sure to have a reliable sampling of the PES. According to the supplementary information, authors run simulations for a given time. How can they be sure that the simulation time was long enough? Could they repeat the simulations for different time scales at least for one system?
To study the convergence of our simulations as a function of the simulation time, we have carried out longer simulations for MIL-53(Al)-F at 100 K, 300 K and 600 K, extending the simulation up to 5 ns compared to the previously reported 1.15 ns simulations while keeping the equilibration time of 100 ps constant. We have chosen this material for this in-depth study such that we can also investigate the effect of temperature on the convergence of the simulations. For these three temperatures, the root-mean-square deviation (RMSD) was calculated for three key properties, using the results from the 5 ns simulation as the reference: (i) the pressure profile before fitting, i.e. containing the raw data obtained from the MD simulation, (ii) the pressure profile after fitting, and (iii) the free energy profile obtained by thermodynamic integration of the pressure profile. The resulting RMSDs are reported in Tables R3.1 (100 K), R3.2 (300 K), and R3.3 (600 K). Across all temperatures, we observe an RMSD of 3-6 MPa on the pressure before fitting when comparing the 1 ns simulation (indicated in gray, close to the 1.15 ns simulations reported in the main manuscript) with the longer 5 ns simulation (assumed to be converged), which slightly increases with increasing temperature (3.5 MPa at 100 K versus 5.2 MPa at 600 K). This may be attributed to the size of the available phase space, which increases at higher temperatures, whereas the slower exploration of the phase space at 100 K, due to the lower velocities, does not lead to a worse sampling. By fitting the pressure, the RMSD drops to about 1-2 MPa for all temperatures, well within the accuracy one can expect to obtain with this method. This RMSD in the pressure results in an RMSD in the free energy of about 0.1 kJ/mol when comparing the 1 ns simulation with the 5 ns simulation. To get more insight in which volume regions of the pressure and free energy profiles are most prone to these small imprecisions, we report in Figure R3.1 the pressure profile before fitting (pane a) and the free energy profile obtained by thermodynamic integration (pane c) as a function of the unit cell volume for the six different simulation times at 300 K. Since the different profiles are virtually indistinguishable on this scale, panes (b) and (d) report the difference in pressure, respectively free energy, with respect to the 5 ns simulation. We observe that the deviations decrease with increasing simulation time. Finally, in Tables R3.4 (100 K), R3.5 (300 K), and R3.6 (600 K), the volumes of the metastable CP and LP states as well as the LP-to-CP and CP-to-LP transition pressures extracted from the pressure profiles with different simulation times are reported. We observe that the deviations due to limited simulation time do only affect these properties to a minor extent. In conclusion, the earlier reported 1.15 ns simulations have converged sufficiently to yield accurate results within a few MPa and a few tenths kJ/mol for the pressure and free energy profiles, respectively, yielding volumes for the metastable states which are accurate within a few Å³.
For DMOF-1(Zn) loaded with benzene molecules and for MIL-53 loaded with xenon atoms, the simulations were also extended and no noticeable differences were observed in the pressure and free energy profiles.
The results of this analysis are taken up in the Supplementary Information section 9.2.    3) Authors adopts a molecular dynamics (MD) approach to predict the thermodynamic behavior of the flexible systems, but they do not mention that the same information could be obtained through lattice dynamics (LD) in a quasi-harmonic approximation, for instance. To confirm their findings and their wider reliability, authors should also try to use LD. In that respect, LD is more suitable to ab initio quantum mechanical methods and not only FF. So, this could be an even stronger evidence of the results obtained with classical simulations.
As requested by the reviewer, we also calculated the vibrational entropy via lattice dynamics calculations. More information on the procedure may be found in the reviews of Butler et al. and Tkatchenko et al. 12,13 This approach is critically dependent on the accurate calculation of the Hessian matrix, which is called the dynamical matrix for periodic systems. It contains the second order derivatives with respect to the geometry of the system. As a result one obtains a series of phonon modes and phonon eigenvectors. To properly include all non-negligible lattice modes a supercell approach might be necessary, as it is important not to introduce artifacts between periodic images.
To obtain thermal corrections at various volume points we applied the quasi-harmonic approximation, where a series of structures are optimized for a number of volume points. The phonons are still harmonic in this case but thermal corrections which are dependent on the volume are taken into account as the phonons become volume dependent.
We applied this procedure both for MIL-53(Al) and xenon@MIL-53(Al) systems using force fields. In case of force field calculations, the Hessian and forces can be calculated in an analytical way and is thus less prone to numerical uncertainty. We used a 1x2x1 unit cell which is necessary to describe all potential energy terms correctly in the force field expression. Calculations on a 2x4x2 unit cell provide very similar results, from which we conclude that the 1x2x1 supercell is sufficiently large. We choose a grid spacing of 25 Å 3 . The results are shown below. For MIL-53(Al) we find that the quasi-harmonic approximation yields fairly good agreement with the MD based results, at least if one uses a force field in which the bonds and bends are described using harmonic terms (which is the case in the submitted version of the paper). To investigate in how far this influences the obtained profiles, we also constructed similar profiles including anharmonic contributions, by adding anharmonic corrections from the MM3 force field to the bond and bend terms. 37 It is clear that the quasi-harmonic profile is in this case in worse agreement with the profile obtained from MD simulations, as only in the latter approach it is possible to correctly sample anharmonicities of the potential energy surface.
For the xenon@MIL-53(Al) system the quasi-harmonic oscillator approximation fails to describe the free energy profile accurately. This can be ascribed to two effects. First of all, it is much more difficult to get smooth normal mode frequencies as function of volume because the xenon guest molecules are only very weakly bound to the framework. This can be observed by the discontinuities in the QHA data points. Second, due to the weak bonding of xenon molecules to the framework, the xenon molecules are able to sample the entire pore of the framework, giving rise to large entropy contributions. However, this high translational degree of freedom is not captured by QHA in which the xenon atoms are restricted to a small region around their equilibrium positions. We added following sentence in the main manuscript on Page 10 : "Alternatively free energy profiles may be obtained in the quasi-harmonic approximation. This was done for some materials under study in this paper, more information can be found in section 10 of the SI" 4) In the manuscript, the first time the authors state that the work is based on molecular dynamics simulations is at page 5. Before that, they generically use "microscopic approach". That's true, but in my opinion, the level of theory adopted in the "microscopic approach" should be clearly specified at the beginning of the paper.
We changed the sentence at the bottom of page 3 towards : "Herein we present a microscopic approach based on classical molecular simulations, to construct the Helmholtz free energy and to uniquely determine the macroscopic response of a material upon stimuli such as mechanical pressure, temperature and adsorbed guest molecules."