Thermodynamic Insight in the High-Pressure Behavior of UiO-66: Effect of Linker Defects and Linker Expansion.

In this Article, we present a molecular-level understanding of the experimentally observed loss of crystallinity in UiO-66-type metal-organic frameworks, including the pristine UiO-66 to -68 as well as defect-containing UiO-66 materials, under the influence of external pressure. This goal is achieved by constructing pressure-versus-volume profiles at finite temperatures using a thermodynamic approach relying on ab initio derived force fields. On the atomic level, the phenomenon is reflected in a sudden drop in the number of symmetry operators for the crystallographic unit cell because of the disordered displacement of the organic linkers with respect to the inorganic bricks. For the defect-containing samples, a reduced mechanical stability is observed, however, critically depending on the distribution of these defects throughout the material, hence demonstrating the importance of judiciously characterizing defects in these materials.

Introduction

Metal–organic frameworks (MOFs) or porous coordination polymers (PCPs) continue to receive abundant attention since their first synthesis about two decades ago.[1−3] These nanoporous materials consist of inorganic moieties interconnected by organic linkers to form ordered crystal lattices.[4−6] Their nanoporous structure makes them very tractable for application domains such as heterogeneous catalysis,[7,8] controlled drug release,[9] and gas storage and separation.[10] The application field of these materials is broadened even further thanks to the concept of isoreticular synthesis,[11,12] allowing to incorporate specific functional groups in the ligand while retaining the topology of the framework structure, hence tuning the material for specific applications.[13−15]

Despite the thriving interest in MOFs, exploration of their behavior under pressure has emerged only recently, and, as a consequence, is not yet well understood.[16−19] In this context, the Cheetham group has been very active in studying mechanical properties in both porous and dense MOFs, and their progress was summarized in two landmark reviews.[20,21] The main focus in these initial studies went to zeolitic imidazolate frameworks (ZIFs)[22] and hybrid perovskite structures,[23] but was very recently extended to also cover the wide range of MOF properties in the presence of defects.[24,25] In many of these papers, the importance of combining experimental and computational techniques to reveal structure–property relationships on a molecular scale was emphasized.[26,27] Also for other materials, a thorough understanding of their mechanical behavior is often a prerequisite for promoting them at the application level.[18] Indeed, MOF powders need to undergo sintering and pelletization procedures before their industrial use.[28] It was recently shown that these procedures may irreversibly change the MOF’s structure.[29−31]

A class of materials that received considerable attention in recent literature are the UiO-66-type materials, which are composed of inorganic Zr6(μ3-O)4(μ3-OH)4 bricks connected through ditopic organic ligands (see Figure ) and were first synthesized by Lillerud and co-workers.[32] Within the MOF field, they exhibit an exceptionally high thermal stability[33,34] and retain their crystal structure under high pressures[28,30,35] as well as in relatively harsh acidic environments.[33,36,37] This exceptional stability may be traced back to the inherent composition of the frameworks, since each zirconium atom is 8-fold coordinated by oxygens and each inorganic brick is 12-fold coordinated by organic linkers in the defect-free materials.

Figure 1

fcu–a topology of the UiO-66 to -68 nets, with indication of the octahedral cage (blue sphere) in the conventional unit cell. Also indicated are the inorganic building blocks (left) and organic ligands (right) needed to describe the 11 structures discussed in this work. Zirconium atoms are shown in cyan, oxygen atoms in red, carbon atoms in gray, and hydrogen atoms in white.

Recently, the mechanical behavior under high pressure of nonfunctionalized UiO-66 and its amine-functionalized analogue was further explored by measuring high-pressure powder X-ray diffraction (PXRD) patterns for pressures up to 3.5 GPa.[35] A significant broadening of the Bragg peaks under pressure was observed, pointing toward a loss of crystallinity in these materials. This pressure-induced decrease in crystallinity was identified as a reversible phenomenon.[35] In this Article we give a molecular-level insight into the observed loss of order under external pressure at room temperature. We do not only investigate the defect-free UiO-66 material, but also explore the influence of expanding the length of the organic linker and of structurally embedded defects on the local order of the material when subject to elevated pressures. The extended materials, UiO-67 and UiO-68, are topologically identical to UiO-66, but are composed of longer biphenyl-4,4′-dicarboxylate (BPDC) and p-terphenyl-4,4″-dicarboxylate (TPDC) linkers, respectively, instead of the benzene-1,4-dicarboxylate (BDC) linker in UiO-66. Intuitively, one might expect a softening of the materials when linkers with more internal degrees of freedom are introduced. Recently, some of the present authors indeed observed flexibility and a loss of crystallinity in the otherwise rigid UiO-66 by introducing conformational flexibility in the organic linkers.[38]

UiO-66-type materials piqued the interest of many scientists since the discovery that their properties could be significantly altered by structurally embedded linker defects. These defects, corresponding to the removal of organic linkers, reduce the inorganic coordination number, as shown via gravimetric characterization,[39,40] X-ray diffraction,[33,41,42] and neutron power diffraction,[43] and verified via density functional theory (DFT) calculations.[44,45] Recently, Goodwin and co-workers showed that correlations between defects can be introduced and controlled, yielding nanoscale defect structures.[46] The presence of these defects may alter the catalytic properties,[47,48] thermal stability,[49] proton conductivity,[50,51] and adsorption behavior[43,52−55] of the host material. While recent work also indicates a decrease in bulk modulus and hence robustness upon the introduction of defects,[54,56] it remains to be investigated how the precise molecular nature of these defects alters the pressure-induced disorder in UiO-66-type materials.[25]

Herein, we tackle these questions by extending a recent, force field based procedure to construct pressure-versus-volume curves at finite temperature.[57] In this procedure, molecular dynamics (MD) simulations in designated thermodynamic ensembles sample the material at various fixed volumes, yet allowing the cell shape to fluctuate. As a result, pressure-versus-volume profiles can be constructed, and the different (meta)stable phases of the material can be extracted. To simulate the material under the high pressures assumed in this work, the availability of accurate and robust force fields is a prerequisite. Some of the present authors recently developed the QuickFF force field protocol to generate first-principles based force fields for MOFs.[58] This protocol was shown to be successful in describing structural features and breathing phenomena in for instance MIL-47(V),[59,60] the phenyl and fumarate versions of MIL-53(Al),[58,61] and MOF-5.[58] For the UiO-66-type materials under study here, the construction of these force fields is particularly challenging since they need to capture the variation in coordination number as well as the variation of linkers. We show herein that our force field protocol is able to capture these features and succeeds in predicting the mechanical instability of the material, coinciding with its experimentally observed amorphization. A critical dependence of this behavior on the number and position of defects and the length of the linkers is observed. To investigate the effect of the molecular structure of the defects on the pressure-induced instability, we propose a classification of defects differing in structural short-range order. Our simulations yield unprecedented molecular-level insight in the pressure-induced behavior of UiO-66-type materials.

Methodology

Materials

All materials discussed in this Article are defined by the conventional unit cell containing four inorganic Zr6(μ3-O)4(μ3-OH)4 bricks (see Figure ). To determine the effect of linker vacancies and expansion of the linkers on the mechanical properties of UiO-66, 11 materials will be considered: three pristine materials and eight materials containing linker defects. A first set of materials is obtained by elongating the benzene-1,4-dicarboxylate (BDC) ligand in UiO-66 with phenyl moieties, obtaining biphenyl-4,4′-dicarboxylate (BPDC) and p-terphenyl-4,4″-dicarboxylate (TPDC), giving rise to UiO-67 and -68, respectively, as depicted in the right pane of Figure . A second set of materials, shown in Figure , is obtained by removing one or two BDC ligands from the pristine UiO-66 unit cell. The type 0 structure is retrieved by withdrawing one of the 24 BDC ligands present in the pristine UiO-66 unit cell. As a result, two of the four inorganic bricks are now only 11-fold coordinated, leading to an average coordination number ⟨CN⟩ of 11.5. To create an average coordination number of 11, an additional ligand needs to be removed. While 23 possibilities exist to remove the second ligand, some of these structures are physically equivalent. These 23 possibilities can be classified in seven distinct classes on the basis of the distance between the two deleted ligands, their relative orientation, and the coordination number of the four inorganic bricks in the unit cell, as outlined in Section S1. For each of these seven classes, which are referred to as type 1 to type 7, a representative structure is shown in the top pane of Figure . While a careful tuning of the synthetic parameters, such as temperature and modulator composition, can reduce the number of linker vacancies, most synthetic procedures yield a coordination number of about 11–11.5.[40,62]

Figure 2

Top: Defect-free structure (with either BDC, BPDC, or TPDC ligands) and the eight types of defects (with BDC ligands) considered in this work, in which the linker vacancies are indicated in red (dotted lines represent periodic images of the solid lines). Bottom: Detail of the Zr6 octahedra present in the top pane, with indication of the positions where ligands are missing (red) and the resulting coordination number (CN). The zirconium atoms are color coded based on their coordination number, and correspond, from light to dark, with a coordination number of 8, 7, and 6.

Force Field Generation

All force fields in this work are generated using the in-house developed QuickFF procedure and are composed of two contributions: the covalent and the noncovalent contribution.[58] On the one hand, the covalent (cov) terms, which mimic the chemical bonds, are analytical expressions in terms of internal coordinates, such as bond lengths, bend angles, dihedral angles, or torsions, and out-of-plane distances (oopd). On the other hand, the noncovalent (noncov) terms model the long-range interactions between nonbonded atoms, and consist of an electrostatic (EI) and a van der Waals (vdW) part. Hence, the potential energy surface , function of the nuclear coordinates, is approximated by the analytical expressionAll force field contributions, except the van der Waals part, are fitted to first-principles data on representative cluster model systems.

Choice of the Cluster Model Systems

Isolated model systems are frequently used to parametrize all-atom force fields via first-principles data.[63,64] To describe the 11 materials in this study, and hence fully take into account the varying coordination number of the inorganic bricks, six types of isolated model systems need to be considered for the inorganic brick. These model systems are constructed by terminating each of the CO2 groups of the six inorganic bricks (see Figure (left) and Figure (bottom)) with phenyl rings, mimicking the environment present in the periodic system. A similar formate-terminated model system was used by the authors of MOF-FF to derive a force field for the pristine UiO-66 material,[64] and the validity of the model systems was assessed in various studies.[44,65,66] Note that these inorganic bricks are hydroxylated. For other applications, such as catalysis or grafting of specific moieties,[67] pretreatment may lead to a certain degree of dehydroxylation.[42,44] For the organic ligands (see right pane of Figure and Figure ), the isolated model systems are obtained by terminating the ligands with an inorganic Zr6O4(OH)4(HCO2)12 brick at each side.

Figure 3

Unique atom types for the pristine UiO-66 to -68 materials.

The initial isolated cluster models were extracted from DFT optimized periodic structures.[68] These isolated models were subsequently optimized with DFT, using the B3LYP functional and keeping the terminating hydrogens fixed. For every structure, the geometry, Hessian and electron density are calculated as outlined in the next paragraphs, and used as input for QuickFF to estimate the unknown force field parameters.[58] For each of the 11 materials, a separate, periodic force field can then be constructed by combining the force field parameters of the relevant isolated cluster models (see Figure ). Note that, while the type 3, type 5, and type 7 defects correspond to the same combination of inorganic bricks, their relative orientation is different, leading to a material with a different organization at the molecular level. The procedure for determining the force field parameters for each of the above-mentioned structures will be outlined below. A more detailed description can be found in Section S2.1.

Determination of the Force Field Parameters

Electrostatic Interactions

The electrostatic interactions between two atoms are described by the Coulomb interaction between spherical Gaussian densities, for which the charges are derived from the DFT electron density of each model system, while the charge radii are based on the fitting procedure of Chen and Martinéz.[69] No exclusion rules are taken into account, and charges from the cluster calculations are transferred to the periodic structure via bond-charge increments,[70] similar to the approach in ref (63).

van der Waals Interactions

The van der Waals interactions are modeled based on the two-parameter MM3 Buckingham potential,[71] for which the two parameters σ and ε are determined from the atomic values via the Lorentz–Berthelot mixing rules. The atomic parameters σ and ε are taken from refs (71) and (72), adopting the 1–2 and 1–3 exclusion rules for bonded pairs from the MM3 rules.

Covalent Interactions

The unknown parameters in the covalent contributions, that is, the force constants K, the rest values q0 and the multiplicities m, are estimated with QuickFF, yielding a force field contribution for each of the bonds, bends, out-of-plane distances and dihedral angles. A harmonic potential is used for all bond, bend and out-of-plane terms as a function of respectively the bond length r, the bend angle θ and the out-of-plane distance d. The torsional terms, except the two terms discussed in more detail below, are described with a simple cosine term with multiplicity mϕ = 2 depending on the dihedral angle.

As already anticipated, these covalent terms are in some specific cases not sufficient to adequately describe the system. A first example is the C3–C4–C4–C3 dihedral pattern in the BPDC linker of UiO-67, describing the torsion angle of one phenyl ring with respect to the other (see Figure ). For an arbitrary ϕ0, it is not always possible to obtain energy minima at ±ϕ0 and ±(π – ϕ0) with a cosine of period 2π/m (m ∈ ). Therefore, the force field should be refined by including material-specific terms during the fitting procedure. In this specific case, the following term is proposed to adequately describe the C3–C4–C4–C3 dihedral anglewhere ϕ0 = 28.423° is obtained from the DFT optimized cluster geometry.

A second example is the C3–C4–C5–C6 dihedral pattern present in the TPDC linker of UiO-68, determining the dihedral angle between the middle and the outer phenyl rings (see Figure ). Here, an extra DFT scan is performed on the TPDC model system to obtain an accurate description of this dihedral motion. A new term is fitted to the DFT profile, excluding the electrostatic interactionswith ϕ0(1) = 40° and ϕ0(2) = 110°. As shown in Figure S5, this new term nicely captures the simulated dihedral profile.

Thermodynamic Ensembles

To fully capture the response of UiO-66-type materials upon application of an isotropic pressure P, MD simulations are performed in two ensembles. In this Article, we will adopt the terminology introduced in ref (57), where the (N, P, σ = 0, T) and (N, V, σ = 0, T) ensembles were first introduced. This terminology is based on the observation that some degrees of freedom present in the cell tensor h and the stress tensor σ can be conveniently separated. When using these two ensembles, different degrees of freedom are kept fixed. On the one hand, in the (N, P, σ = 0, T) ensemble, the number of particles N is fixed, while the internal pressure P, the internal deviatoric stress σ and the internal temperature T are controlled. This ensemble is also often called the flexible NPT ensemble, since it allows both the cell volume V and cell shape h0 to fluctuate. On the other hand, in the (N, V, σ = 0, T) ensemble, the volume V is kept fixed instead of the average pressure. In this type of simulation, the cell shape h0 will fluctuate such that, on average, the internal deviatoric stress σ equals the applied deviatoric stress σ = 0. While the (N, P, σ = 0, T) ensemble can be efficiently used to simulate the response of a material under an isotropic pressure P, the (N, V, σ = 0, T) ensemble allows one to determine the pressure the material can withstand at a given volume V. As a consequence, the pressure-versus-volume behavior of the material is unveiled, providing information such as its equilibrium volume, bulk modulus and, if present, the pressure at which a phase transition takes place. Note that, following this approach, the pressure-versus-volume behavior is retrieved by fixing the volume instead of controlling the pressure. The latter approach has been applied in refs (73) and (74) to investigate the mechanical stability of ZIFs. While both approaches are expected to yield the same information for large simulation cells, it was outlined in ref (57) that fluctuations in the instantaneous pressure may induce phase transitions at artificially low pressures. While these fluctuations can be controlled by employing sufficiently large simulation cells, their effect may be mitigated altogether by employing our fixed volume approach. Moreover, note that, while both methods may unveil the pressure-induced mechanical instability of the material, they do not provide direct information on the amorphous phase, since nonreactive force fields do not capture the bond-breaking behavior associated with amorphization.[73]

Pressure-versus-Volume Behavior

As outlined above, the (N, V, σ = 0, T) ensemble can be used to obtain the pressure-versus-volume behavior of UiO-66 and its analogues. To do so, first a volume grid is defined on the interval [V, V] with a grid spacing ΔV, where V, V, and ΔV depend on the material under study. Then, (N, P, σ = 0, T) simulations are carried out starting from a sufficiently high volume and at a sufficiently high pressure, such that during the simulation the instantaneous volume V reaches each of the predefined volume grid points within a given threshold. Each of these configurations is then used to initialize a separate (N, V, σ = 0, T) simulation at a volume V corresponding to a volume grid point.

For each of these (N, V, σ = 0, T) simulations, the average isotropic pressure ⟨P(V)⟩ is calculated. When assuming mechanical equilibrium, this internal isotropic pressure counteracts the external pressure the material can withstand at the given volume V. When repeating this for each of the volume grid points in the interval [V, V], the pressure-versus-volume curve is obtained in this range, yielding information on the equilibrium volume at each pressure as well as on the bulk modulus K = −V ∂ ⟨P⟩/∂V of the material.[75] Here, ⟨P⟩ as a function of the volume is first fitted by means of a seventh-order polynomial, where the order of the polynomial is chosen such that no overfitting issues arise, that is, such that the Vandermonde matrix is not rank-deficient. As a result, the derivative above can be carried out analytically, suppressing any noise.

Computational Details

The building blocks for the UiO-66-type materials were assembled using Zeobuilder.[76] These isolated clusters were then optimized via DFT calculations in Gaussian 09,[77] using the B3LYP exchange-correlation functional[78−81] and keeping the outer hydrogens fixed. The 6-311G(d,p) Pople basis set was used[82] for all atoms, except for zirconium, which was described using the LANL2DZ basis set, including an effective core potential.[83] A vibrational frequency analysis ensured that the optimized structures correspond to minima on the potential energy surface.

The covalent force field parameters were obtained with QuickFF.[58] For the noncovalent force field parameters, atomic charges were obtained according to the Minimal Basis Iterative Stockholder (MBIS) partitioning scheme, a new iterative variant of the Hirshfeld atom-in-molecules (AIM) scheme recently developed by Verstraelen et al.[84] and implemented in HORTON.[85] The all-electron density necessary for the MBIS scheme was obtained via a GPAW calculation.[86−88] Finally, the van der Waals interactions are based on the MM3 model of Allinger et al.[89]

The MD simulations reported in this work were all carried out in either the (N, P, σ = 0, T) or the (N, V, σ = 0, T) ensemble using Yaff, a freely available in-house developed software package,[90] using the conventional unit cell with four inorganic bricks. The influence of using a larger supercell on the observed mechanical behavior was also tested, by constructing the pressure-versus-volume curve for the defect-free UiO-66 using a 2 × 2 × 2 supercell containing 3648 atoms. In this case, the long-range (electrostatic and van der Waals) interactions were calculated with LAMMPS for the sake of computational efficiency.[91] For this particular material, it was shown that the conventional unit cell is sufficiently large to discuss its mechanical stability. More information on this assessment can be found in Section S5.

The equations of motion were updated via a Verlet scheme, with a time step of 0.5 fs to ensure energy conservation. The electrostatic interactions were efficiently calculated 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.32 Å–1.[92] The van der Waals interactions were also calculated with a smooth cutoff at 15 Å. The temperature during these simulations, fixed at 300 K, was controlled via a single Nosé–Hoover chain consisting of three beads and with a relaxation time of 100 fs.[93−96] This thermostat was coupled to both the particles and the barostat. To control the pressure, a Martyna–Tobias–Tuckerman–Klein (MTTK) barostat was employed with a relaxation time of 1000 fs.[97,98] As shown previously, this combination of relaxation times ensures a complete yet efficient sampling of the accessible phase space.[57] All simulations were equilibrated for 50 ps, followed by a 500 ps production run.

VMD was used to visualize the MD trajectories,[99] while the pore size distribution was calculated with Zeo++.[100−102] Symmetry analyses were performed with PLATON, using the default error margin.[103]

Results and Discussion

Force Field Validation

Before applying the force fields to study the pressure-induced behavior of the UiO-66-type materials, their validity and robustness at finite temperature need to be assessed. Both geometry optimizations and (N, P, σ = 0, T) simulations at 300 K were carried out. As shown in Section S2.4, the force fields reproduce the internal coordinates of the UiO-66 to -68 materials very well when compared to the DFT cluster calculations. Comparison with periodic DFT calculations, including dispersion interactions, shows that the derived force fields slightly overestimate the unit cell lengths by 0.7–1.0%.

The true interest lies in the validation of experimental data at finite temperature and pressure. Therefore, (N, P, σ = 0, T) simulations at 300 K and 100 kPa were carried out for each of the 11 materials, and the internal coordinates were averaged over 500 ps following a 50 ps equilibration run. For the pristine UiO-66 and UiO-67, these results can be compared with accurate single-crystal X-ray diffraction (SCXRD) data,[41,42,56,104,105] while, to the best of our knowledge, no such accurate data are available for the UiO-68 due to the more difficult synthetic conditions. Furthermore, defects are not well resolved in X-ray diffraction data, due to the random nature of these defect clusters, so that on average XRD data of defect materials will more closely resemble the defect-free material compared to one specific defect structure. In Table , comparison with SCXRD data of ref (41), which provides data for both UiO-66 and -67, is made for some interesting bond lengths, angles and torsions in both defect-free materials. It should be noted that the UiO-66-type materials are very sensitive to the synthesis and activation procedures, resulting in small deviations in the internal coordinates when comparing different experimental results. More details on the synthesis and activation procedures for the experimental sample discussed here can be found in ref (41). As shown in Table , most interatomic distances and all interatomic angles are reproduced within an error of 2%. The largest deviation, 4.71%, is obtained for the bond length between a zirconium atom and the neighboring carboxylic oxygen, but is still acceptable. Moreover, this deviation can be traced back partially to the DFT input data used to generate the force field, which yields an optimized distance of 2.277 Å at 0 K (see Table S21). These MD results also confirm that the inorganic zirconium brick barely changes between UiO-66 and UiO-67, validating the isolated building block approach and the choice of cluster models.

Table 1

Comparison between Our (N, P, σ = 0, T) MD Simulations and Single-Crystal X-ray Diffraction Data of a Selected Set of Internal Coordinates for UiO-66 and -67 at 300 K and 100 kPa (see Figure for the Definition of the Atom Types for the Pristine Materials)

interatomic distance [Å]
	UiO-66			UiO-67
	this work	SCXRD[41]	rel. diff. [%]	this work	SCXRD[41]	rel. diff. [%]
Zr···Zr	3.509(2)	3.513b	–0.11	3.513(3)	3.512b	+0.03
Zr···Ooh	2.269(2)	2.259(7)	+0.44	2.271(3)	2.254(5)	+0.75
Zr···Oox	2.086(3)	2.065(3)	+1.02	2.085(3)	2.059(2)	+1.27
Zr···Oca	2.3140(4)	2.210(5)	+4.71	2.318(3)	2.218(1)	+4.51
Oca···Cca	1.2719(1)	1.259(4)	+1.02	1.2729(3)	1.269(1)	+0.31
Cca···C1	1.5134(1)	1.495(4)	+1.23	1.5126(2)	1.494(2)	+1.24
C1···C2	1.4105(1)	1.379(3)	+2.28	1.4103(2)	1.386(2)	+1.75
C2···C3a	1.3939(1)	1.389(5)	+0.35	1.3951(3)	1.386(2)	+0.66
C3···C4				1.4179(3)	1.392(2)	+1.86
C4···C4				1.5077(2)	1.489(4)	+1.26

For UiO-66, the atom types C2 and C3 coincide.

No experimental errors reported in ref (41).

In addition to the internal coordinates, also the unit cell parameters were determined from (N, P, σ = 0, T) MD simulations. As shown in Table , cubic unit cells are retrieved for UiO-66 to -68, with unit cell parameters which are slightly overestimated by about 2%. This overestimation arises from the DFT calculations that were used for the fitting of the force field parameters, and is a well-known problem for these systems.[106] Moreover, all pristine unit cells belong to the Fm3̅m space group, in accordance with experimental results,[41] but different from periodic DFT results at 0 K which predict space groups with lower symmetry because of the fixed position of the linkers at 0 K.[33,107−109] The same analysis for the defect structures is carried out in Table S23, showing that the creation of linker vacancies distorts the cubic unit cell, resulting in an appreciable decrease in symmetry.

Table 2

Unit Cell Properties for UiO-66 to -68 as Calculated via (N, P, σ = 0, T) MD simulations at 300 K and 100 kPa, Compared to SCXRD Data When Available[41]

	cell length [Å]		cell volume [Å3]		space group
material	this work	SCXRD[41]	this work	SCXRD[41]	this work	SCXRD[41]
UiO-66	21.117(1)	20.743(3)	9 416(2)	8 925(4)	Fm3̅m	Fm3̅m
UiO-67	27.328(2)	26.883(3)	20 407(5)	19 428(7)	Fm3̅m	Fm3̅m
UiO-68a	33.434(3)	32.7767(5)	37 374(10)	35 212(2)	Fm3̅m	−

For UiO-68, the SCXRD parameters for the amine-functionalized variant of ref (104) are reported due to the absence of single crystals of nonfunctionalized UiO-68.

For many applications based on gas adsorption or relying on the mechanical rigidity of the material, a good reproduction of its nanoporous structure is a prerequisite. Using Zeo++, the pore size distributions of the three pristine and the eight defect materials are determined based on (N, P, σ = 0, T) simulations at 300 K and 100 kPa. As shown in Figure S7, each structure exhibits two peaks in the pore size distribution, corresponding to the expected locations of the tetrahedral and octahedral pores, and the size of both types of pores increases when going from UiO-66 to -68.[32] Moreover, when removing linkers, both peak heights decrease. Indeed, when creating a linker defect, two octahedral and two tetrahedral pores merge, creating a pore which largest included sphere is identical to the sphere present in the octahedral pore. Hence, the net result of creating a linker vacancy is the removal of two tetrahedrally sized pores and one octahedrally sized pore, such that the tetrahedral peak decreases twice as fast as the octahedral one. This is indeed observed in Figure S7.

Mechanical Behavior of the Pristine Materials

Using the validated first-principles force fields, the pressure-induced behaviors of the UiO-66-type materials are studied by constructing pressure profiles as a function of the constrained unit cell volume for the defect-free UiO-66, -67, and -68 materials via (N, V, σ = 0, T) simulations at 300 K. In Figure , the resulting pressure profiles and their corresponding free energy profiles, obtained by thermodynamic integration,[110] are shown. The equilibrium volumes at 300 K can be read off readily from the pressure profile by determining the intersection of each of these curves with the horizontal P = 0 MPa. These equilibrium volumes, which correspond to the minima of the free energy profiles, are listed in Table . By determining the derivative of the pressure profile, the bulk moduli of the stable structures can be obtained, yielding 22.2, 13.3, and 8.1 GPa for UiO-66 to -68, respectively (see Table ). While the obtained bulk moduli are appreciably smaller than DFT values (39.5–41.0 GPa for UiO-66, 17.4–22.1 GPa for UiO-67),[28,46,56,109] they are in line with the experimental bulk modulus of 17(1.5) GPa for an 11-fold coordinated UiO-66.[35] From Table , it is clear that the bulk modulus decreases with increasing linker length. This observation was expected, since an increase of the linker length also leads to an increase in internal pore volume, weakening the material.[111]

Figure 4

Top: Internal pressure ⟨P⟩ as a function of the constrained unit cell volume V for UiO-66 to -68, resulting from (N, V, σ = 0, T) simulations at T = 300 K. Bottom: The corresponding free energy profiles F as a function of the constrained unit cell volume V, obtained via thermodynamic integration.

Table 3

Equilibrium Volumes, Bulk Moduli, and Loss-of-Crystallinity Pressures for UiO-66, -67, and -68, Based on (N, V, σ = 0, T) Simulations at T = 300 K

material	cell volume [Å3]	bulk modulus [GPa]	loss-of-crystallinity pressure [GPa]
UiO-66	9 419	22.2	1.83
UiO-67	20 441	13.3	0.45
UiO-68	37 400	8.1	0.20

Two regimes may be distinguished in the profiles of the three UiO materials, as evident from Figure . Around equilibrium, a quasi linear pressure-versus-volume dependence is retrieved, corresponding to a parabolic free energy profile. This indicates that the materials satisfy Hooke’s law at volumes close to equilibrium. For UiO-66, a deviation from this linearity is observed around 8400 Å3, causing a maximum in the P(V) profile at a pressure of 1.83 GPa. At this volume, a bend in the pressure profile is observed, after which the pressure decreases with decreasing volume. In this regime, −∂P/∂V is negative and an unstable branch is encountered. The local maximum in the pressure profile can hence be associated with the onset of mechanical instability and correlates well with the pressure at which the experimental loss of crystallinity was observed. Moreover, when pushed into this mechanically unstable regime, a reduction of short-range order (vide infra) is observed, indicating a short-range loss of crystallinity. Hence, the maximum of this curve will be referred to as the loss-of-crystallinity pressure. Once a pressure higher than this maximum is applied, a sudden drop in volume is observed. While these results on itself do not elucidate on the exact nature of the amorphous phase, and may also be explained by a phase transformation to a crystal with a lower symmetry, they coincide with the experimental loss of crystallinity described by Yot et al., who observed a similar drop in the unit cell parameter of an 11-fold coordinated UiO-66 at about 1.4 GPa.[35] For even lower volumes, one again expects to find a stable branch.

To validate these results and identify the soft mode responsible for the observed mechanical instability, the constant-pressure approach proposed in ref (73) has been applied (see Section S4.3 for the computational details). To study the effect of an elevated pressure on the stability of the defect-free UiO-66 material, 21 extended (N, P, σ = 0, T) simulations were carried out on the conventional cell with a fixed pressure between 0 and 2000 MPa, using a step of 100 MPa. For a cubic material subject to a hydrostatic pressure P to be stable, the following three Born stability criteria should be met:[112]where C11, C12, and C44 are the three independent elastic constants for a cubic material. As shown in Figure , the second Born criterion is violated first, near a pressure of about 1.8 GPa. This pressure indicates the onset of mechanical instability, which is caused by compression, since C11 – C12 – 2P < 0, and corresponds very well with the loss-of-crystallinity pressure of 1.83 GPa determined using our fixed-volume approach.

Figure 5

The three Born stability criteria (eq ) for the pristine UiO-66 as a function of the applied pressure P during (N, P, σ = 0, T) simulations at T = 300 K.

The behavior for the two other pristine materials, UiO-67 and -68, is similar but less evident. The drop of the P(V) profile takes place close to the equilibrium volume when systematically decreasing the volume, and the profile behaves almost constant when decreasing the volume even further. The loss-of-crystallinity pressures for the three pristine materials are reported in Table .

In ref (35), it was shown that the sudden decrease of the unit cell parameter in UiO-66 is accompanied by a broadening of the Bragg peaks, demonstrating a pronounced loss of crystallinity. To verify whether the same behavior can be obtained using our newly developed force fields, the average structures at fixed volumes ranging between 8000 and 9700 Å3 were determined from a number of (N, V, σ = 0, T) MD simulations, corresponding to the volume range for which the P(V) curve was obtained. For each average atomic structure, the symmetry of the unit cell was determined with PLATON,[103] and the number of symmetry operators was calculated. At volumes near equilibrium, the pristine materials can be assigned the Fm3̅m space group, containing 192 symmetry operators. As shown in the top pane of Figure , the number of symmetry operators of UiO-66 sharply drops when the volume of the material is forced to be smaller than ∼8400 Å3. Comparison of this profile with Figure reveals that this sudden decrease in symmetry occurs near the volume for which the maximum in the pressure profile is obtained. Hence, our (N, V, σ = 0, T) simulations indicate that the decrease in unit cell volume is associated with a sudden loss of symmetry, confirming the experimental results.[35] Also for UiO-67 and -68, a distinct loss of symmetry is observed at those volumes for which their pressure curves start to deviate from the linear behavior. These bends occur at appreciably lower pressures of about 0.45 and 0.20 GPa for UiO-67 and UiO-68, respectively. For UiO-67, Hobday et al. experimentally observed a loss of crystallinity at about 0.3 GPa, in line with our results.[56]

Figure 6

Log10 of the number of symmetry operators as a function of the constrained unit cell volume for UiO-66 to -68, resulting from (N, V, σ = 0, T) simulations at T = 300 K. The space groups corresponding to both the low- and high-volume limit are also indicated.

To get more insight in this loss of symmetry and detect at what distance from the inorganic brick the crystalline order is broken, a noninteracting dummy atom is introduced at the center of a zirconium brick during an (N, V, σ = 0, T) simulation at 300 K. Subsequently, the radial distribution function (RDF) of every atom with respect to this dummy atom is determined. This procedure is carried out at two volumes: (i) at 9400 Å3, where an intact crystalline order is expected, and (ii) at 8000 Å3, where the loss of crystallinity has taken place. The two RDFs are shown in Figure . For the RDF at a volume of 9400 Å3, sharp peaks can be resolved over a large range of distances, and each peak can be assigned to a unique atom type. Even at distances as high as 17 Å from the dummy atom, zirconium atoms of the other inorganic bricks can be distinguished. A completely different picture emerges for the RDF at a volume of 8000 Å3. Here, the inorganic brick and carboxylate group closest to the dummy atom are still nicely resolved. However, from a distance of about 6 Å onward, corresponding to the nearest carbon atom of the phenyl ring, the peaks start to spread out and overlap with each other. Hence, the observed loss of symmetry corresponds to a short-range loss of crystallinity, which already manifests itself at the connection between the inorganic brick and the first organic ligand. An analogous behavior was observed in ZrCDC, a material isoreticular to UiO-66, based on a cyclohexanedicarboxylate ligand.[38] Note that this also implies that our simulation cell, with cell lengths larger than 20 Å, will be able to capture this loss of crystallinity, since the loss of crystallinity already occurs at much smaller distances. This is also explicitly validated using a 2 × 2 × 2 supercell (see Section S5).

Figure 7

Radial distribution functions of UiO-66 with respect to the center of an inorganic brick, resulting from (N, V, σ = 0, T) simulations at T = 300 K and V = 8000 Å3 (top) or V = 9400 Å3 (bottom). Zirconium atoms are cyan, oxygen atoms red, hydrogen atoms black, and carbon atoms gray.

Effect of Linker Vacancies on the Mechanical Behavior

The aforementioned results correspond to the pristine materials. However, for UiO-66, experimentally synthesized materials are known to contain linker vacancies, which may impact the high stability of this material. To investigate the effect of structurally ordered defects in the UiO-66 material, pressure-versus-volume curves for the eight classes of defect-containing materials shown in Figure were constructed. In Figure , the polynomial fits to these pressure profiles are shown, while the extracted parameters are listed in Table .

Figure 8

Internal pressure ⟨P⟩ as a function of the constrained unit cell volume V for the defect-free UiO-66 and the different defect-containing materials, resulting from (N, V, σ = 0, T) simulations at T = 300 K, with indication of the initial decrease in loss-of-crystallinity pressure when introducing the first linker vacancy. Insets: Schematic depiction of the type 3 and type 5 structures.

Table 4

Equilibrium Volumes, Bulk Moduli, and Loss-of-Crystallinity Pressures for the Defect-Free UiO-66 and the Different Classes of Defects, Based on (N, V, σ = 0, T) Simulations at T = 300 K

material	cell volume [Å3]	bulk modulus [GPa]	loss-of-crystallinity pressure [GPa]
defect-free	9 419	22.2	1.83
type 0	9 389	19.9	1.55
type 1	9 361	17.4	1.29
type 2	9 360	18.2	1.37
type 3	9 385	18.7	1.51
type 4	9 330	18.2	1.39
type 5	9 383	15.5	1.17
type 6	9 337	18.9	1.38
type 7	9 390	17.2	1.35

As expected, introducing linker vacancies lowers the equilibrium volume only slightly, but has a profound effect on both the bulk modulus and the loss-of-crystallinity pressure. After removing a first linker, the bulk modulus drops from 22.2 to 19.9 GPa, whereas the removal of a second linker leads to a bulk modulus in the range of 15.5–18.9 GPa, confirming the trends observed in recent DFT results.[54] Likewise, a first linker defect decreases the loss-of-crystallinity pressure from 1.83 to 1.55 GPa, as indicated by the black arrow in Figure , while removal of a second linker yields a loss-of-crystallinity pressure between 1.17 and 1.51 GPa. Note that the ranges of both bulk moduli (15.5–18.9 GPa) and loss-of-crystallinity pressures (1.17–1.51 GPa) agree very well with the experimental values of 17(1.5) GPa and 1.4 GPa, respectively.[35]

While most of the structures containing two defects give rise to a similar pressure behavior, with a bulk modulus of 17.2–18.9 GPa and a loss-of-crystallinity pressure between 1.29 and 1.39 GPa, two types of defects have a more pronounced effect on the stability of the material. For the type 3 defect, obtained by removing two linkers in such a way that a 1D channel is created in the material (see right inset of Figure ), the mechanical stability remains exceptionally intact, with a bulk modulus and loss-of-crystallinity pressure of 18.7 GPa and 1.51 GPa, respectively, which are both close to the properties of the type 0 defect, with only one missing linker. This can be explained by considering that, although a 1D channel is formed, the walls constituting this channel remain reinforced by surrounding linkers, and no weak directions are created. In contrast, for the type 5 defect, a relatively strong deterioration of the mechanical properties is observed, yielding a lower bulk modulus of 15.5 GPa and a lower loss-of-crystallinity pressure of 1.17 GPa. This effect is due to the equal orientation of the two removed linkers on two different planes in a type 5 defect (see left inset of Figure ). It is hence possible to compress half of the unit cell perpendicular to this direction by a shearing motion, compromising the stability of the material. No other type of structure with two defects shares this property. Hence, this analysis shows that not only the number of linker vacancies plays an important role for the stability of the material, but also their position throughout the unit cell may be of paramount importance.

Conclusions

In this work, accurate force fields were developed for the highly intriguing UiO-66 family of materials, using the QuickFF protocol for the covalent force field terms. For an accurate description of UiO-67 and UiO-68, specific force field terms were introduced to describe the peculiar shape of the dihedral patterns of the ligands. The obtained force fields were validated based on both DFT and experimental single-crystal X-ray diffraction data. It was shown that the internal coordinates, unit cell properties and the porous structure of each MOF, calculated via (N, P, σ = 0, T) simulations, were in good agreement with already published data, which assures the applicability of the force fields to study the pressure behavior of the UiO-66-type materials.

Subsequently, pressure and free energy profiles as a function of the unit cell volume were obtained via (N, V, σ = 0, T) MD simulations on the defect-free UiO-66 to -68 at 300 K. These simulations point toward the experimentally observed loss of crystallinity and accompanying drop in unit cell volume for UiO-66. For the isoreticular UiO-67 and UiO-68, similar effects are observed albeit at lower pressures. Applying an external pressure clearly induces short-range disorder, an aspect which was further investigated by considering radial distribution functions at volumes corresponding to pressures above and below the loss-of-crystallinity pressure. When the applied pressure is higher than the loss-of-crystallinity pressure, the radial distribution function reveals that at about 6 Å the crystalline order of the material is already lost, which is a far smaller distance than the periodicity used in the simulations, indicating a short-range loss of crystallinity.

Since it is generally accepted that inclusion of linker defects may substantially alter the material properties, a series of structures with one or two defects were constructed to elucidate the influence of linker vacancies on the mechanical properties. These defect structures were classified into eight unique types, of which one corresponds to an 11.5-fold coordination, while seven correspond to an 11-fold coordinated inorganic brick with distinct combinations of linker defects. While these defects have only a minor effect on the equilibrium volume, they introduce a substantial and gradual decrease in both the loss-of-crystallinity pressure and the bulk modulus. For the structures with two linker vacancies, the obtained mechanical properties are in excellent agreement with the experimental findings of Yot et al.[35] The here applied procedure gives unprecedented insight in the influence of the relative orientation and position of the two linker defects on the stability of the material. The most profound effect on the stability is obtained when the two linker vacancies share the same orientation but lie in neighboring lattice planes. In contrast, when the two defects are created by removing two perpendicular ligands which have no inorganic bricks in common, the structural integrity of the material is preserved to a large extent, sharing similarities with the structure containing only one linker defect.

In conclusion, our force field based approach succeeds in yielding a molecular level insight in the structural short-range order of framework materials upon application of an external pressure. Extending the organic linker or introducing linker vacancies enables to tune the mechanical behavior and stability of the materials. The true challenge consists in tuning the materials toward specific applications by the intentional creation of linker vacancies. The first-principles force field based approach introduced here is a strong tool to aid the engineering of materials toward specific applications, for instance for catalysis or gas adsorption, while ensuring to a large extent the retention of their structural integrity.[25,113−116]