Jessica K.
Bristow
^{a},
Jonathan M.
Skelton
^{a},
Katrine L.
Svane
^{a},
Aron
Walsh
*^{ab} and
Julian D.
Gale
*^{c}
^{a}Centre for Sustainable Chemical Technologies and Department of Chemistry, University of Bath, Claverton Down, Bath, BA2 7AY, UK
^{b}Department of Materials, Imperial College London, Exhibition Road, London SW7 2AZ, UK. E-mail: a.walsh@imperial.ac.uk; Tel: +44 (0)1225 385432
^{c}Curtin Institute for Computation, Department of Chemistry, Curtin University, PO Box U1987, Perth, WA 6845, Australia. E-mail: j.gale@curtin.edu.au
First published on 30th September 2016
We report the development of a forcefield capable of reproducing accurate lattice dynamics of metal–organic frameworks. Phonon spectra, thermodynamic and mechanical properties, such as free energies, heat capacities and bulk moduli, are calculated using the quasi-harmonic approximation to account for anharmonic behaviour due to thermal expansion. Comparison to density functional theory calculations of properties such as Grüneisen parameters, bulk moduli and thermal expansion supports the accuracy of the derived forcefield model. Material properties are also reported in a full analysis of the lattice dynamics of an initial subset of structures including: MOF-5, IRMOF-10, UiO-66, UiO-67, NOTT-300, MIL-125, MOF-74 and MOF-650.
A current dilemma for the development of forcefields for MOFs is the choice between transferability and accuracy. Large scale screening procedures offer a powerful tool for guiding experimental work but often involve multiple approximations, creating a level of uncertainty when used to calculate complex properties. Ab initio derived forcefields for prediction of charges and force constants offer accurate reproduction of the properties of individual or families of MOFs with similar topologies, but lack transferability. However, one must also consider that the purpose of using a forcefield for materials with large primitive unit cells is to remove the need to conduct expensive higher level calculations in the first place. A final consideration is associated with the diversity of MOFs. Due to differing compositions and topologies, the nature of the interaction between metal and ligand ranges between ionic and covalent. To derive one forcefield capable of reproducing all such interactions is a challenging task.
A vast number of forcefields for MOFs have already been reported; here we give a summary of the most prevalent and recent developments as a brief but broad overview. When parameterising a forcefield for MOFs there are many different approaches one can take. Firstly, existing transferable forcefields, such as CHARMM,^{1,2} MM3,^{3,4} GAFF^{5} and UFF,^{6} which have been extensively derived for common organic and, to a lesser extent, inorganic compounds, offer an abundant source of reasonable parameters for the individual components of MOFs. It is therefore only the interaction between metal and ligand for which additional potential parameters must be derived. UFF4MOF^{7} is an example of such a forcefield, reported by Addicoat et al., where the forcefield is an extension to the Universal forcefield (UFF). The UFF consists of multiple parameters capable of adequately reproducing the structures of organic molecules and inorganic clusters with little specific fitting. Indeed its incorporation in many user-friendly visualisation programs has increased the popularity of the UFF forcefield over many other more specifically parameterised transferable forcefields such as CHARMM. UFF4MOF employs additional atom types with parameterised valence coordination, equilibrium bond distances, effective charges and bond angles of the UFF to achieve a more accurate reproduction of the structure of clusters and periodic models of different MOFs. A more specific modification of UFF to reproduce the interaction of IRMOF-10 with CO_{2} was reported by Borycz et al.^{8}
Another approach is to derive bonding force constants and partial charges for each individual MOF using DFT calculations. MOF-FF is an example of a DFT-derived forcefield and was developed by Schmid et al.^{9} MOF-FF is capable of reproducing the structure and properties of many MOF topologies with initial parameterisation using the MM3 forcefield. The forcefield is then further fitted for individual MOFs based on data obtained from first principles calculations. The transferability for many families of MOFs has therefore not been extensively tested. Quick-FF was published by Van Speybroeck et al.^{10} and offers a method for rapid application of a forcefield to a MOF based on force constants extracted from a first principles Hessian. Quick-FF is still based on the MM3 non-bonding functional form but bonding parameters must be input by the user for the MOF of interest. The application of Quick-FF to MOF-5 and MIL-53 was reported, and was shown to reproduce both the structural parameters and the breathing behaviour of MIL-53.^{10} Here derivation of charges and bonding parameters remains dependent on first principles calculations. BTW-FF, our own previously reported MOF forcefield, which has been shown to reproduce the structure and properties of many different MOF topologies, also used partial charges obtained from DFT, using a Bader analysis of the charge density of each MOF considered.^{11}
Finally, large scale screening procedures can be based on a primitive mathematical description of the bonding in a framework or using a transferable forcefield, such as UFF, for all structures. High-throughput computational screening offers a valuable method for an approximate initial analysis of properties such as volumetric gas uptake or surface area. Screening also offers a means of structure prediction for hypothetical frameworks based on possible bonding considerations.^{12–15}
Comparisons of forcefields for simple MOFs, such as MOF-5, are now standard practice and offer little evaluation of the difficultly in derivation of the forcefield. Furthermore, the transferability of parameters across variations in the topology and crystallinity are rarely extensively tested. In particular, a recent focus is on defects and disorder in MOFs. Alterations of the frameworks must still render accurate and reliable results with the same forcefield parameters.
We highlight the calculation of phonon properties to be a relatively sparsely populated area of study for MOF forcefields; indeed, many forcefields, including previous work of ours, have merely compared vibrational Γ-point frequencies and plotted IR spectra. Phonon properties are critical for the analysis of dynamic stability, particularly if soft-modes are present, energetic stability (via free energies) and finite-temperature properties.^{16–20} Calculating phonon properties, such as phonon dispersion curves, with DFT is costly, and is only affordable and practical for specific MOFs of interest. Routine screening of phonon properties for large numbers of MOFs is currently only affordable with forcefield methods. One must also remember that experimental structures are often determined using X-ray diffraction. This can lead to inaccurate structure refinement of hydrogen positions and assignment of space groups that represent average structures. For some MOFs this may lead to loss of information regarding subtle distortions, such as non-planarity of carboxylate groups. Analysing the phonon stability of the structure, particularly for DFT calculations, is expensive and often avoided during electronic structure analysis. A transferable forcefield, not specifically fitted for each individual MOF, may allow small distortions of a framework to be identified, though the extent to which this is possible may depend on whether polarisability is included in the model or not. Optimising MOFs with a forcefield prior to using DFT may therefore be beneficial.
A final note regarding the importance of phonon property calculations is that MOF forcefields are often fitted at 0 K to an experimental structure determined at room temperature. Furthermore, temperature dependent properties (such as bulk moduli) from a 0 K optimisation are also often compared to room temperature experiments. Incorporating the consideration of temperature through free energy minimisation is a desirable alternative solution to having to make such approximate comparisons.
Here we present a new forcefield derived with the intention to bridge the gap between accuracy and transferability, while also incorporating an extensive analysis of phonon properties. The forcefield, named VMOF (vibrational metal–organic framework), is derived with the intention to be transferable, accessible and accurate when reproducing the structure and dynamical properties of MOFs. VMOF is a development of our previously reported BTW-FF forcefield for MOFs. In this paper we report the foundations and derivation of VMOF, along with a comparison of initial structure parameters and mechanical properties calculated for a range of MOFs. The main focus of the paper is then on discussing the accuracy of the forcefield for reproducing phonon properties obtained from DFT. We report densities of states (DOS), infra-red (IR) spectra and temperature dependent thermodynamic properties, such as free energies, vibrational entropy, and constant volume heat capacities. We show this new forcefield to be capable of accurately reproducing properties for an initial subset of MOFs. Finally, we perform quasi-harmonic calculations that, to the best of our knowledge, have not been previously reported for the given structures, and report additional temperature dependent structural properties.
(1) |
The MM3 Buckingham functional (eqn (1)) form consists of defined constants, A (1.84 × 10^{5}), B (12) and C (2.25), and was proposed as a “softer” energy function to that used in MM2. MM3 has been shown to accurately reproduce hydrogen and carbon positions in many aromatic compounds.^{3,4} The two remaining parameters per atom type, epsilon and the van der Waals radius, were fitted to reproduce phonon stable metal oxide structures by deriving these parameters for both the metal and the inorganic oxygen. The reference to inorganic oxygen here describes oxygen atoms that are coordinated only to metal atoms. Epsilon and van der Waals terms for the carboxylate oxygen and hydroxyl oxygen atom types were fitted in a relaxed fitting procedure to reproduce structural and mechanical properties of all the MOFs being tested. A further feature is the use of combination rules, where; ε_{ij} = and d_{ij} = , to reduce the number of parameters that require fitting, thus increasing the transferability. The long-range cut-off of the MM3 Buckingham potentials was set to 12 Å.
The derivation of parameters for the organic ligands was considered separately. Intramolecular bonding parameters for the ligands are taken directly from the CHARMM library and are implemented as harmonic functions. We consider the intramolecular bonding parameters between neighbouring atoms, angles between three connected neighbours and torsions between four connected atoms. A small modification was made to the CHARMM parameters for the 4-body torsion across the carboxylate head of all ligands considered, the derivation of which will be discussed later.
The total energy (U) can be written as;
The charges of the ligands are derived within GULP using the charge equilibration model of Gasteiger,^{27,28} while formal charges are used for the metal nodes and inorganic oxygen atoms. Gasteiger charges were selected since the charges are geometry independent and depend only on connectivity. Whilst other charge equilibration schemes suffer from charge delocalisation errors, Gasteiger charges do not.^{28} Initial charge parameterisation involved taking the average charge of each atom type in a subset of common MOF ligands including; 1,4-dicarboxylate (BDC), 1,3,5-tricarboxylate (BTC), 4,4′-biphenyl dicarboxylate (BPDC), 2,6-azulenedicarboxylate (AZ), 4,4′-biphenyl tricarboxylate (BPTC) and 2,7-pyrene-dicarboxylate (PDC). Once derived for a specific atom type, the same charges are used for all the structures modelled. Although Gasteiger charges are approximate, the advantage is that they are efficient to compute and well suited for our purpose.
The quasi-harmonic approximation^{31} allows a greater number of thermodynamic properties to be computed along with their temperature-dependence. The practical reality of the quasi-harmonic approximation is to minimise the internal energy at constant volume at a given number of lattice expansions and contractions away from the global minimum structure. An equation of state is then fitted across the calculated temperature-dependent Helmholtz free energies, which changes the minimum free energy volume according to a defined temperature. The temperature dependence of the phonon frequencies is then expressed in terms of Gibbs free energy. The theory is still dependent on calculations conducted under the harmonic approximation, but the consequent volume dependence of the vibrational frequencies is an anharmonic effect. In general it is possible to use the analytical derivatives of the free energy with respect to strain to optimise all cell parameters independently within the ZSISA^{32} (zero static internal stress) approximation. While this can be routinely achieved using some forcefield implementations,^{33} the requirement to compute third-order derivatives makes this particularly demanding with quantum mechanical approaches.
Throughout the calculation of thermodynamic properties a consistent Brillouin-zone phonon sampling (q-point) mesh of 32 × 32 × 32 was used, with a symmetry tolerance of 10^{−3} Å for determining the space group symmetry during atomic displacement generation.
To remove imaginary modes, observed only during compression in the DFT calculations, the corresponding Hessian eigenvectors were used to displace the atoms, thereby lowering the symmetry, and the structure relaxed accordingly. Removing imaginary modes with the forcefield required the use of a rational function optimisation approach to ensure that the Hessian has the correct final structure after optimisation.
a | b | c | M–O | C _{11} | C _{12} | C _{13} | C _{33} | C _{44} | C _{55} | B _{0} | |
---|---|---|---|---|---|---|---|---|---|---|---|
ZnO | |||||||||||
Exp.^{34} | 3.250 | 3.250 | 5.207 | 1.992 | 209.6 | 121.1 | 105.1 | 210.9 | — | 42.5 | 183.0 |
DFT | 3.225 | 3.225 | 5.211 | 1.972 | 187.8 | 123.1 | 108.6 | 209.0 | — | 27.9 | 143.4 |
% error | (0.77) | (0.77) | (0.08) | ||||||||
VMOF | 3.219 | 3.219 | 5.014 | 1.967 | 242.6 | 108.2 | 100.7 | 199.2 | — | 78.4 | 143.2 |
% error | (0.95) | (0.95) | (3.71) | ||||||||
ZrO_{2} | |||||||||||
Exp.^{35} | 5.070 | 5.070 | 5.070 | 2.195 | 533.5 | 97.9 | — | — | 64.3 | — | 243.7 |
DFT | 5.064 | 5.064 | 5.064 | 2.193 | 575.6 | 139.8 | — | — | 80.1 | — | 279.3 |
% error | (0.12) | (0.12) | (0.12) | ||||||||
VMOF | 5.123 | 5.123 | 5.123 | 2.218 | 630.8 | 131.2 | — | — | 125.7 | — | 297.7 |
% error | (1.05) | (1.05) | (1.05) | ||||||||
TiO_{2} | |||||||||||
Exp.^{36} | 4.594 | 4.594 | 2.959 | 1.980, 1.949 | 366.0 | 225.0 | — | — | 189.0 | — | 282.0 |
DFT | 4.585 | 4.585 | 2.941 | 1.971, 1.944 | 274.5 | 208.8 | — | — | 236.3 | — | 251.9 |
% error | (0.20) | (0.20) | (0.61) | ||||||||
VMOF | 4.414 | 4.414 | 3.168 | 1.980, 1.934 | 362.4 | 337.2 | — | — | 213.7 | — | 335.1 |
% error | (3.92) | (3.92) | (7.06) | ||||||||
Al_{2}O_{3} | |||||||||||
Exp.^{37} | 4.764 | 4.764 | 13.001 | 1.858 | 497.3 | 162.8 | 116.0 | 500.9 | 146.8 | — | 240.0 |
DFT | 4.755 | 4.755 | 12.962 | 1.966 | 482.2 | 157.0 | 121.0 | 493.3 | 162.6 | — | 269.4 |
% error | (0.19) | (0.19) | (0.30) | ||||||||
VMOF | 4.870 | 4.870 | 12.899 | 1.848 | 564.3 | 224.3 | 152.2 | 463.3 | 123.9 | — | 291.3 |
% error | (2.23) | (2.23) | (0.78) |
The forcefield model that we have chosen to adopt to increase transferability involves formal charges at the metal node. The unit cell parameters and elastic constants of the metal oxides reproduced by the forcefield are generally reasonable (Table 1), though with a tendency to overestimate the hardness of materials. Formal charge models usually include the shell model for polarisation of the oxide ions, which can effectively soften the mechanical properties, but we maintain the use of a rigid ion model for consistency with the CHARMM parameters for the organic ligands.
A final comparison of the mechanical properties between experiment and DFT further highlights the variation in calculated values of the elastic constants and bulk moduli of the metal oxides (Table 1). The mechanical properties calculated with VMOF remain within a reasonable range with both experiment and DFT values.
Fig. 1 Structures of MOF-5, IRMOF-10, UiO-66, UiO-67, MOF-650, MOF-74, MIL-125 and NOTT-300. Metal polyhedra are highlighted for Zn (silver), Zr (green), Ti (blue) and Al (grey), with atoms coloured black for carbon, white for hydrogen and red for oxygen. Compositions and symmetries are given in Table 2. |
Prior to the calculation of thermodynamic properties, several observations regarding the vibrational stability of IRMOF-10 and MIL-125 were made during optimisation. Firstly, for IRMOF-10 we initially calculated a significant number of imaginary vibrational modes. To relax the structure into the ground-state and remove all sources of instability, all imaginary modes were simultaneously relaxed following an initial displacement along the corresponding phonon eigenvectors. The final structure was re-converged and no imaginary modes were found. We attribute the initial structural instability to the BPDC ligand; following optimisation we observed a rotation about the central C–C bond connecting the two aromatic rings. The final torsion across this bond was 30–31° between rings (Fig. 2). We propose that the planar experimental configuration may be a thermally averaged structure, and that the true ground-state actually involves twisted ligands. For IRMOF-10 we calculate the structure with twisted ligands to be 0.275 eV (0.046 eV per ligand) more stable than the planar structure. UiO-67 is formed of the same BPDC ligands, which are experimentally characterised as twisted with near identical angles to those in IRMOF-10 between torsion planes. Electronic structure calculations were recently reported by Hemelsoet et al., highlighting the flexibility of the BPDC ligand in UiO-67. This study reported the difference in relative occurrence in torsion angles between the aromatic rings between 0–90° during a molecular dynamics simulation. As IRMOF-10 is formed of weaker intermolecular interactions, an increased flexibility of torsion angles would be expected.^{41} Furthermore, we calculate the BPDC^{2−} ligand in the gas phase to possess a torsion angle of approximately 33° in its ground-state configuration (further details in the ESI†).
Fig. 2 Overlaid planar (blue) and twisted (black) 4,4′-BPDC ligands showing the change in geometry of the ligand after following the imaginary phonon modes in the initial structure of IRMOF-10. |
Similar structural instabilities were observed for MIL-125, which is reported as belonging to the tetragonal space group number 139. Initially we obtained 17 associated imaginary modes, and re-optimised to produce a structure with two imaginary modes. This structure was found to have a broken symmetry with the hydrogen on the hydroxyl groups flipped into the pore of the MOF rather than being held in the pore windows. The remaining two imaginary modes could not be removed with further optimisation, and the calculations became too expensive to continue. However, the lower symmetry structure was 0.407 eV per primitive unit cell more stable, with no external pressure on the cell.
The selected MOFs studied here were chosen to ensure a variety of different topologies and bonding interactions were present, thus testing the broad applicability of the forcefield. The metal cations and ligands comprising these MOFs are given, along with the experimentally determined space groups, in Table 2.
Name | Metal | Ligand | Space group |
---|---|---|---|
MOF-5 (IRMOF-1) | Zn^{2+} | 1,4-Benzene dicarboxylate | 225 (Fmm) |
IRMOF-10 | Zn^{2+} | 4,4′-Biphenyl dicarboxylate | 225 (Fmm) |
MOF-650 | Zn^{2+} | 2,6-Azulene dicarboxylate | 225 (Fmm) |
UiO-66 | Zr^{4+} | 1,4-Benzene dicarboxylate | 225 (Fmm) |
UiO-67 | Zr^{4+} | 4,4′-Biphenyl dicarboxylate | 225 (Fmm) |
MIL-125 | Ti^{4+} | 1,4-Benzene dicarboxylate | 139 (I4/mmm) |
MOF-74 | Zn^{2+} | 2,5-Dihydroxyterephthalic acid | 2 (P) |
NOTT-300 | Al^{3+} | 3,3′,5,5′-Biphenyltetracarboxylic acid | 98 (I4_{1}22) |
A comparison of optimised unit cell parameters between DFT and forcefield methods are given in Table 3. All structures are reproduced by the forcefield with low errors on the unit cell parameters, thus demonstrating the accuracy of the forcefield and its ability to reproduce different structural features of MOFs, despite the simplicity of its derivation.
Name | DFT unit cell parameters | VMOF unit cell parameters |
---|---|---|
MOF-5 (IRMOF-1)(Zn^{2+}) | 25.894, 25.894, 25.894 | 25.935 (0.16), 25.935 (0.16), 25.935 (0.16) |
IRMOF-10 (Zn^{2+}) | 34.385, 34.385, 34.385 | 34.417 (0.09), 34.417 (0.09), 34.417 (0.09) |
MOF-650 (Zn^{2+}) | 30.695, 30.695, 30.695 | 30.766 (0.23), 30.766 (0.23), 30.766 (0.23) |
MOF-74 (Zn^{2+}) | 6.740, 15.142, 15.142 | 6.764 (0.35), 15.031 (0.73), 15.031 (0.73) |
UiO-66 (Zr^{4+}) | 20.798, 20.798, 20.798 | 20.909 (0.53), 20.909 (0.53), 20.909 (0.53) |
UiO-67 (Zr^{4+}) | 27.094, 27.094, 27.094 | 26.878 (0.80), 26.878 (0.80), 26.878 (0.80) |
MIL-125 (Ti^{4+}) | 18.852, 18.843, 17.921 | 18.859 (0.01), 18.859 (0.08), 18.043 (0.68) |
NOTT-300 (Al^{3+}) | 14.836, 14.836, 11.871 | 14.862 (0.18), 14.862 (0.18), 11.500 (3.23) |
Name (metal) | DFT (GPa) | VMOF (GPa) |
---|---|---|
MOF-5/IRMOF-1 (Zn^{2+}) | 16.9 | 8.8 |
IRMOF-10 (Zn^{2+}) | 8.6 | 5.1 |
MOF-650 (Zn^{2+}) | 12.5 | 6.8 |
UiO-66 (Zr^{4+}) | 40.4 | 19.0 |
UiO-67 (Zr^{4+}) | 21.9 | 11.7 |
MIL-125 (Ti^{4+}) | 25.1 | 18.5 |
MOF-74 (Zn^{2+}) | 28.1 | 14.9 |
NOTT-300 (Al^{3+}) | 47.8 | 25.2 |
The forcefield predicts smaller bulk moduli than DFT. The discrepancy in bulk moduli between DFT and forcefield is likely to be due to an inaccurate description of the long-range dispersion interactions. The converged unit cell volumes calculated with DFT are generally smaller than those predicted with the forcefield, leading to larger bulk moduli. We do, however, highlight that experiment often finds MOFs to have a softer bulk modulus than those predicted with electronic structure methods. Yot et al. reported experimental bulk moduli, as measured using high pressure XRD methods, for UiO-66 and MIL-125 to be 17.0 and 10.0 GPa, respectively, which are closer to the forcefield values than those calculated by DFT.^{42} The trend in mechanical strength for each framework is the same between the two methods, with increasing mechanical strength from IRMOF-10 < MOF-74 < MOF-650 < MOF-5 < UiO-67 < MIL-125 < UiO-66 < NOTT-300.
Higher charged cations form primarily ionic interactions between metal and ligand. It can therefore be rationalised that UiO-66, UiO-67 and MIL-125 would possess stiffer bulk moduli (greater resistance to compression) than the Zn-isoreticular MOFs, which possess large internal voids and weaker van der Waals interactions between Zn^{2+} and organic ligands.
The decomposition of the DOS into elemental contributions is shown in the partial DOS plots (Fig. 5 and 6). Several features are evident in both spectra calculated with DFT and FF methods that are chemically well established. Firstly, the range of modes involving the metal cations all occur at low frequencies (<500 cm^{−1}). Also observed in this region for all MOFs is a small contribution from C, H and O from the rocking motions of the ligands. At finite temperature it is these low frequency modes that are populated, and therefore control the MOF dynamics (e.g. the shape of the thermal ellipsoid). Therefore, the motions of the MOF will occur primarily at the metal nodes, as well as subtle rotations at the MOF-ligand connections. In the mid-frequency range, between 1300–1500 cm^{−1}, the high contribution of C and O to the density of states is due to motions of the asymmetric and symmetric C_{carb}–O_{carb} stretches of the carboxylate groups in the MOFs. Finally, modes above 3000 cm^{−1} are associated with C–H and O–H stretches at the ligand and within some metal nodes, respectively. The important conclusion from the density of states plots is that we see good agreement between DFT and forcefield methods, suggesting the vibrational properties of the MOFs are well reproduced by VMOF.
It is common to characterise the vibrational properties of a material by assigning specific IR frequencies. However, we demonstrate that a vast amount of vibrational information is not accounted for by doing this for MOFs. The comparison of DOS and IR spectra show the difference in detail and highlight the importance of considering vibrational modes that are not IR active when parameterising a forcefield. The DOS also shows the significant number of soft vibrational modes that MOFs possess, which give rise to structural instability with temperature and pressure. We note a shift of the frequencies of the C–H stretch between the two methods. As the C–H stretch occurs as one of the highest frequency modes, it is contributing the most to the zero-point vibrational energy and is likely to be the biggest contribution to the C–H stretch calculated error between methods. The reliability of the forcefield is not likely to be affected by the disagreement in zero point energy between the methods for the study of complex processes such as phase changes. It is unclear if it is DFT or the FF model that contains the greatest error on the C–H stretch. The forcefield parameters for the C–H interaction remain unchanged from the CHARMM forcefield and therefore have not been derived specifically for BDC incorporated into a MOF. On the other hand, a scaling factor is often used on the vibrational frequencies in DFT calculations that would have the greatest effect on the stretching mode of the C–H bond. Such scaling factors can correct for anharmonicity, while forcefield parameters can be derived explicitly to reproduce experimental anharmonic frequencies.
The constant volume (Helmholtz) free energy, vibrational entropy and constant-volume heat capacities as a function of temperature are plotted for all MOFs considered with both first principles and forcefield methods (Fig. 7). The forcefield is shown to reproduce the calculated thermodynamic properties from DFT very well. Little deviation across all structures is observed, further supporting that the forcefield can accurately reproduce the vibrational properties of the subset of MOFs studied.
Whilst conducting tests for the quasi-harmonic approximation we noted a particular sensitivity of the fitted forcefield parameters to the free energy equation of state with framework compressions. Consequently, calculated properties varied depending on the volume sampled for the expansions and contractions away from the local minimum structure. For MOF-5 the initial volume sampled was ±3% away from the energy minimum in 0.05% steps. We observed, with both DFT and forcefield methods, that beyond approximately 2% compression multiple imaginary vibrational modes emerged. When following these imaginary modes, a subtle twist at the carboxylate head relative to the benzene ring (approximately 3–5°) was observed. The same observation was made for IRMOF-10, suggesting that with compression “softer” MOFs undergo this subtle rotation of the benzene rings leaving a non-planar torsion between the ring and carboxylate heads of the ligand. We note that we could remove imaginary modes for all structures during compression with the forcefield, but the same process with electronic structure methods became too expensive for IRMOF-10, which possessed 2 imaginary modes at the highest compression of 3%.
To investigate this phenomenon further we modelled the same rotation of BDC in MOF-5 away from the initial planar 0° structure. Calculations were conducted on the periodic systems with the VMOF forcefield. To model the same rotation with DFT, we cut a representative cluster with the chemical formula Zn_{4}O(BDC-H_{2}). The PBE0 functional was used in the NWChem program with the Dunning correlation consistent cc-pVTZ basis sets.^{39,40,44} The rotation leaves the carboxylate heads in the same position in a direct bonding interaction with the metal centres, and only moves the benzene ring; following the rotation, there is a non-planar torsion between ring and carboxylate head (see ESI†). Rotating the ligand in this manner in MOF-5 is shown to cost little energy with both DFT and VMOF up to 10° (Fig. 8). The potential energy surface is shallow, with a 5° change being comparable to k_{B}T. Therefore it can be expected that the group will be rotationally active at room temperature. We highlight the similarity in the calculated trend between DFT and VMOF for the modelled rotation, supporting our observations of the instability of the planar structure with compression with DFT.
Once the physical significance of structural distortion with compression was investigated, we were able to rationalise trends observed from quasi-harmonic calculations. Due to the expense of the quasi-harmonic approximation calculations with first principles, we report only forcefield thermodynamic properties for UiO-67, MIL-125, MOF-650 and MOF-74 (Fig. 10).
There is a limited amount of data in the literature for experimental verification of the predicted thermal expansion profiles. Our calculated thermal expansion profile for MOF-5 however resembles the profile measured by Zhou et al., where an increasing thermal expansion coefficient is observed above 100 K. At 300 K, Zhou et al. measure the thermal expansion coefficient of MOF-5 to be approximately 15 × 10^{−6} K^{−1}, the corresponding thermal expansion coefficient calculated with VMOF and DFT are 23 × 10^{−6} and 36 × 10^{−6} K^{−1}, respectively.^{45}
The use of the quasi-harmonic approximation allows the calculation of a more extensive range of properties, such as the Grüneisen parameter, bulk modulus, heat capacity and thermal expansion coefficient, as a function of temperature. For all structures we found that the QHA properties that showed the greatest sensitivity were the bulk modulus and thermal expansion. These properties are derived directly and indirectly from the curvature of the potential energy surface, respectively, and therefore exaggerate the difference between the methods. We highlight that VMOF, which was not specifically parameterised to reproduce negative thermal expansion of MOFs, manages to exhibit this phenomenon and accurately in comparison to DFT for MOF-5, IRMOF-10, UiO-66 and NOTT-300. Due to the expense of the quasi-harmonic approximation calculations with first principles, we report only forcefield thermodynamic properties for UiO-67, MIL-125, MOF-650 and MOF-74 (Fig. 10).
Whilst fitting an equation of state for IRMOF-10, we observed a sensitivity of the calculated properties with DFT to the extent of compression considered. Properties such as the bulk modulus, Grüneisen parameter and thermal expansion were calculated to vary significantly when including high-pressure points (i.e. large compressions). We also observed a significant shift in minimum of Helmholtz energy with change in volume at the defined temperatures. Such a shift and variation in calculated properties suggests the zero point energy contribution to have a large effect on the structure.
Several interesting features more specific to each structure are observed from the QHA calculations (Fig. 9 and 10). Firstly, a comparison of thermal expansion coefficients show the vast difference in structural changes with temperature. The “softer” Zn-isoreticular MOFs, such as MOF-5 and IRMOF-10 with both DFT and FF methods, are calculated to possess the largest negative thermal expansion coefficients at low temperature, suggesting the extent of thermal expansion reflects the mechanical strength of the materials. Indeed, the “hardest” MOFs, UiO-66 and NOTT-300, show little variation in thermodynamic properties with temperature. The temperature dependent bulk moduli are shown to differ between DFT and forcefield methods (Fig. 9), following the calculated trend in static bulk moduli. Specifically, VMOF yields softer mechanical properties at finite temperatures than electronic structure methods. We note that the trend and shape of each profile with temperature is reproduced by the forcefield, and that the temperature dependent Grüneisen parameter and heat capacity at constant pressure are reproduced accurately and appear unaffected by the discrepancy in bulk moduli. Due to the accuracy of the quasi-harmonic approximation being restricted to 1/2–2/3 the melting point temperature (depending on the material), we cannot model the full behaviour of the heat capacity at high temperatures. The temperature dependence of the thermodynamic properties of MOF-650 appears to show a significantly different trend to all other structures, which all show similar behaviour. MOF-650 has a large internal void with cell parameters exceeding 30 Å. The azulene ligand comprising MOF-650 is also rigid and would allow little structural flexibility when compared to IRMOF-10, which is of similar size. It is likely to be these two factors that result in positive thermal expansion and Grüneisen parameter at low temperature. Finally, the Grüneisen parameter shows a similar trend for each structure, excluding MOF-650. The increase in Grüneisen parameter with temperature reflects the increase in mechanical strength following contraction of the cell parameters.
Footnote |
† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6cp05106e |
This journal is © the Owner Societies 2016 |