Efficient Construction of Free Energy Profiles of Breathing Metal-Organic Frameworks Using Advanced Molecular Dynamics Simulations.

In order to reliably predict and understand the breathing behavior of highly flexible metal-organic frameworks from thermodynamic considerations, an accurate estimation of the free energy difference between their different metastable states is a prerequisite. Herein, a variety of free energy estimation methods are thoroughly tested for their ability to construct the free energy profile as a function of the unit cell volume of MIL-53(Al). The methods comprise free energy perturbation, thermodynamic integration, umbrella sampling, metadynamics, and variationally enhanced sampling. A series of molecular dynamics simulations have been performed in the frame of each of the five methods to describe structural transformations in flexible materials with the volume as the collective variable, which offers a unique opportunity to assess their computational efficiency. Subsequently, the most efficient method, umbrella sampling, is used to construct an accurate free energy profile at different temperatures for MIL-53(Al) from first principles at the PBE+D3(BJ) level of theory. This study yields insight into the importance of the different aspects such as entropy contributions and anharmonic contributions on the resulting free energy profile. As such, this thorough study provides unparalleled insight in the thermodynamics of the large structural deformations of flexible materials.

Introduction

Metal–organic frameworks (MOFs) have attracted enormous attention during the last 20 years due to their versatility and their potential for various applications such as gas storage, separation, catalysis, and controlled drug delivery.[1−7] An important feature of some of these MOFs is their so-called breathing behavior, i.e., their ability to undergo large structural deformations upon external stimuli such as temperature, mechanical pressure, gas adsorption, and light absorption, while maintaining their structural integrity.[8−12] These structural transitions are accompanied by large volume changes up to 40%. The so-called MIL-53 materials, exhibiting a wine-rack topology, are generally considered as prototype flexible frameworks.[13] The structure is composed of inorganic metal hydroxide chains [M3+(μ2–OH)]∞ connected by terephthalate linkers, forming one-dimensional channels. A variety of metal nodes, such as Al3+, Ga3+, In3+, Sc3+, V3+, Cr3+, and Fe3+, and organic linkers have been employed to synthesize the MIL-53 framework.[13−19] The specific nature of these building blocks has a decisive role in the framework breathing.[11,13,14,20−26] In this paper, we focus on the extensively investigated MIL-53(Al) framework, which has the ability to switch reversibly between a closed pore (cp) and large pore (lp) phase in terms of temperature[20] and pressure[24] (see Figure ). Moreover, a very large hysteresis in terms of temperature was observed, with a transition from the lp to the cp phase between 125 and 150 K, while the transition from the cp to the lp phase occurs between 325 and 375 K. The experimental results point toward a slow transition from the lp to the cp phase and vice versa, which clearly indicates that the transition is activated and hence a rare event at these temperatures.

Figure 1

Schematic representation of the free energy profile of MIL-53(Al). The free energy profile exhibits two local minima centered around the closed pore (cp) and large pore (lp) structure with respect to the breathing behavior of MIL-53(Al). Moreover, the molecular structures of MIL-53(Al) in the cp and lp phases are shown.

To date, there is no comprehensive theory that fully predicts the necessary conditions to induce breathing. Nevertheless, some empirical rules were formulated that should be met to allow the structure to breathe.[10−12,27] In literature, various (semi)empirical thermodynamic models have been developed to describe structural transitions in flexible MOFs and the influence of adsorbed guest species on these transitions.[28−33] A key property in the various thermodynamic models is the free energy profile of the empty host, i.e., the MOF in absence of any guest molecules (schematically shown in Figure ). Although indirect information related to the free energy profile can be observed from experiment, such as in the seminal work of Coudert et al. based on adsorption isotherms,[28] the complete free energy profile cannot be extracted from experiment. In this respect, the use of computational techniques to construct free energy profiles of the empty host is indispensable. These free energy profiles form the key to functionalize the materials with various metal nodes and linker patterns for exhibiting the desired flexible behavior.

The prediction of the free energy profile at finite temperature is a challenging task[28] as it requires a good estimate of both the potential energy surface on the one hand and the thermal corrections and entropy contributions on the other hand. While a series of ab initio calculations have been performed to estimate the energy difference between the (meta)stable states of flexible MOFs,[34−39] only a few papers report free energy profiles for flexible materials. Cockayne used a quasi-harmonic approximation to construct the free energy profile as a function of the cell volume for MIL-53(Cr).[40] Haigis et al. studied the hydrothermal breakdown of MIL-53(Ga) by means of first-principles molecular dynamics (MD) using well-tempered metadynamics.[41] Some of the present authors constructed free energy profiles of various MOFs using MD simulations and thermodynamic integration based on force field methods,[32,42−44] an approach also used to study the thermodynamics of the negative gas adsorption in DUT-49(Cu).[45] In any case, when using an MD approach, advanced sampling methods need to be adopted to simulate the phase transformation, as it is a rare event. Herein, a thorough investigation is performed to construct free energy profiles with different molecular simulation methods. Various sampling methods to generate free energy profiles of flexible materials at finite temperature are tested. Our study includes free energy estimation methods based on metadynamics, umbrella sampling, variationally enhanced sampling, thermodynamic integration, and free energy perturbation.

A second ingredient which is decisive in the free energy profile is the description of the potential energy surface. The most accurate result can be obtained with a Hamiltonian which contains all interactions necessary to describe the detailed electronic structure. Herein, the specific way that noncovalent interactions are taken into account plays a vital role.[34] Nevertheless, due to the computational expense associated with first-principles calculations, phase transitions are often described using force fields. In this paper, we use a two-step procedure, where first the efficiency of various advanced MD methodologies for free energy computations is assessed based on a force field description of the potential energy surface. At second instance, the most efficient free energy method will be adopted to generate from first principles the free energy profile at finite temperature. The advantage of the proposed method lies in the fact that entropies and thermal corrections may be assessed beyond the harmonic oscillator approach at high accuracy. To the best of our knowledge, such a procedure has not been performed before.

Theoretical Background

From basic statistical physics, the free energy F at a temperature T is defined from the partition function Z aswith N the number of particles characterized by their positions r and momenta p, h the Planck constant, and β = 1/kBT. The Hamiltonian contains the kinetic energy and the potential energy of all particles. Depending on the boundary conditions, e.g., fixed volume or fixed pressure complemented by fixed temperature and fixed number of particles, one defines either the Helmholtz free energy or the Gibbs free enthalpy. In practice, we do not compute absolute free energies, but free energy differences. Computation of absolute free energies would require an integration over the entire phase space which is only feasible for simple systems.

A free energy calculation consists of three components: (i) an appropriate Hamiltonian to describe the potential energy of the system, (ii) a sampling protocol to generate a representative ensemble of configurations under the given boundary conditions, and (iii) methods to estimate the relevant free energy differences.[46−48] Hereafter, we explain each of the three components that contribute to the free energy calculation.

The potential energy surface may be constructed with first-principles methods such as density functional theory (DFT). Although DFT calculations are computationally attractive for large systems, it is not straightforward to use DFT as a standard method for free energy estimations, since large portions of the configuration space need to be sampled which requires numerous evaluations of the potential energy of the system. An attractive alternative description of the Hamiltonian, especially to reach large time and length scales, is based on force fields. The latter neglect the detailed electronic structure of the system and give an analytical expression for the potential energy in terms of the atomic positions. In this work, various advanced sampling methods are tested at first instance using in-house developed force fields derived from ab initio training data.[49] In some of our previous contributions we showed that our force field approach is able to capture the flexibility of MOFs.[32,42−44] Afterward the most efficient sampling method is further used for the construction of a free energy profile based on first principles such as DFT methods.

The sampling protocol to explore the phase space can rely on MD simulations, Monte Carlo simulations, or a combination hereof. In this work, MD-based methods are used. However, obtaining a sufficient sampling of the configuration space becomes troublesome for systems exhibiting metastable states which are separated by high free energy barriers as the probability to overcome these barriers is very low. A variety of sampling methods has been developed to enhance the sampling probability of highly activated regions. For some recent reviews on the topic we refer to Refs. [48 and 50]. Here, we only discuss those aspects which are necessary for the current work. One can distinguish between techniques which enhance the sampling of all degrees of freedom[51−54] and those methods which enhance the sampling along certain degrees of freedom, the so-called collective variables. The latter variables are a function of the high-dimensional microscopic coordinates and need to be chosen judiciously to sample the various interesting regions of the phase space, including those separated by high free energy barriers. The particular choice of collective variables is crucial to sample all metastable states.[55,56] The sampling along these collective variables can be enhanced by locally increasing the temperature,[57,58] by introducing a bias potential,[59,60] or by constraining the collective variable at intermediate (unstable) states.[61,62] In this work, we are interested in constructing a free energy profile for the breathing behavior of MIL-53(Al). This structural transition is characterized by a substantial volume change. As such, the volume is introduced as the collective variable.

We applied various free energy methods to estimate free energy differences. A summary of the applied methods is shown in Figure together with some key expressions. More information on the underlying statistical physics is given in Section 1 of the Supporting Information. In general, one can distinguish between equilibrium and nonequilibrium methods. In the first class, sufficiently long simulations are performed in relevant regions of the phase space. Methods in this category are free energy perturbation,[61] thermodynamic integration,[62] multistate Bennet acceptance ratio,[63] and umbrella sampling.[60] Herein, we apply the free energy perturbation method (FEP), thermodynamic integration (TI), and umbrella sampling (US). For FEP and TI, the sampling protocol relies on sampling regions of the phase space corresponding to different volume states. To this end, various fixed volume simulations are performed in the (N, V, σ = 0, T) ensemble. More information on the various ensembles is given in Ref. [42]. In FEP, the free energy differences are estimated by calculating the ensemble average of the system in state A corresponding to volume V (see panel 6 in Figure ). In TI, free energy differences are obtained by integrating over the averaged free energy derivative in terms of the volume, which is the negative of the pressure (see panel 5 in Figure ). In this case also a series of (N, V, σ = 0, T) simulations are performed for different volume states, from which the average isotropic pressure ⟨P(V)⟩ is obtained. More details on this method can be found in the work of Rogge et al.[42] In the umbrella sampling method, an external potential is added to the true Hamiltonian to enhance the sampling in low probability regions. In practice, various simulations are run in the (N, P, σ = 0, T) ensemble where a harmonic bias potential is systematically added (see panel 2 in Figure ) that forces enhanced sampling in the phase space around the reference volume V. Free energy differences are subsequently obtained with the aid of analysis methods such as the weighted histogram analysis method (WHAM)[64] and the dynamic histogram analysis method (DHAM).[65] The latter methods reconnect the local information from the biased MD runs to the unbiased system to obtain the global free energy profile of the unbiased system. For more details on these analysis methods we refer to Refs. [64−66].

Figure 2

Schematic representation of the free energy methods considered in this work. Each panel represents a different free energy method. Within each panel, the top figure shows the simulation result and the bottom figure concerns the estimated free energy profile. The color coding is black for the unknown free energy, blue for the simulation results and the estimated free energy, and red for the sampling methods. An extended scope on the theory of the different free energy methods can be found in Section 1 of the Supporting Information.

The other two applied methods belong to the so-called nonequilibrium methods, where equilibrium free energies are estimated from nonequilibrium simulations. Pioneering in this respect is the remarkable Jarzynski equality where a free energy difference between two states can be obtained from the average work done on the system by forcing the transition in a finite time.[67] Within this class of methods we distinguish, among others, local elevation,[68] coarse MD,[69] adaptive bias,[70] and adaptive force biasing.[71] Another nonequilibrium method which has received considerable attention is metadynamics.[72] This method can be considered as the finite temperature extension of the Wang–Landau approach.[73] Moreover, several notable extensions of metadynamics were developed: the well-tempered metadynamics method[74] and the variationally enhanced sampling method.[75] An extensive review of related methods is given in the recent paper of Valsson et al.[76]

Within this work both metadynamics (MTD) and variationally enhanced sampling (VES) were applied. In MTD, a history dependent potential is added to the energy surface (see panel 3 in Figure ). Gaussian hills are added to the bias potential with a preset height h and width w. Their locations V stem from the visited volume states during the MD simulation. The total number of hills depends on the update time T and the total simulation time t. An MTD simulation is converged when the sum of the Gaussian hills exactly compensates the unknown free energy profile, i.e., when all volume states become uniformly distributed under the combined action of the potential energy of the system and the bias potential. The free energy profile is systematically estimated from the negative bias potential. VES is a newly introduced enhanced sampling concept, which is particularly useful in sampling high-dimensional free energy landscapes and systems that are not easily described by a small number of collective variables. The variational principle permits us to specifically target the region we are interested in through the use of a “target distribution”. No approximate form for the bias potential is adopted; on the contrary, the bias is constructed variationally by minimizing a functional of the bias potential, which is equivalent to minimizing the Kullback–Leibler divergence between the sampled distribution and the target distribution. This is accomplished by an efficient stochastic optimization algorithm. In practice, we employ a preset number of localized Gaussian contributions with some width w and varying height h, which is updated at each iteration of the optimization algorithm (see panel 4 in Figure ). Once convergence is reached, the bias potential is in a quasi-stationary state. The advantage of this method is the flexibility of the sampling method. There is a large flexibility in choosing the target distribution which allows us to tailor the sampling. For more information on VES, we refer to Refs. [75−77].

Computational Details

Force Field Based Molecular Dynamics Simulations

All the MD simulations are performed using the in-house software code Yaff[78] in either the (N, V, σ = 0, T) ensemble or the (N, P, σ = 0, T) ensemble. This notation for the thermodynamic ensembles was introduced recently in the work of Rogge et al. and corresponds to simulations where either the volume or the pressure is controlled, whereas the unit cell shape in both cases is allowed to fluctuate.[42] This is achieved by bifurcating the isotropic pressure P and the anisotropic stress σ on the one hand and the cell shape h0 and volume V on the other hand. The temperature and pressure are controlled using the Nosé–Hoover[79] chain thermostat with three beads and a time constant of 100 fs and the Martyna–Tuckerman–Tobias–Klein barostat[80,81] with a time constant of 1 ps. Depending on the specific free energy method, either the volume or the pressure is controlled as shown schematically in Figure . All simulations were carried out at a temperature of 300 K and at a pressure of 1 MPa in the case of the (N, P, σ = 0, T) ensemble. The potential energy of the material is described using an in-house developed force field.[49] The MIL-53(Al) material is modeled using a supercell consisting of two unit cells doubled along the metal-oxide chain without including symmetry. The electrostatic interactions are computed using an Ewald summation with a real-space cutoff of 15 Å, a splitting parameter of 0.213 Å–1, and a reciprocal space cutoff of 0.320 Å–1. For the van der Waals interactions a smooth cutoff at 15 Å is applied. The MD simulations are performed using the velocity Verlet integration scheme with a time step of 0.5 fs. The total simulation time amounts to 500 ps and 1 ns for the equilibrium and nonequilibrium free energy methods, respectively. No equilibration time was taken into account since the input structures originate from dynamical simulation snapshots.

For FEP and TI, independent (N, V, σ = 0, T) simulations are performed, where the unit cell volume is chosen between 721 Å3 and 1551 Å3 with a grid spacing of 5 Å3 (ΔV in Figure ). For US, (N, P, σ = 0, T) simulations are performed, where we apply an umbrella potential with a force constant k = 300 K kB/1000 Å–6 centered around unit cell volumes in the range [721 Å3, 1551 Å3], again with a step size of 5 Å3. For the MTD simulations the bias potential has an update time T of 0.6 ps with a Gaussian hill characterized by a height h of 1 kJ/mol and a width w of 50 Å3. For the VES simulations we use Gaussian kernel functions which are centered around the unit cell volumes in the range [721 Å3, 1551 Å3] separated by 50 Å3 (ΔV in Figure ). The width w of each Gaussian kernel function is fixed at 50 Å3, and the height h is updated via a stochastic optimization scheme with an update parameter μ of 1 kJ/mol and an update time T of 8 ps. An overview of all key parameters specific for each free energy method is summarized in Table SI1 of the Supporting Information, and their optimal numerical values are provided in Table SI2. The latter have been determined after an extensive computational testing of the system and correspond to the reference free energy profiles. A further refinement of the parameters would no longer induce noticeable changes in the free energy profile within realistic convergence criteria.

The sensitivity of each method to a change of the parameter setting is a measure of the efficiency of the method. We gave preference to a procedure where we gradually make the parameter settings less stringent. We systematically evaluate the deviation with the reference free energy profile induced by the parameter change. The deviations are related to an incomplete sampling of the phase space and are measured by means of two sampling error expressions, which are introduced and discussed in Section 7 of the Supporting Information. Less stringent parameter settings are obtained by reducing the simulation time to a minimum of 50 ps. Its impact on the reproducibility of the free energy profile is measured and determines the efficiency of one specific free energy method. Furthermore, other parameters which are changed and which decrease the computational effort are the following: for the equilibrium methods, the volume grid spacing (ΔV) is enlarged from 5 Å3 to 50 Å3, while for the nonequilibrium methods, we have altered the ratio of the update time and update parameter. In the case of MTD, we have run simulations for and in the case of VES for .

Static Periodic Density Functional Theory Simulations

In order to obtain more insight into the influence of the potential energy surface description on the free energy profile, we performed a series of periodic plane-wave density functional theory (DFT) calculations. For the first-principles calculations we used the projected augmented wave[82] (PAW) method implemented in the Vienna Ab Initio Simulation Package[83] (VASP). The exchange and correlation behavior of the electrons is modeled with a generalized gradient approximation functional constructed by Perdew, Burke, and Ernzerhof (PBE).[84] Furthermore, we employ two schemes to include van der Waals interactions: the Grimme DFT-D3 interactions[85] with Becke–Johnson damping[86] (PBE+D3(BJ)) and the MBD scheme of Tkatchenko[87−89] (PBE+MBD). On the basis of the stable geometries for the cp and lp phase, a normal-mode analysis was performed to estimate the free energy difference between the two (meta)stable states. The normal-mode analysis was performed using TAMkin.[90]

First-Principles Based Molecular Dynamics Simulations

Finally, to construct a free energy profile for the structural transformation including all anharmonic contributions, we performed umbrella sampling simulations in combination with a DFT-based description of the potential energy surface. To that end, we have interfaced our in-house molecular dynamics software code Yaff with VASP, where the former serves as the MD engine and the latter as the first-principles calculation tool. At each step, VASP yields the energy and forces to the Yaff engine. Subsequently, first-principles based MD simulations in the (N, P, σ = 0, T) ensemble are performed with the same thermostat and barostat settings as in the force field simulations. The potential energy surface is constructed with the PBE+D3(BJ) scheme (as outlined in Section ). For the US simulations, harmonic potentials with a force constant k = 300 K kB/250 Å–6 centered around volumes in the range of [750 Å3, 1600 Å3] are used. The simulation time of each umbrella simulation is minimally 0.75 ps.

An important effect which needs to be taken into account when constructing a free energy profile with the aid of the PAW method is the Pulay stress.[91,92] The Pulay stress arises because the plane wave basis set is not complete with respect to the changes in the volume. A detailed discussion of the Pulay stress and its impact on the optimization of flexible frameworks was given in the work of Vanpoucke et al.[92] Herein, we compute Pulay stresses at a discrete number of volume points. Based on the electronic energy profile, E(V), and pressure, PVASP(V), as a function of the volume, we compute a Pulay stress profile: . We find a Pulay stress which follows a 1/V relation (see Figure SI1 in the Supporting Information), in line with Ref. [91]. By means of integration, we can compute a free energy FPulay(V) related to this Pulay stress. Subsequently, we correct our free energy profile, FVASP(V), obtained from the umbrella sampling simulations, with the volume dependent free energy due to the Pulay stress to obtain the correct free energy profile: F(V) = FVASP(V) – FPulay(V). Note that the latter shift is far from a negligible correction, since the Pulay stress differs approximately 0.25 GPa between the cp and lp phase of MIL-53(Al) at the PBE+D3(BJ) level.

Results and Discussion

Assessment of Free Energy Methods at the Force Field Level

Determination of the Metastable States, Transition Structure, and Free Energy Differences

Using the various free energy methods outlined in Figure , free energy profiles as a function of the volume were constructed for MIL-53(Al). These results are summarized in Figure . All applied methods predict the bistable behavior of MIL-53(Al); i.e., free energy minima are found at the cp and lp configurations. Averaging over the free energy profiles obtained with the various methods, the cp and lp states are characterized by a unit cell volume of 819 Å3 and 1449 Å3, respectively. The values are in reasonable agreement with the experimental unit cell volumes of 886.9 Å3 and 1430 Å3 at 295 K.[20] Moreover, we can determine the transition structure at the local maximum of the free energy profile, characterized by a unit cell volume of 1270 Å3. A complete overview of the volumes of the two phases, cp and lp, and the transition state predicted by the different methods is provided in the Supporting Information (see Table SI3).

Figure 3

Overview of the free energy profiles for the breathing behavior of MIL-53(Al) at 300 K obtained via the enhanced sampling methods considered in this work and using a force field description of the potential energy surface. The error bars represent one standard deviation.

Based on these definitions for the cp, lp, and transition structure, we can determine the corresponding free energy differences between any two of these phases. The force field calculations indicate that the cp structure is more stable than the lp structure at 300 K, with a free energy difference (ΔF) of 26.57 kJ/mol. For the transition from the cp toward the lp structure we find a free energy barrier of 28.25 kJ/mol (ΔF‡) and a barrier of 2.34 kJ/mol for the opposite transition (ΔF‡). The reported free energy differences stem from averaging over the different free energy methods. The individual free energy differences obtained with the various methods are summarized in Table SI4 of the Supporting Information.

At first instance it is interesting to qualitatively compare our results with the experimental results of Liu et al. and Yot et al.[20,24] Yot et al. observed a pressure-induced transition from the lp to the cp phase at room temperature which was furthermore irreversible. This suggests that the cp phase is a stable minimum in the free energy surface at room temperature. This is confirmed by our force field derived free energy profile. Liu et al. observed a transformation from the lp to the cp phase at approximately 125 to 150 K and from the cp to the lp phase around 325 to 375 K, where the latter transformation is governed by slow kinetics. Based on the force field description, we do not observe a global minimum corresponding to the lp state in this temperature range (see Figure SI2 of the Supporting Information). Thus, the free energy profiles qualitatively show the expected bistable behavior. The quantitative mismatch with experiment at the force field level may be ascribed to the particular description of the nonbonding interactions, as shown earlier by the authors.[49,93] With the force field, an electronic energy difference between the lp and cp state of 42 kJ/mol is predicted.[49] To assess the influence of the description of the potential energy surface further, we calculated the energy difference between the two phases using periodic DFT calculations. In literature some first-principles energy differences have already been reported. Those energy differences are summarized in Table SI5 of the Supporting Information. The reported energies are electronic DFT energies including dispersion corrections. Hence, no zero-point energies were taken into account, which would lead to a reduction of the lp–cp energy difference with approximately 1–2 kJ/mol. Large variations in the reported energy differences are noticeable. Walker et al. found values between 34 and 42 kJ/mol depending on the use of PBE or B3LYP as DFT functional and the particular choice of dispersion corrections.[34] Moreover, they report a value of 72 kJ/mol at the vdw-DF level. Ling and Slater obtained a much lower energy difference of 14.4 kJ/mol using a HSE06 functional with D3 corrections.[38] Herein, we find energy differences of 26.8 and 25.5 kJ/mol employing PBE+D3(BJ) and PBE+MBD calculations. These results show that the energy differences are very sensitive to the applied electronic level of theory and the employed dispersion scheme and thus should be interpreted carefully. An in-depth study of the level of theory used for the description of the potential energy surface is beyond the scope of the current paper.

To predict the free energy differences at finite temperature between the two phases based on static calculations, one could at first instance use the harmonic oscillator approximation. It must be acknowledged that this entails a serious approximation as the material possesses several anharmonic motions. Using this harmonic oscillator approximation, we find a free energy difference of 18.02 kJ/mol at 300 K based on the PBE+D3(BJ) scheme.

Summarizing, the analysis so far shows that force field based enhanced MD simulations succeed in predicting a bistable free energy profile. This approach also fully captures the anharmonic corrections to the free energy profile. Nevertheless, the exact energy difference between the two phases is sensitive to the used dispersion scheme, i.e., MM3 in our force field.[94] In general, DFT calculations including dispersion corrections seem to predict smaller free energy differences with respect to the force field results, but the former are prone to the use of the harmonic oscillator approach in the evaluation of the entropic contributions. Ultimately, at the end of this paper a fully dynamical free energy profile is constructed based on a first-principles determined potential energy surface and with proper inclusion of all anharmonic contributions. Such a simulation is, however, computationally very expensive. Therefore, we first assess the various free energy schemes and select the most efficient one which is then used in the construction of the free energy profile fully from first-principles.

Assessing the Imprecision and Efficiency of Each Free Energy Method

In this section a careful balance will be made between the efficiency of each free energy method and the precision of the obtained profile. In general, errors may be present due to inaccuracy or imprecision. Inaccuracy results from specific choices to model the system, such as the description of the potential energy surface or the unit cell size. Imprecision relates to an incomplete sampling of the phase space. The definitions as given here are inspired by the work of Lejaeghere et al.[95] Within this work we investigate the sampling error or the precision of the obtained free energy profiles. Formally, statistical mechanics methods provide a framework for a precise computation of the free energy within a given set of system-specific choices, provided we can sample the entire phase space. In reality, computer simulations are finite which introduces sampling errors or imprecision on the obtained results. In general, sampling errors can originate from two contributions: incomplete sampling in the direction of the collective variable and incomplete sampling of the other degrees of freedom. In Table , a summary is given of the various parameters introduced in the free energy methods with an indication of whether they mostly contribute to the sampling error along or orthogonal to the collective variable. The procedure allowing a quantification of the impact of an adapted parameter setting on the free energy profile has already been outlined in the Computational Details (Section ).

Table 1

Parameters which Contribute the Most to the Sampling of the Collective Variable (∥) and the Orthogonal Degrees of Freedom (⊥) for the Different Free Energy Methodsa

	∥ sampling	⊥ sampling
TI	ΔV	t
FEP	ΔV	t
US	ΔV, k, t	t, k
MTD	t	Tu , w, h
VES	t	Tu  , w, μ

The parameters in bold are the ones by which increasing their value gives less sampling, whereas the other parameters yield more sampling upon increasing their value. The parameters are defined in Figure .

To determine the sampling error originating from all degrees of freedom apart from the collective variable, n free energy profiles are constructed for each parameter set. For the equilibrium free energy methods, i.e., for TE, FEP, and US, 50 free energy profiles are constructed by bootstrapping the simulation data (see Figure ). For MTD and VES, 20 free energy profiles are computed from independent simulations. The sampling error is estimated from the standard deviation on the n free energy profiles. In Figure , the average free energy profiles are shown together with this standard deviation represented by the error bars.

The sampling error along the collective variable can be determined by comparing the obtained free energy profiles with a fixed reference free energy profile, via a squared differences metric. The most obvious choice for the reference free energy profile corresponds to the simulations with the smallest grid spacing of 5 Å3 and the longest simulation time. Increasing the grid spacing reduces the overlap between the different stages and, hence, the sampling of the phase space in the direction of the collective variable. This is expected to reduce the precision of the obtained free energy profile.

In contrast, for the nonequilibrium methods, the sampling along the collective variable depends on the total simulation time or the total number of updates of the bias potential energy. Hence, the reference free energy profile stems from the longest simulation time (1 ns). A more complete discussion on the computation of errors is provided in Section 7 of the Supporting Information.

Rather than directly comparing the sampling errors, we introduce an error threshold and estimate the computational effort to reach this threshold. Two thresholds, 0.7 and 0.5 kJ/mol, are introduced. We compare those thresholds with the two empirically determined sampling errors, which are averaged over the volume range of interest, i.e., from 800 Å3 to 1500 Å3. The approach is inspired by the work of Bruckner and Boresch.[96]Figure summarizes the results of this efficiency analysis. The figure shows for each method the total required simulation time to reach a certain error threshold. In general, the performance of the nonequilibrium methods is better compared to the equilibrium methods in terms of the total simulation time. However, as equilibrium methods can easily be parallelized, they are much more efficient when run in a high-performance computing environment. Although the nonequilibrium methods can be parallelized via a multiple walker approach, such simulations cannot be run independently.[97] The equilibrium methods allow us to reach the lowest error threshold. The relative short simulation times of the various equilibrium runs in the TI, FEP, and US methods are partly due to the fact that the initial structures of the metastable states were rather well-known based on our earlier simulations on this system. Prior knowledge of the stable states is no prerequisite for the nonequilibrium methods, where the simulation time is only slightly affected by this initial choice.

Figure 4

Minimal computational effort required for each free energy method to meet the threshold error of 0.5 and 0.7 kJ/mol. The errors used in this plot are based on a superposition of errors due to incomplete sampling in the direction of the collective variable and incomplete sampling along the other degrees of freedom. The different blocks represent various independent and parallelizable simulations for the equilibrium free energy methods. Numerical values are provided in Table SI6 of the Supporting Information.

For the nonequilibrium methods, neither MTD nor VES is able to reach the 0.5 kJ/mol threshold given our parameter search space. The VES converges faster than MTD. A possible strategy to enhance the precision of the MTD results is to further alter the properties of the Gaussian hills, which would of course increase the overall simulation time.

For the equilibrium methods, FEP is the least efficient method. To reach the error threshold of 0.7 kJ/mol, 165 simulations (i.e., at the smallest grid distance of 5 Å3) of 50 ps each are required, adding up to a total simulation time of 8250 ps. With the parameter set employed here, FEP did not reach the error threshold of 0.5 kJ/mol. The FEP method is particularly inefficient since the grid spacing needs to be sufficiently small. Only a grid spacing of 5 Å3 allows for the construction of accurate free energy profiles, to ensure enough overlap between the various volume points. From our analysis, the umbrella sampling methods are the most efficient. More in particular, WHAM is more efficient than DHAM. However, Zhu and Hummer—correctly—noted that an important limitation of WHAM is its convergence criteria[98] (see Section 8 of the Supporting Information). In our experience, this shortcoming implies that prior knowledge of some crucial features of the final free energy profile is required to appropriately use (and trust) WHAM. An alternative approach is to use DHAM rather than WHAM to obtain free energy profiles, where no such convergence criteria are required. Whether DHAM also results in better free energy estimates for nonequilibrated simulation trajectories, as claimed by its authors,[65] was not tested in this work. The latter would require the harmonic constant of the umbrella potentials to be altered and the total simulation time to be further reduced.

For the system at hand, where the collective variable is well-defined and with prior knowledge of some specific features of the final free energy profile such as the volume range of interest, US in combination with WHAM is selected as the method of choice for the generation of the free energy profiles in conjunction with an ab initio description of the potential energy surface. For more complex transformations, where the choice of the collective variables is not so straightforward, other techniques might be necessary.

Since umbrella sampling was identified as the most efficient free energy estimation method to describe the breathing behavior of MIL-53(Al), this technique is used to obtain a more accurate free energy profile using a full quantum mechanical description of the potential energy surface at the PBE+D3(BJ) level of theory and taking into account all anharmonic contributions. This procedure was applied at three different temperatures (100, 300, and 500 K). The resulting free energy profiles are shown in Figure together with the electronic energy profile. The unit cell volumes of the cp, lp, and transition states at the different temperatures are summarized in Table SI3.

Figure 5

Free energy profile for the breathing behavior of MIL-53(Al) obtained from first-principles umbrella sampling simulations. We compare the electronic energy, including dispersions, with free energy profiles obtained at 100, 300, and 500 K.

The free energy differences between the cp and lp phases at various temperatures and using either the harmonic oscillator approximation or dynamically determined free energy profiles are shown in Table . We clearly see that the free energy difference between the lp and cp states (ΔF) is drastically reduced upon increasing the temperature. The transition barrier from the cp to the lp state (ΔF‡) also reduces, while the reverse barrier (ΔF‡) remains virtually unaffected. When comparing the free energy profile obtained at 300 K with the force field based results, it is clear that the relative stability of the cp state is too high compared to the DFT result. Moreover, based on the ab initio results in Figure , the volume of the free energy minimum corresponding to the cp state is approximately 870 Å3 which is in good agreement with experiment,[20] but a steady shift to larger volumes is visible with increasing temperatures. This volume shift is also present in the lp state; however, it is less pronounced. This volume shift accompanying a temperature increase generates a positive thermal expansion. Our DFT results yield a thermal expansion of 187.4 × 10–6 K–1 for the cp phase and 29.1 × 10–6 K–1 for the lp phase which are in line with the experimental observed thermal expansions of 120.5 × 10–6 K–1 and 18.6 × 10–6 K–1, respectively.[20]

Table 2

Free Energy Differences at Different Temperatures between the (Meta)Stable States and Transition Structures Based on First-Principles Umbrella Sampling and Free Energy Differences Obtained from Geometry Optimization at the DFT Level of Theory Complemented with Thermal Corrections in the Harmonic Oscillator Approximationa

	ΔEcp→lp‡ [kJ/mol]	ΔElp→cp‡ [kJ/mol]	ΔElp–cp [kJ/mol]
DFT
el. ener.	29.04	2.25	26.78

The electronic energies (el. ener.) correspond with electronic DFT energies including dispersion corrections. The last column (Δ) denotes the difference between the dynamic and static approaches, giving a crude indication of anharmonic effects. Force field predictions of the free energy differences at 300 K are included in the last row. These values are obtained by averaging over the various free energy methods.

A further comparative analysis with experiment reveals some other features which merit some attention.[20,24] Liu et al. observed a transformation from the lp phase to the cp phase upon decreasing the temperature to 125–150 K and a reverse transformation when increasing the temperature to 325–375 K.[20] Hence, the breathing motion is governed by a large hysteresis. The bistability observed in the free energy profiles is in agreement with this large hysteresis. The global free energy minimum for a volume of approximately 850 Å3 at 100 K is in line with the first transformation. The global free energy minimum of approximately 1460 Å3 at 500 K is in line with the second transformation. Nevertheless, when assuming a collective behavior, all free energy barriers, even those being of the order of kBT, should completely disappear for the system to undergo a transition from one state to the other.[32] At least two facts contribute to the observed mismatch. First of all, we study the breathing behavior of one unit cell; hence, we do not take into account large-scale fluctuations. Such large scale effects can be studied with the aid of multiscale thermodynamic models.[29] Second, as already indicated in Section and Ref. [40], the exact energy difference between the cp and lp phase strongly depends on the applied level of theory and dispersion scheme. To assess the accuracy of various dispersion schemes, a thorough benchmark study would be required where even higher order dispersion schemes with roots in many electron methods need to be used. This is beyond the scope of current computational resources, certainly in combination with the dynamic sampling methods used here.

However, within the given DFT scheme, anharmonic corrections appear to exert a large influence on the obtained free energy differences. For example, at 500 K, the dynamic methods predict almost equally stable cp and lp phases, whereas simple harmonic oscillator based approximations of the free energy still predict the cp phase as the most stable structure at this elevated temperature. In general, the stability between the two phases is determined by two factors: the dispersion interactions, particularly between the linkers which stabilize the cp phase, and second the entropy contributions, which have the tendency to stabilize the lp phase at higher temperatures. The latter seem to be crucially dependent on the specific method to determine the thermal corrections. The lp state is much more stable using the dynamic sampling methods that properly include the anharmonic corrections. Rather than using NMA, intermediate methods such as quasi-harmonic approximation (QHA) can be used. Via constrained optimization, anharmonic effects can be partially taken into account. At the force field level, we find a free energy difference of 23.63 kJ/mol using QHA with respect to 30.48 kJ/mol using NMA (see Figure SI3 of the Supporting Information).

Conclusions

In this work, we studied the construction of free energy profiles for breathing metal–organic frameworks. More in particular, we investigated the distinct breathing behavior of the MIL-53(Al) framework, which transforms from a closed pore to a large pore structure characterized by a large change in unit cell volume under the influence of temperature and pressure. Estimating free energy profiles at finite temperature is a challenging task as a good estimate of both the potential energy surface on the one hand and the thermal corrections and entropy contributions on the other hand are required. To account for both aspects, a two-step procedure was adopted, where first the most efficient free energy method at a force field level is determined and second the selected method was used to construct a more accurate free energy profile at a high level potential energy description.

Several free energy estimation methods were assessed for their ability to generate the free energy profile of the MIL-53(Al) framework taking into account the full dynamic behavior of the material. Five free energy methods were tested, namely, free energy perturbation, thermodynamic integration, umbrella sampling, metadynamics, and variationally enhanced sampling. All these methods succeed in predicting the bistable behavior. With the force field description of the potential energy surface, the free energy difference between the closed pore and large pore state amounts to 26.57 kJ/mol at room temperature, corresponding with an electronic energy difference of 42 kJ/mol, with the closed pore state being the most stable. This free energy difference is in principle too large to match experimental observations. A large influence may be expected from the particular method to estimate the dispersion corrections. Indeed, DFT predicts energy differences between the two phases of 26.78 and 25.48 kJ/mol at the PBE+D3(BJ) and PBE-MBD schemes, respectively. However, as both dispersion corrections and anharmonicity contribute to the free energy profile, one should ideally use a dynamic method in conjunction with a DFT based description of the potential energy surface. To achieve this goal, we carefully assessed the efficiency of the five free energy methods at the force field level. For the given system, umbrella sampling was selected as the most efficient free energy method.

Finally, umbrella sampling was used to predict the free energy profile at three temperatures (100, 300, and 500 K) based on a DFT PBE+D3(BJ) description of the potential energy surface. A large number of parallel simulations in the (N, P, σ = 0, T) ensemble were run to achieve this goal. The analysis shows that anharmonic corrections are very important to estimate the free energy profile of MIL-53(Al). Upon increasing the temperature, entropy contributions become more important, yielding a more stable large pore phase compared to methods based on the harmonic oscillator approximation. At 100, 300, and 500 K, free energy differences between the large and closed pore phase of 18.67 kJ/mol, 9.74 kJ/mol, and −0.67 kJ/mol are predicted.