The impact of lattice vibrations on the macroscopic breathing behavior of MIL-53(Al)

IR spectra

Experimentally, only the IR spectra of the (meta)stable CP and LP phases have been measured as an in situ measurement of the IR spectrum during the phase transformation is not yet possible. These experimental reports describe several IR active vibrations that distinguish the CP phase from the LP phase [24], [40], [41]. It is observed that the vibrational frequencies of most of these vibrations increase after transition from the CP to the LP phase and thus cost more energy. In our theoretical approach, it is possible to predict the volume-dependency of these IR fingerprints. The calculated IR spectrum of MIL-53(Al) as a function of the volume is represented in Figure 2a. The fingerprint vibrations are labeled in correspondence with our earlier work [24] (Labels from A to I, further explained in the caption of Figure 2). To facilitate the recognition of specific trends, the volume-frequency and volume-IR-intensity relations of these specific modes are represented in Figure 2b and c, respectively.

Fig. 2:

(a) Theoretical IR spectra of MIL-53 (Al) at different unit cell volumes within the wavenumber range 0–1800 cm⁻¹. The inset zooms in upon the range 280–780 cm⁻¹. The region of the aluminum-oxide stretch vibrations is highlighted with a gray background. Fingerprint vibrations, identified in earlier reports [24], [40], [41], are indicated with a capital letter and follow the same labeling as in Ref. [24]. (A1) Symmetric stretching of the aluminum-oxide backbone. (A2) Antisymmetric stretching of the aluminum-oxide backbone. (B) Out-of-plane wagging of hydrogen atoms on the aromatic ring. (C) Bending of the carboxyl group and within the aromatic ring. (D) Rocking of hydroxyl group. (E) Bending within the aromatic ring. (F) Stretching between the aromatic ring and the carboxyl group. (G) Symmetric stretching of the carboxyl group. (H) Stretching between the aromatic ring and the carboxyl group. (I) Antisymmetric stretching of the carboxyl group. (b) Volume-frequency relations of the IR fingerprint vibrations. The wavenumbers are shifted by their average value over the computed volume range. Dotted lines indicate the equilibrium volumes of the CP and LP phases obtained from electronic structure calculations. The volume ranges of the three different volume-frequency regimes are indicated. (c) Volume-IR-intensity relations of the IR fingerprint vibrations. Dotted lines indicate the equilibrium volumes of the CP and LP phases obtained from electronic structure calculations.

Starting with the IR intensities as a function of the volume (Figure 2c), no clear trend can be observed. Moreover, the intensity difference between the vibrational modes in the CP and the LP phases is small, except for mode A1 (symmetric stretching of aluminum-oxide backbone), which doubles in intensity when going from the CP toward the LP phase.

A more interesting behavior is observed for the volume-frequency relations. The vibrational frequencies of the majority of the fingerprint vibrations follow a certain trend as a function of the volume. Three regimes can be distinguished (see Figure 2). In the first regime, at low unit cell volumes, the frequency strongly decreases as a function of the volume. A second regime starts around the volume where the CP phase is stable (843 ų from electronic structure calculations). From this point onwards, the vibrational frequencies of most normal modes steadily increase until the volume region of the stable LP phase is encountered (1426 ų from electronic structure calculations). Finally, the frequency drops again in the third regime at high unit cell volumes. The strong decrease in frequency as a function of the volume in regimes 1 and 3 can be explained by the weakening of covalent interactions. This can be illustrated by looking at the example of the aluminum-oxide bond (Figure 3a–c). At very low and high volumes, the bond distance strongly increases as a function of the volume indicating a decrease in the strength of the covalent bond. At the same time, the frequencies of all stretching vibrations of the aluminum-oxide bond decrease (the corresponding frequency region is highlighted by a gray background in Figure 2a). Within regime 2 the strength of the covalent interactions remains more or less unaltered as the bond distances remain almost constant. Nevertheless, an increase in the vibrational frequency is observed as a function of the volume. This is a result of dispersion interactions, which are important in MIL-53(Al) [15] (see Supporting Information Section S2).

Fig. 3:

(a) Visualization of the aluminum-oxide bonds. (b) Volume-frequency relations of all vibrational modes describing stretches of aluminum-oxide bonds. (c) Length of the aluminum-oxide bonds as a function of the volume. (d) Visualization of the hydroxyl bond. (e) Volume-frequency relations of all hydroxyl stretch vibrations. (f) Length of the hydroxyl bond. Dotted lines indicate the equilibrium volumes of the CP and LP phases obtained from electronic structure calculations.

Some modes differ from the general volume-frequency trend, for example the wagging mode of the CH-bonds (mode B in Figure 2), the antisymmetric stretching mode of the carboxyl group (mode I in Figure 2), and the bend and stretch modes of the hydroxyl group around 970 cm⁻¹ (D in Figure 2) and 3750 cm⁻¹ (Figure 3e), respectively. Limiting the discussion to the last two examples, it is observed that the hydroxyl bend vibration shows a decreasing frequency for increasing volumes at small unit cell volumes and afterwards the frequency increases, but instead of going down again it keeps on increasing. The hydroxyl stretch vibration shows the opposite behavior: at small volumes, the frequency increases upon increasing volume, but it starts decreasing at larger volumes. Those modes are extremely localized, which makes them only dependent on the geometry of the hydroxyl group. Looking at the volume-frequency relations of the hydroxyl stretch vibrations and the bond length of the hydroxyl group as a function of the volume (Figure 3e and f, respectively), there exists indeed a clear correlation between the bond length and the vibrational frequency of the hydroxyl group.

For wavenumbers between 300 and 500 cm⁻¹, the IR spectrum as a function of the volume shows more abrupt changes. This is caused by a change in IR activity of specific modes, rather than deviations from the general volume-frequency trend (see Supporting Information Section S3). In our previous work [24], it was argued that some of the changes in IR activity originate from rocking of the hydroxyl group, which is present at large volumes, but suppressed at small volumes. This claim is confirmed in our simulations.

Terahertz vibrations

The terahertz region (<300 cm⁻¹) is especially interesting for soft porous crystals as it contains the modes which influence the dynamics of the framework most profoundly. In this frequency region, the IR spectrum of MIL-53(Al) contains only a few IR active modes, so we have to restrict our discussion to the volume-frequency relations to obtain a complete view of the vibrational spectrum (Figure 4). It is clear that the behavior of low-frequency modes differs from the behavior of high-frequency modes as it concerns mainly collective vibrations. Here, the modes showing the most interesting volume-frequency curves are highlighted. First of all, the vibrations describing linker rotations (indicated in blue in Figure 4) and trampoline-like motions (red curves in Figure 4) are discussed. Similar terahertz vibrations were recently reported around the same frequencies in MIL-140A, a Zr-based MOF with BDC linkers [42]. We find in our study that such modes display a large increase in frequency for decreasing unit cell volume, as these vibrations are heavily constrained at small unit cell volumes. Among the vibrations exhibiting linker rotations, two types can be distinguished, which do also appear in the MIL-140A framework [42]. Namely, those where opposite linkers rotate in the same direction and those where opposite linkers rotate in the opposite direction. The former are heavily constrained in dense structures yielding high vibrational frequencies, while the latter are less restricted in small unit cells. This causes a splitting of their volume-frequency relations around 1200 ų. It has to be noted that, according to our simulations, the vibrational frequencies of all linker rotation modes show a sudden decrease when increasing the volume above 1520 ų. This decrease in frequency is induced by a decrease in the dihedral angles of the Al-O-C-O units, illustrated in Figure 5. For large volumes this dihedral angle approaches zero, which lowers the energy barrier for linker rotations and thus its corresponding vibrational frequency. This sudden decrease in frequency results in a large increase of the vibrational entropy (see Section “Thermodynamic properties”) which may trigger the CP-to-LP phase transformation.

Fig. 4:

Volume-frequency relations of the terahertz vibrations in MIL-53(Al). The collective vibrations with the most significant frequency differences over the computed volume range are highlighted and visually represented. Blue curves: vibrations inducing linker rotations. Green curves: soft modes combining translations and rotations of the aluminum-oxide backbone and the linkers. Yellow curves: vibrations inducing rotations of the aluminum-oxide backbone. Red curves: trampoline-like vibrations.

Fig. 5:

(a) Visualization of the Al-O-C-O dihedral. (b) Volume-frequency relations of the terahertz vibrations exhibiting rotations of the linkers. (c) Different Al-O-C-O dihedrals.

Furthermore, there are three modes showing rotations of the metal-oxide backbone (yellow curves in Figure 4). This rotation is also constrained at small volumes, hence the increase in frequency upon lowering the volume, but the frequency also increases at large volumes. This probably results from linkers which become overstretched when the backbone rotates at large volumes. Finally, two modes have a completely different volume-frequency relation, showing a continuous increase in frequency upon increasing volume (green curves in Figure 4). Both modes consist of a combination of linker and metal-oxide chain movements. The mode that has the lowest frequency over the entire volume range can be described by a rotation of the metal-oxide backbones, which induces a trampoline-like motion of the linkers (see green box in Figure 4 for a visualization). The other mode shows a translation of the metal-oxide backbones, inducing a rotation of the linkers (see Supporting Information Section S4 for a visualization). The volume-frequency relations of these two modes are strongly coupled with the volume-dependency of the diagonal elements of the Lagrangian strain tensor η [43] (see Figure 6). The latter can be obtained by taking matrix products of the cell matrices at a reference unit cell volume, h₀, and a deformed unit cell volume, h(V):

Fig. 6:

The strain within MIL-53(Al) along the three Cartesian directions as a function of the volume for structures with the CP phase (843 ų) as reference volume (a) and the LP phase (1426 ų) as reference volume (b). The volume-dependency of the strain along the aluminum-oxide chain, η_xx, is small and proportional to the volume-frequency relation of the soft mode exhibiting rotations of the metal-oxide backbone. The volume-dependency of the strains along the pore diagonals, η_yy and η_zz, is large and proportional to the volume-frequency relation of the soft mode exhibiting translations of the metal-oxide backbone.

The vibrational frequency as a function of the volume of the mode exhibiting rotations of the aluminum-oxide chain is strongly linearly correlated with the strain along this chain, η_xx, as a function of the volume. Moreover, a linear relation is observed between the volume-dependency of the vibrational frequency of the mode characterized by translations of the aluminum-oxide backbone and the strain along the pore diagonals, η_yy and η_zz, as a function of the volume. This correspondence indicates that the particular volume-frequency relations of these two terahertz vibrations originate from the strain developed within the framework. As such, this procedure allows us to reveal the relevant vibrational modes that are most affected by the breathing-induced strain in the material, a macroscopic property, through the strain tensor. As both modes become very soft at small volumes, they are very important for the dynamics in the CP phase of MIL-53(Al). The vibration inducing rotations of the aluminum-oxide backbone, for example, can trigger the transition from a closed pore phase to a very closed pore phase as reported in MIL-53(Sc) [44]. In the remainder, we will refer to these two modes as soft modes.

Helmholtz free energy profile

Starting from the analytic expressions for the volume-frequency relations, obtained after fitting a ninth-order Taylor expansion to the data points, the Helmholtz free energy profile as a function of the volume can be constructed at different temperatures (Figure 7a). Those results are compared with the electronic energy profile without vibrational contributions such as the zero-point energy (ZPE) (black curve).

Fig. 7:

Helmholtz free energy as a function of the volume at different temperatures (solid lines) and the electronic energy as a function of the volume (black dashed line) using the QHA (a) (this work) and dynamic methods (b) (Figure reproduced from Ref. [45] with permission of the American Chemical Society).

Within the considered temperature window, we recognize two minima in the Helmholtz free energy profile corresponding to the presence of the (meta)stable CP and LP phases. At 1 K, the equilibrium volumes of the CP and the LP phases are located around 860 and 1461 ų, respectively. Using another degree for the fitted polynomials can shift these values by 5 ų, introducing an uncertainty on the predicted values. Note that, even at 1 K, the Helmholtz free energy profile differs substantially from the electronic energy profile. The main effect of the ZPE is a shift of the equilibrium volumes of both the CP and LP phases by about 20 ų. The value for the CP phase volume corresponds well with the experimental value measured by Liu et al. [18] at 77K, while the LP phase volume is slightly overestimated in our simulations. The free energy difference between the (meta)stable phases is nearly unaffected by the ZPE and has a value of 26 kJ·mol⁻¹ at 1 K. This free energy difference falls in the range of the results obtained by other computational studies [15], [25], [45].

By looking at the temperature response of the Helmholtz free energy profile, a slight increase in equilibrium volumes is observed for both the CP and the LP phases for temperatures above 50 K. More pronounced is the decrease of the free energy difference for increasing temperatures. It reaches values of 19 and 14 kJ·mol⁻¹ at 300 and 500 K, respectively. This behavior is expected as the LP phase is entropically favored over the CP phase [15], [18], [25]. However, even for temperatures as high as 500 K, the LP phase has not become the most stable phase, while the experimental CP-to-LP transition temperature is found around 350 K [18]. This deviation with the experimental observations is caused by an exaggerated stability of the CP phase due to the applied computational methodology [15]. The PBE-D3(BJ) level of theory, used in this work, overestimates the dispersion interactions, which results in an unphysical stabilization of the CP phase. This point was extensively studied in a recent contribution by some of the presenting authors [15]. When the random-phase approximation is applied, a higher level of theory for a more accurate description of correlation effects, the electronic energy difference between the CP and the LP phase reduces to about 7 kJ·mol⁻¹ [15]. Given the excessive computational cost and the current computational resources, a vibrational analysis using the latter approach is not feasible. Nevertheless, starting from the energy difference at the RPA level of theory and considering the vibrational frequencies as predicted in this work using PBE-D3(BJ), the CP and LP phases become equally stable around 300 K (see Supporting Information Section S5), which is more in line with experimental observations [18]. Furthermore, comparison of the lattice vibrations at the PBE-D3(BJ) level of theory with a hybrid functional did not yield significant changes (see Supporting Information Section S5). Finally, in our earlier study a good agreement was found between the experimental IR and Raman spectra and the computational spectra at the PBE-D3(BJ) level of theory [24]. This hints at a reliable prediction of the lattice dynamics at the current level of theory. To conclude this paragraph, it is worth mentioning that the Helmholtz free energy profile in Figure 7a is constructed via the QHA, which still neglects certain anharmonic effects. By including these effects via dynamical methods, which are computationally more expensive, an increased entropic stabilization of the LP phase is observed [45] (see Figure 7b). Although, it has to be noted that the results obtained from these ab initio molecular dynamics simulations are prone to large sampling uncertainties.

Energetic and entropic contributions

The specific temperature dependence of the Helmholtz free energy profile can be more thoroughly understood when considering the contributions of the electronic energy, the vibrational energy and the entropy to the free energy. The free energy difference between the LP and the CP phases of these contributions and the total free energy difference as a function of temperature are presented in Figure 8. At low temperatures, the total free energy difference between the LP and the CP phases is mainly the result of the difference in electronic energy. Furthermore, this plot illustrates that the difference in the vibrational energy between the LP and the CP phases is limited over the complete temperature range, whereas the difference in entropic contribution decreases significantly, canceling out the increase in the electronic energy difference. As temperature alters the equilibrium LP and CP volumes, also the electronic energy difference between these two phases will change even though the electronic energy in itself is temperature independent at a fixed volume. The steady decrease in the difference of the entropic contribution originates from the activation of the vibrations inducing linker rotations, which have very low frequencies at large unit cell volumes. As a consequence, the entropic stabilization of the LP phase is mainly the result of this type of vibrations, as already suggested by Liu et al. [18].

Fig. 8:

Total free energy difference between the LP and CP phases as a function of temperature, ΔF_tot (black curve), and its three contributions: entropy −TΔS (blue), vibrational energy ΔE_vib (red), and electronic energy ΔE_el (green).

Contributions from vibrational modes

Up to this point, we could identify the vibrations exhibiting linker rotations as important contributors to the Helmholtz free energy profile, but they are not the only ones shaping the free energy landscape. This becomes clear when analyzing the vibrational free energy profile of each mode separately (Figure 9). At 1 K, the vibrational free energy is restricted to the ZPE. Several terahertz modes undergo a large frequency change as a function of the volume and, therefore, they yield ZPE differences over the calculated volume range as large as 1.2 kJ·mol⁻¹. Nonetheless, their impact on the total free energy profile remains small as also multiple localized vibrations give rise to a substantial ZPE difference. At 300 K, the picture looks completely different. Now, the vibrational free energy differences are dominated by the terahertz vibrations. First in line are the vibrations exposing linker rotations (blue curves) which generate free energy differences of at least 4 kJ·mol⁻¹, favoring open pore structures. However, these modes are counteracted by the soft modes (green curves), showing combined collective behavior of the linkers and the metal-oxide backbones. The magnitude of the free energy difference produced by these two vibrations are of comparable size to the free energy difference of the vibrations exhibiting linker rotations. As such, it can be expected that these modes are of crucial importance in the LP-to-CP phase transition at low temperatures and that they partially hinder the transition from the CP to the LP phase at elevated temperatures. Other modes worth mentioning are those inducing trampoline-like motions of the linkers (red curves) and metal-oxide backbone rotations (yellow curves). These modes attribute considerably to the stabilization of the LP phase at high temperatures.

Fig. 9:

Contributions of the different vibrational modes to the Helmholtz free energy at 1 K (a) and 300 K (b). The terahertz vibrations discussed in Section “Terahertz vibrations” are highlighted with the same colors as in Figure 4.

Derived properties

Besides the Helmholtz free energy profile, knowledge of the vibrational modes allows to predict other thermodynamic properties, such as the specific heat capacity, the bulk modulus, and the volumetric thermal expansion coefficient. The temperature dependence of these properties will be investigated in this subsection for both the CP and the LP phases. For a discussion of the pressure-versus-volume (P(V)) equation of state, we refer to the Supporting Information Section S6. As the bulk modulus and the volumetric thermal expansion coefficient are derived from the Helmholtz free energy profile, uncertainties on the latter due to fitting may accumulate. Therefore, the predicted values, obtained by fitting the Rose-Vinet equation of state [46] to the Helmholtz free energy profiles around the equilibrium volumes, will be presented together with an uncertainty region determined by the fitting of different orders of Taylor expansions to the volume-frequency relations.

Specific heat capacity

The specific heat capacity is a measure for the ability of a material to store heat. A large specific heat capacity will allow the material to absorb more heat and, as a consequence, increase its resistance against undesired thermal effects. The specific heat capacity at constant pressure for both the CP and LP phases are shown in Figure 10a in the temperature range 0–500 K. The curves of the (meta)stable phases coincide almost over the complete temperature range, except for very low temperatures. In that region, the CP phase has a somewhat lower value as the heat capacity is dominated by the lowest vibrational frequency. The specific heat capacity increases monotonically with temperature from 0 Jg⁻¹K⁻¹ at 0 K towards 1.4 Jg⁻¹K⁻¹ at 500 K. These values remain unaffected by applying different fitting polynomials. Our predictions overestimate the experimental value obtained by Kloutse et al. only by 0.1 Jg⁻¹K⁻¹ [47].

Fig. 10:

Thermodynamic properties of the CP (blue) and LP (red) phases as a function of temperature. The blue and red solid lines represent the results from a Rose-Vinet fit [46] (except for the specific heat capacity). The shaded area represents the uncertainty due to the used fitting polynomial. (a) Specific heat capacity at constant pressure. (b) Bulk modulus. (c) Equilibrium volume. (d) Volumetric thermal expansion coefficient. All calculated properties are compared with experimental reference data from Refs. [13], [18], [47], [48].

Influence of lattice vibrations on thermodynamic properties

In previous paragraphs, we have argued that the thermodynamic properties are mainly affected by terahertz vibrations. As we have derived analytical expressions for the vibrational frequencies as a function of the volume, it becomes possible to determine the influence of each lattice vibration on the thermodynamic properties (see Supporting Information Section S8). In Figure 11, the vibrational contributions to the specific heat capacity at constant volume, the bulk modulus, and the volumetric expansion coefficient of the CP and LP phases are presented. Starting with the heat capacity at constant volume, we observe only minor differences between both phases. As the LP phase consists of more vibrations with very-low frequencies, it has a slightly larger total heat capacity at low temperatures, but once the contribution from these modes to the heat capacity is saturated, the total heat capacities of both phases have equal magnitudes.

Fig. 11:

Contribution of the lattice vibrations to the thermodynamic properties of the CP and LP phases as a function of temperature. The terahertz vibrations discussed in Section “Terahertz vibrations” are highlighted in the same color as in Figure 4. Heat capacity at constant volume of the CP phase (a) and the LP phase (b). Bulk modulus of the CP phase (c) and the LP phase (d). Volumetric thermal expansion coefficient of the CP phase (e) and the LP phase (f).

The lattice vibrations have a more pronounced effect on the bulk modulus and the volumetric thermal expansion. The bulk modulus of the CP phase depends mainly on the contributions from the soft modes (green curves) and a few other low-frequency vibrations, although their magnitudes remain small below 500 K. For the LP phase structure, we observe a large negative contribution to the bulk modulus from the vibrations exhibiting linker rotations (blue curves) on the one hand and a large positive contribution from the vibrations describing rotations of the aluminum-oxide backbone (yellow) and one of the soft modes (green). The magnitudes of these contributions increase with temperature reaching values of −0.7 GPa and 0.4 GPa, respectively. The large negative bulk modulus of the vibrations inducing linker rotations has a considerable impact on the total bulk modulus. Nevertheless, the latter increases for increasing temperatures. This is a result of the curvature of the electronic contribution to the Helmholtz free energy profile, which increases for increasing equilibrium temperatures, as the equilibrium volume alters.

Finally, we discuss the mode contributions to the volumetric thermal expansion. For the CP phase, we notice that the two soft modes (green) give a large negative contributions to the volumetric thermal expansion coefficient giving rise to NTE even at very low temperatures. At higher temperatures, there are more vibrational modes with a positive contribution to α_V, for example the vibrations inducing linker rotations (blue) and trampoline motions (red), yielding PTE. Also for the LP phase, the soft modes tend to lower the value of α_V, but they are counteracted by both the vibrations showing rotations of the linker (blue) and the metal-oxide backbone (yellow). Especially the modes with linker rotations have a large impact on the PTE behavior over the complete temperature range. We assume that an overestimation of the contribution of these low-frequency vibrations are the main cause of the too large α_V of the LP phase and, hence, the too large increase in equilibrium volume leading to an overestimation of the bulk modulus.

It is clear that terahertz vibrations have a considerable effect on the prediction of thermodynamic properties. Therefore, it is essential to simulate the vibrational frequencies of terahertz modes accurately, as small changes in the frequencies or their derivatives with respect to the volume can substantially affect the predicted thermodynamic properties. Although, even with high computational standards, a correct prediction of low-frequency vibrations remains hard when the potential energy surface is flat and/or anharmonicities are present [51], [52].