Configurational Entropy Driven High‐Pressure Behaviour of a Flexible Metal–Organic Framework (MOF)

Abstract Flexible metal–organic frameworks (MOFs) show large structural flexibility as a function of temperature or (gas)pressure variation, a fascinating property of high technological and scientific relevance. The targeted design of flexible MOFs demands control over the macroscopic thermodynamics as determined by microscopic chemical interactions and remains an open challenge. Herein we apply high‐pressure powder X‐ray diffraction and molecular dynamics simulations to gain insight into the microscopic chemical factors that determine the high‐pressure macroscopic thermodynamics of two flexible pillared‐layer MOFs. For the first time we identify configurational entropy that originates from side‐chain modifications of the linker as the key factor determining the thermodynamics in a flexible MOF. The study shows that configurational entropy is an important yet largely overlooked parameter, providing an intriguing perspective of how to chemically access the underlying free energy landscape in MOFs.


Pawley Profile Fit Analysis
Pawley profile fits of the data were performed using TOPAS V6. [3] The standard deviations (σ) of lattice parameters and volumes were obtained and the error of the pressure was estimated as p = ± 0.0020 GPa. Based on the sharp diffraction peaks the presence of diffraction domain sizes in the nanoregime can be excluded.
Cu2bdc2dabco was indexed in the space group P4/mmm (No. 123) with two minor additional reflections at 2θ = 2.63° and 2.81°, which were observed in previous studies [5] and therefore fitted manually. All cell parameters and their standard deviations are listed in Table S1.
Cu2(DB-bdc)2dabco crystallises like its parent MOF in a tetragonal space group. Although Schwedler et al. indexed the Zn analogue in the space group P4/ncc (No. 130), [6] the Pawley fits using Topas V6 point out that for the Cu analogue the space group P4/n (No. 85) suits best for the lp form. The space groups differ in their diffraction conditions (00l : l = 2n for P4/ncc but not for P4/n) [7] which leads to a doubling of the cell size and reflections that are not observed in the experimental diffraction pattern (e.g. 2θ = 3.8, 5.6, 5.9, 6.4° etc.). The phase transition comes with a symmetry lowering to the monoclinic space group P2/c (No. 13). Similar as reported by Henke et al. for the flexible Zn2(BME-bdc)2dabco, [8] we observe that the first reflection (corresponding to the (110) reflection) shifts from 2θ = 2.25° to 2θ = 2.65° and the (200) reflection from 2θ = 3.18° to 2θ = 2.63° which is attributed to an elongation of cell parameter a and contraction of cell parameter b. Henke et al. also reported that the (001) reflection at 2θ = 2.53° does not shift, which is different for our material. Here, the (001) reflection is not seen for the np phase because of the reflection conditions of the space group P2/c which are h0l: l = 2n and 00l: l = 2n. [7] Instead the (010) reflection is seen at 2θ = 2.30°. Note, that it cannot be deduced from the PXRDs as a bulk method, if both phases are observed simultaneously in one particle or exist next to each other in different particles. From p = 0.15 GPa both phases are present and therefore a biphasic fit is performed. All cell parameters and their standard deviations are listed in Table S2 and Table S3.        Table S2 because the data were fitted with a biphasic fit from p = 0.15 -0.425 GPa. Note that this space group embeds features of two pores per unit cell. Thus, to compare the volumes reported here with the volumes mentioned in the main text, these have to be divided by a factor of two.

NiDMG as Reference Material
To verify the setup and the pressure applied, Nickel dimethylglyoxime (NiDMG) was used as a reference material because of its known pressure behaviour (see Refs. [21,22] ). The capillary was prepared in the same way as the MOFs except that is was not done in the glovebox. 22 HPPXRD patterns were collected between p = ambient and 0.40 GPa including one after releasing the pressure. The step size was ∆p = 0.02 GPa. NiDMG was successfully indexed in the space group Ibam (No. 72) and with the cell parameters a = 16.58, b = 10.44, c = 6.46 Å and V = 1120.43 Å 3 as reported by Refs. [21,23,22] . A complete list with all cell parameters and their standard deviations can be found in Table S5. The bulk modulus was K(NiDMG) = 8.31 ± 0.08 GPa which is in accordance with the bulk modulus reported by Takeda

Computational Details
The simulation cells used in all simulations contain four formula units and hence four pores of the respective MOFs amounting to 216 atoms for Cu2(bdc)2dabco and 424 atoms for Cu2(DB-bdc)2dabco. All results shown herein are given per formula unit/pore. Details on how the structures were prepared can be found in the section below. Pressure ramp NPT simulations were performed with the LAMMPS implementation of the anisotropic MTTK barostat (from https://github.com/stevenvdb/lammps/tree/newbarostat) using LAMMPS's native pressure ramp functionalities. A thermostat relaxation time of t = 100 fs and a barostat relaxation time of t = 1 ps was used for all runs. For each ps timeframe a structure was written. For the series of NV(σa = 0)T simulations the average lp and np structure cell parameters were taken from the pressure ramp runs and the volume was interpolated linearly between these forms using 50 volume windows. For each target volume the structure that differs in volume as little as possible was chosen and scaled to the target volume.
The NV(σa = 0)T simulations were performed with the MTTK barostat using the same settings as above but utilizing a volume constraint. Initial structures were equilibrated for t = 200 ps after which the stress tensor was recorded in t = 1 ns production runs. PXRD patterns were calculated every t = 10 ps using FOX [24] with a Lorentzian line shape width of 0.001 ° and the wavelength of Sn kα at λ = 0.4246 nm also used in the experiment. The p(V) equation of state was fitted for Cu2(DB-bdc)2dabco and interpolated using cubic splines for the unfunctionalised MOF. These continuous functions were used for any further analysis of the p(V) EoS described below. The bulk modulus K was computed from the EoS at the roots of p(V) where the slope is negative by numerically estimating the derivative ⁄ from the p(V) fit and applying the formula The Helmholtz free energy A was obtained via numerically integrating the p(V) EoS starting from the root of p(V) corresponding to the large pore by using From the average of the internal Energy U(V) of each simulation, the entropic energy contribution TS(V) and the entropy S(V) was computed via the fundamental relation U = A-TS. All thermodynamic quantities were arbitrarily shifted so the minimum corresponds to the origin, since their absolute differences do not carry any meaning without having sampled a proper path that connects them. The vibrational entropy contributions in the harmonic approximation were computed as follows: After a tight optimization of the initial structural models, the Hessian matrix of second derivatives was computed by finite differencing of the forces using MOF-FF in LAMMPS. Imaginary frequencies were assured to be absent. The entropy at T = 300 K was computed from the Hessian using phonopy [25] v2.1.4. using a 12x12x12 q-point mesh, which we assured to be reasonably well converged. The calculations for Cu2(DB-bdc)2dabco were performed for a single snapshot for each Group-1 structure at a volume of V = 937.5 Å 3 (np) and V = 1125 Å 3 (lp). The averaged entropy penalty term -TΔS obtained in this way amounts to 21.1 kJ•mol -1 with a standard deviation of 2.9 kJ•mol -1 . All simulation data are available upon reasonable request. Run scripts and the analysis routines are accessible at the CMC groups software repository at https://github.com/cmc-rub/supporting_data/tree/master/XX-Keupp-XX_2020 Whereas the model preparation of Cu2(bdc)2dabco is relatively straight forward, the disorder present in longer side chains poses a problem in the simulation of the free energy profiles for these systems: We have recently pointed out that the simulations of functionalised MOFs as the herein presented Cu2(DB-bdc)2dabco strongly depend on the initial linker configuration, since within the accessible simulation timescales of several nanoseconds a linker flip is observed only very rarely. Therefore, different initial configurations do not transform into each other, but the result we obtain is that of an infinitely large crystallite composed solely of this single configuration. Initial structural models derived from experiment are not readily available as the single crystal data reported by Henke et al. [1,26] cannot uniquely be transformed into an initial model to use for the simulations because the side chains are not resolved (even at T = 87 K) and the positions of the oxo-moieties are disordered. Thus, there are many different configurations that could be valid structures based on the comparison to the XRD data. We therefore constructed many trial initial structural models from the blueprint shown in Figure S15 using the weaver code [27] followed by a lattice optimization step: 32 individual starting structures with random phenyl orientations (i.e. random torsion angles as indicated by the blue and magenta colours in Figure S15) were generated, optimized and used to run the NPTs from which the p(V) EoS is computed.
The evaluation of the trial structures of Cu2(DB-bdc)2dabco is done by comparing PXRD patterns computed along and averaged over the NV(σa = 0)T trajectories at a cell volume of V = 937.5 Å 3 . Figure S16 shows all 32 PXRDs compared to the experimental highpressure pattern at p = 0.4 GPa. To assess the similarity between two PXRD patterns, we used the Rietveld fit inspired value function:  Figure S16. A similarity above f = 1500 is considered a full match and a value below f = 1000 is considered a full mismatch. In between we linearly interpolated the two colours, but for a zero-one decision, we used a threshold of f = 1300. The averages of Group-1 and Group-2 are thus calculated from all structures featuring an f smaller or larger than 1300, respectively. The averaged PXRD patterns are shown in Figure 4 in the main text. Interestingly, after this classification into Group-1 and Group-2, where Group-1 are the structures with matching PXRD patterns and Group-2 are those that do not match, their distinct behaviour can also be observed in the PXRD patterns along the pressure ramp trajectories: Whereas the Group-1 structures show a step as observed in the experiment, the Group-2 structures do not. Instead one can follow their PXRD pattern shifting continuously as the pressure increases.

Cell Parameters
The cell parameters were averaged for each NV(σa = 0)T simulation. Figure S19 shows the cell lengths and Figure S20 the cell angles, respectively. All Group-1 structures compare very well with the experimental cell angles, whereas in the Group-2 structures qualitative differences are found. In general Group-2 type structures change their aspect ratio to a lesser extent or not at all, whereas all Group-1 structures have a very similar cell angles vs volume profile. Group-1 structures differ in itself in the point at which the cell lengths a and b start to differ from each other. In any case this point is close to the lp cell lengths.
Cell angles as a function of cell volume differ more from the experimental values. In general, the amplitude by which cell angles differ from 90 ° is often higher and even the lp form volumes of some of the Group-1 structures adopt cell angles different from 90 °. We attribute the flexibility of the cell angles to the respective linker configurations of the MOF combined with our small simulation cell comprised of only four formula units of the MOF material. Neighbouring pore environments will allow for certain distortions of a given pore and restrict others, an effect we can only capture in larger simulation cells. Probably in experimentally derived cell angles by means of PXRD fits, these effects average out towards a small deviation in the cell angles of only 1.5 °.   From the p(V) EoS (cf. Figure 3) ∆A was computed by numerical integration, and from the average of the total energy, the internal Energy ∆U was obtained. Since the free energy difference of two different structures cannot be inferred in this way, each individual structure was shifted to its own origin. For ∆U it is different, since the very same force field terms give rise to a comparable total energy. Therefore, all structures were shifted to the global minimum of the global minimum structure. The plots above show the data derived from polynomial fits of p(V) of the raw data points every V = 12.5 Å 3 . Experimentally it is observed that Cu2(DB-bdc)2dabco at ambient conditions remains at its lp form even after removal of solvent or guest molecules. This is properly predicted by the global free energy minimum of all structures at V = 1125 Å 3 . As seen from the p(V) plots, a volume step at elevated pressures between p = 0.1 and 0.3 GPa is present for some of the simulated structures, but none of them is a metastable phase at ambient pressure. This is also in accordance with the fact, that in experiment after compression, the phase transition is entirely reversible. Interestingly, the internal energy of the Group-1 structure models has its global minimum at a volume corresponding to a np form (approx. V = 950 Å 3 ). It is the entropy that adds an average of 42 kJ mol -1 per simulation cell to the free energy to render the np form a high-pressure phase that cannot be trapped in a metastable state at ambient temperatures. Also note that the internal energies of the lp form are almost indistinguishable when comparing Group-1 and Group-2, which may hint towards a strong disorder in this system, since Cu2(DBbdc)2dabco is synthesized in the lp form at elevated temperature.

Structural Details via Collective Variables
In order to analyse the individual structures on an atomistic level, a selection of collective variables shown in Figure S22 was computed along all trajectories of all NV(σa = 0)T simulations. The quest is to find those collective variables, where the two sets of structures that do (Group-1) and do not (Group-2) fit the np PXRD pattern differ. The histograms of the collective variables shown in Figure S23 were collected along the NV(σa = 0)T trajectories of the given volume. The distributions of the lp form (at V = 1125 Å 3 ) collective variables are almost identical, with only slight variations in the COCC, OCCC and CCCC profiles that point towards the side chains being a little more elongated and stiff in case of the Group-1 structures. For the np forms this trend is amplified as the peaks for the Group-1 structures get sharper whereas the peaks for the Group-2 structures get shallower, indicating a slightly larger degree of coiling for the Group-2 structures. In contrast to that, the collective variables OMMO, MOCC and ph_CCCO are very different in the np forms comparing the Group-1 and the Group-2 structures with the largest difference residing in the distribution of the OMMO collective variable. The OMMO torsion angle describes a rotation of two paddlewheel units stacked in z-direction in a counter-rotating fashion, leading to the loss of the alignment of the linkers as one looks at the structure along the z-direction. We believe that this effect is what stiffens the structures throughout the simulations. This change is accompanied with a distorted paddlewheel unit as the MOCC torsion shows a larger change and hence a wider distribution for the Group-2 structure in the np forms. Finally, for the ph_CCCO torsion there is even a qualitative difference comparing the two sets of structures: In the np form, the Group-1 structures show four maxima in the 0 ° region, whereas the Group-2 structures show only two peaks, which are even narrowed down and shifted towards 0 ° during the lp to np transition. Whereas all other collective variables statistics are very similar in the lp phase for the two sets of structures, the ph_CCCO torsion profile is very different even in the lp phase, indicating that it is the main reason for the apparent differences in the phase transition behaviours simulated herein. Note that this collective variable was sampled uniformly for each of the eight phenyl moieties of the 32 investigated starting structures in the construction process.

Variability of the Linker Orientations of the Cu2(DB-bdc)2dabco Simulations
The entire system is comprised of eight different sites where the functionalized bdc linker is inserted. For each of these different sites, a uniform random orientation angle for the linker is drawn and the linkers is rotated accordingly. The number of distinct configurations depends on the question how these sampled orientations change during the simulation. If we assume that there are four distinct stable orientations, it will amount to 4^8= 65.536 different configurations. The data to analyse this is readily available as the ph_CCCO collective variables were computed for all trajectories. Below is a figure displaying the individual contributions of all 16 torsion profiles that give rise to the total histogram also shown in the Figure S23. The data presented here is at V = 1125 Å 3 for the Group-1 structures of Cu2(DB-bdc)2dabco. We observe the general trend that this torsion deviates from the 0 / 180°angle, but it is not a clean cut, but takes into account both this in-plane configuration and the ± 30° deviation from it within a single simulated structure. Figure S24. Histograms of the ph_CCCO collective variable for the Group-1 type structure. The coloured lines denote one single collective variable, the black lines correspond to the histogram of all individual collective variables. The latter have been scaled by a factor of five for visibility.