Assessing negative thermal expansion in mesoporous metal-organic frameworks by molecular simulation †

Most conventional materials display expansion in response to heating, so there is considerable interest in identifying materials that display the opposite behavior, negative thermal expansion (NTE). The current study investigated the temperature-induced contraction of seven mesoporous metal-organic frameworks (MOFs) of varying topology and composition, which exhibit outstanding porosity, using molecular simulation. We found exceptional NTE for the most porous MOFs, as well as a correlation between the coefﬁcient of NTE and porosity. The large molecular subunits of the MOFs were further studied to ﬁnd they intrinsically display NTE, corresponding to terahertz vibrational modes. As a result, NTE has a considerable effect on the mechanical properties of these MOFs and is an important consideration for understanding the mechanical stability of new extremely porous materials.


Introduction
Conventional dense materials are expected to expand in response to heating, however, a number of materials display the opposite behavior, shrinking with increasing temperature. 1 This phenomenon, known as negative thermal expansion (NTE), is considered counterintuitive and is observed in only few inorganic solids and certain classes of framework materials. 2 The technological importance of NTE relies on the application of this thermal volume contraction to negate the positive thermal expansion of conventional materials. 3 For example, this control of thermal expansion is crucial in materials used for dental fillings. 4 A number of physical features can cause NTE. These include phase transitions, transverse vibrational modes and rigid unit modes. 5 For example, ZrW 2 O 8 shows a remarkable contraction between 0.3-1050 K owing to coupled rotations of rigid unit vibrational modes. 1 NTE appears to occur more frequently in framework materials. In particular, several families of zeolites exhibit NTE. 6,7 Metal-organic frameworks (MOFs) are porous materials, subject to increasing attention, that can produce a wide variety of porous framework structures through the combination of organic ligand units and metal nodes. 8 Many MOFs exhibit NTE, including the prototypical MOF-5 and UiO-66. [9][10][11] MOFs have advantages over traditional dense solids due to their high porosity, with  Table 1. Classical molecular simulation was employed to predict structural changes from 40-400 K and identify the origins for NTE in these materials. We also considered the intrinsic temperature-induced structural changes of the large molecular subunits that make up the frameworks and thus contribute to their NTE. Finally, the changes in mechanical properties originating from these pronounced structural deformations were computed to provide essential insight into the stability of these materials.  field, which also represents the basis for the functional form of MOF-FF.
Coulomb interactions were truncated at 12 Å by employing the damped shifted force method from Fennell and Gezelter. 26 The MOF-FF parameterization methodology was carefully validated by reproducing the NTE behavior of MOF-5 and HKUST-1, featuring the same inorganic structural building units as the investigated mesoporous MOFs. 27,28 The framework response to temperature was simulated by N, P, T molecular dynamics employing a Nosé-Hoover chain thermostat 29 and Martyna-Tuckerman-Tobias-Klein (MTTK) barostat. 30,31 These temperature and pressure controls have been previously assessed for use in describing flexible MOF systems. 32 Systems were heated from 40 to 400 K in 20 K intervals. The timestep for the simulations was 1 fs and the thermostat and barostat dampening parameters were 100 fs and 2000 fs, respectively. Each temperature included an equilibration period of 0.5 ns followed by a production period of 1.0 ns, where the average volume and dynamics were investigated. The complete mechanical properties were calculated from the trajectories by analysis of the fluctuations of unit cell vectors. 33 The Hessians H of the optimized structures as 2 × 2 × 2 supercells were calculated with a double-sided finite difference scheme based on the analytical forces using a distortion of 0.001 Å. Based on these Hessians, phonon calculations were conducted using the phonopy software package. 34 The largest ligand was extracted from each MOF and the same force field parameters were assigned to the molecule as for the entire MOF structure. Molecular dynamics were conducted in the N,V, T ensemble, and the dynamics of the oxygen atoms of the ligand were constrained in one plane. Trajectory periods and the thermostat were identical to that of the framework simulations. The average radius of gyration of the oxygen atoms was calculated to provide a representative value of the planar size of the ligands.
Representative input files for molecular simulations are available online in our data repository at https://github.com/ jackevansadl/supp-data.

Negative thermal expansion
Dynamics of the investigated MOFs were simulated between 40-400 K, with 20 K increases in temperature, and average cell parameters were computed. As displayed in the Electronic Supplementary Information (ESI), each of the MOFs exhibit considerable contraction in response to heating, and no significant changes to the unit cell shapes were observed.
This contraction produces pronounced NTE, Figure 2a, with frameworks displaying a change in volume between 1-3 %. In general, NTE materials show total volume contractions <1.5 %, making the magnitude of contraction observed by these materials exceptional. 35 The volumetric thermal expansion coefficient (α V ), defined in Equation 1, provides a measure of the expansion rate.
The volumetric thermal expansion coefficient was calculated at different temperatures within the range investigated, and we found a relatively constant rate of contraction, Figure 2b. However, there was a noticeable change in contraction rate observed at approximately 300 K for DUT-49, which is discussed in Section 3.3. Given that NTE in MOF materials is considered a phonon-mediated phenomenon 35 it is unsurprising that NTE persists over this entire temperature range. Notably, this has also been experimentally observed for MOF-5, which shows NTE from 80-500 K. 36 The associated volumetric thermal expansion coefficients were also calculated over the entire temperature range, presented in Table 2, which provides a comparable metric for these materials. We find that a majority of the investigated MOFs show NTE coefficients > 50 MK −1 , and MOF-399 has an amazingly high coefficient of 87.6 MK −1 . There are very few materials with "colossal" thermal expansion, such as Ag 3 [Co(CN) 6 ], which is defined by 37 Similarly, the linear thermal expansion coefficient (α L ), defined in Equation 2, can be computed to provide a measure of the expansion rate with respect to each cell length.
Careful analysis of the linear contraction of each cell parameter, Table 2, shows that most of the investigated frameworks exhibit an isotropic contraction with little variation in cell shape even though no symmetry constrains were applied. However, the only framework considered in this study that has a triclinic unit cell, MOF-210, shows a distinct difference in linear contraction in the a-direction producing an anisotropic response.  Interestingly, there appears to be a correlation between the porosity features of the MOFs and the magnitude of NTE, Figure 3, suggesting the use of large pores and high levels of porosity as one way to design new and exceptional NTE materials. This was previously observed for studies of NTE in a series of isoreticular MOFs (pcu topology), 9 however, here we find it occurs across different topologies, although the network certainly influences NTE, as observed for the comparison of DUT-49 and PCN-68 and reported by Bouëssel du Bourg et al. for other MOF structures. 38 This relationship mirrors previous studies in Cd(CN) 2 frameworks, which showed that guest molecules within the framework can reduce NTE and even produce positive expansion. 39 Recently, Schneider et al. reported on a similar phenomenon in the MOF HKUST-1, whereby infiltration with TCNQ reduces NTE. 40 Large empty cavities provide the open space for transverse vibrational modes often responsible for NTE behavior. We suggest that large empty cavities can accommodate larger amplitude vibrations, thereby producing a greater degree of NTE.
The magnitude of NTE is also linked to the mechanical properties of the material through Equation 3, which can be derived from Equation 1, 41 where S is the entropy and B is the bulk modulus.
There is also a relationship between the porosity of a material and the bulk moduli, as reported previously in ceramics and zeolite materials. 42,43 Often these highly porous materials have extremely small bulk modulus, 4.97 and 2.07 GPa for DUT-60 and MOF-399, 14 respectively, which are very low for crystalline solids. Thus, the unique softness of mesoporous MOFs may also be responsible for the remarkable contraction observed.

Origins of expansion
The phenomenon of NTE in MOFs and related network materials is generally a result of phonon vibrations. 44 Specifically, NTE results from the often complex combination of transverse vibrational modes throughout the periodic lattice and rigid unit modes. 45,46 To confirm the presence of transverse phonon modes in the investigated MOFs, we computed the normal modes of vibration and subsequently generated total phonon density of states for MOF-399, DUT-60 and NU-110, the three best performing frameworks, using finite-displacements and the harmonic approximation. 34 The total phonon density of states are displayed in Figure 4a.
The MOF structures show several phonon modes with low frequency, <2 THz, which often give rise to NTE. Phonon modes of this frequency are attributed to collective vibrational modes, which can be visualized using the Molden 47,48 software package with the normal mode data files available online from our data repository at https://github.com/jackevansadl/ supp-data. For DUT-60 we can confirm these low-energy vibrational modes correspond to "trampoline"-like vibrations of the ligands that are also NTE-contributing modes in MOF-5. 49 Phonons that soften in response to contraction of the lattice have negative Grüneisen parameters, which directly contribute to the NTE observed in the material. 50 Please note, the softening described here relates to phonons shifting to lower frequencies. A comparison of the phonon density of states at different volumes, corresponding to the average MOF unit cell volume at 40 and 300 K highlight which of these low-energy modes give rise to NTE (Figure 4b). For each MOF, we find one or several peaks that display softening upon contraction, which contribute to the NTE observed in these systems.

Intrinsic ligand response
The analysis of phonon modes highlighted the contribution of "trampoline"-like dynamics of the large ligands to NTE in mesoporous MOFs. To decouple the effect of these transverse vibrational modes from other contributions to NTE, such as rigid unit modes, constrained dynamics of the largest molecular subunits of the MOFs (for MOFs constructed from two ligands the largest ligand was chosen) were investigated, Figure 5a.
The ligands show considerable contraction following the increase in temperature. DUT-6, the smallest ligand tested, showed less than 1 % contraction, while the largest ligand, NU-110, showed greater than 2 % contraction in length. This investigation clearly demonstrates the intrinsic contraction from constrained molecular units.
Interestingly, we find the magnitude of correlation is not entirely correlated with the size of the ligand. The largest molecular subunit of DUT-60, which has an average radius of 12 Å at 40 K, shows a similar amount of contraction as NU-110, with an average radius 20 Å at 40 K. This suggests that the chemical constituents of the subunit strongly influences the amount of contraction. The ligand of DUT-60 and MOF-399 is comprised entirely of phenyl groups, whereas NU-110 contains additional alkyne groups to produce a larger ligand. It appears the alkyne functionality produces less contraction for a larger molecular unit.
Furthermore, simulations of the DUT-49 molecular unit demonstrates a change in contraction rate at approximately 300 K, which was also observed for the entire framework as discussed in Section 3.1. This demonstrates how the individual dynamics of the ligands in the DUT-49 system can be transferred directly to the properties of the total framework. This direct transfer of dynamics may be related to the unique flexibility exhibited by DUT-49. 51,52

Mechanical consequences of negative thermal expansion
The maximum contraction of 3 % observed in these systems has a minor effect on their porosity. However, as reported by previous studies, this contraction has an important effect on the mechanical properties. For example, some MOF materials have shown pressure-induced softening before the point of instability. 38,53 To investigate the effect that temperature-induced contraction has on the mechanical stability of these materials, cell fluctuations from the molecular dynamics simulations were analyzed to compute the elasticity tensor, C i j . This elastic tensor contains key . The values at 0 K correspond to density functional theory simulations reported previously. 14 quantities that characterize the mechanical behavior of the structures in the elastic regime. For example, the directional Young's modulus is derived from the elasticity tensor. We found that all materials show temperature-induced softening, as illustrated in Figure 6a for MOF-399, which demonstrates considerable softening of the Young's modulus.
In particular, the eigenvalues of the C i j matrix are required to be positive for the structure to be elastically stable. 54 The maximum and minimum values of these eigenvalues indicate the most rigid and softest deformation mode. As we do not impose symmetry constraints, these are useful metrics in non-cubic crystals. 55 Notably, increasing temperature was found to soften the stiffest deformation mode, however, it does not change the softest mode, displayed for each MOF in the ESI. For example, MOF-399 shows a 50 % decrease in the stiffest deformation mode (Figure 6b). Therefore, NTE is of great consequence to the mechanical stability of these often fragile systems. Although MOF-399 has the potential to demonstrate record-breaking porosity, this has yet to be reported experimentally, suggesting mechanical or chemical fragility. Previously, we employed density functional theory calculations to determine the maximum eigenvalue of the stiffness matrix of this material, which was found to be 6.33 GPa. 14 However, the finite temperature calculations shown here demon-strate that at room temperature this material will be significantly less mechanically robust than first thought. We note that certain MOFs can be more resilient to temperature, such as PCN-68 which shows little variation in mechanical properties over the temperature range.
Recently, the decrease in stiffness with increasing temperature was experimentally reported for the MOF HKUST-1, which shows a decrease in Young's modulus in the temperature range 25-100 • C. 56 This provides further design rules to pushing the limits of ultrahigh porosity materials. Namely, this temperature-induced softening should be considered to produce the most mechanically robust porous materials over the relevant experimental temperature range.

Conclusions
We used molecular simulation to investigate NTE for a number of mesoporous MOFs. Given that MOFs have framework structures with large internal voids, perhaps unsurprisingly, they all showed contraction following a temperature increase, or NTE. The MOFs demonstrate impressive values of NTE > 50 MK −1 , but more importantly, the magnitude of contraction for the temperature range 40-400 K is extremely large, with MOF-399 exhibiting more than 3 % volume contraction.
The origin of NTE in DUT-60, MOF-399 and NU-110 was investigated by calculation of phonon modes. We identified a number of low-energy phonon modes that contribute to the NTE observed in these systems. In DUT-60 these are observed to correspond to transverse "trampoline"-like vibrations of the ligands within the lattice. Subsequently, we considered temperature-induced structure changes intrinsic to the large molecular subunits, which construct the studied MOFs, and found that they alone produce considerable NTE.
Finally, the effect of NTE on the mechanical properties was tested. There was significant softening of the most rigid modes by NTE, and this softening represents an important consideration for the mechanical stability and working temperature range of these highly porous materials.
This work demonstrates the potential for large NTE exhibited by structures with high porosity. Futhermore, the extent of NTE in these materials highlights another important characteristic to consider when designing ultrahigh porosity materials. For example, to design new materials with record-breaking surface area, NTE must be limited to ensure that a robust and experimentally realizable material is produced.

Conflicts of interest
There are no conflicts to declare.