A general forcefield for accurate phonon properties of metal–organic frameworks

Abstract

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.

---

1 Introduction

Metal–organic frameworks (MOFs), formed from metal cations and anionic organic molecules, are versatile materials with a range of functional properties. There are a vast number of possible candidate structures from known metals and ligands, which may self-assemble into a variety of 3D MOF structures. Developing a cheap, transferable and accurate method for initial property screening of potential MOF candidates for applications such as gas absorption, explosives detection and use in solar energy conversion is therefore desirable to reduce the cost and time of experimental work. Given the tractable, but still computationally expensive, nature of density functional theory (DFT) and higher levels of first principles theory, the development of forcefields capable of reproducing structural, mechanical and vibrational properties of MOFs would be highly advantageous to the large community of computational chemists interested in the thermoelastic properties of MOFs.

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⁵ and UFF,⁶ 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⁷ 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₂ was reported by Borycz et al.⁸

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.⁹ 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.¹⁰ 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.¹⁰ 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.¹¹

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.

2 Methodology

2.1 First principles reference calculations

Reference quantum mechanical calculations were conducted to have a standard method of validation of the new forcefield. The Vienna ab initio Simulation Package (VASP)²¹ code was used to perform Kohn–Sham density functional theory (DFT) calculations using the PBEsol exchange–correlation functional.²² The projector-augmented wave method²³ was used for the interaction between core and valence electrons of all atoms in the system. During optimisation, all forces were converged to values of less than 0.001 eV Å⁻¹ with a plane wave basis set cut-off of 600 eV. Γ-Point sampling of the Brillouin zone was considered sufficient for the MOFs owing to the unit cell dimensions of the systems, excluding MOF-74, which required a 4 × 4 × 4 k-point grid. A D3 van der Waals correction²⁴ was included and found to be necessary to remove phonon instabilities for some MOF structures. Reference calculations for the binary oxides: ZnO, Al₂O₃, TiO₂ and ZrO₂, were performed with the same convergence criteria in VASP, with the chosen polymorphs being wurtzite, corundum, rutile and baddelite, respectively. The k-point grid and plane wave cut-off were converged separately for each metal oxide, with final values being given in the ESI† along with the optimised unit cell parameters.

2.2 Forcefield calculations

VMOF was derived using the General Utility Lattice Program (GULP) code, which has extensive capabilities suited to both inorganic and organic materials.^25,26 VMOF considers the metal node and organic component as essentially separate entities interacting only by modified MM3 Buckingham potentials, in addition to the Coulomb terms.

The MM3 Buckingham functional (eqn (1)) form consists of defined constants, A (1.84 × 10⁵), 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;

where, k_r, k_θ and k_Ψ are interatomic force constants, r the distance between pairs of atoms, θ and Ψ are angles, q represents point charges and ε₀ is the vacuum permittivity. Note that the harmonic bonding terms in GULP possess a multiplication factor of ½, and so CHARMM force constants were appropriately scaled.

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.²⁸ 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.

2.3 Property calculations

2.3.1 Mechanical properties. Bulk moduli (B₀) were calculated from the relevant components of the elastic constant and compliance tensors, which were determined from the analytical second derivatives of the energy with respect to strain on the system. The elastic compliance tensor is just the inverse of the elastic constant tensor. The reported bulk moduli calculated with the forcefield follow the Hill convention, i.e. they are the averages of the Reuss and Voigt definitions. First principles bulk moduli for each structure were calculated at 0 K by fitting a third-order Birch–Murnaghan equation of state to energy-volume data calculated at a series of expansions and contractions about the equilibrium structure. For each expansion and contraction the atomic positions were optimised. First principles elastic constants were calculated using the finite differences method to construct the second derivatives of the energy with respect to the atomic positions. Only symmetry inequivalent displacements were considered.

2.3.2 Vibrational properties. IR frequencies and intensities, as well as DOS, were post-processed from both forcefield and first principles methods to apply a consistent broadening factor of 10 cm⁻¹ with a frequency sampling resolution of 1 cm⁻¹. First principles calculations of dynamical matrices were performed using density functional perturbation theory within VASP to obtain ε^∞ and phonopy was used to calculate the eigenvalues. IR frequencies and intensities were then plotted with identical resolution and broadening factors.

2.3.3 Lattice dynamics and thermodynamic properties. To calculate the thermodynamic properties, such as Helmholtz and Gibbs free energies, vibrational entropy and volumetric heat capacity, and the dependence with temperature, we use phonopy.^29,30 Phonopy is a python package for setting up post-processing finite-displacement phonon calculations that can be integrated with multiple first principles codes, and includes an extension for the quasi-harmonic approximation. Further details are given in the ESI.†

The quasi-harmonic approximation³¹ 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³² (zero static internal stress) approximation. While this can be routinely achieved using some forcefield implementations,³³ 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⁻³ Å 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.

3 Results

3.1 Forcefield fitting

3.1.1 Universal forcefield and metal oxide parameters. When considering the initial derivation of parameters for each metal, besides fitting the epsilon and van der Waals radius for each species, we have also examined the influence of changing the Universal MM3 constants. This was found to be beneficial to the overall quality of the results. After extensive testing, we found that the modification of the MM3 constants by changing the B parameter to 11.5 and C parameter to 2.55, reproduced DFT and experimental structural and mechanical properties of the MOFs and metal oxides more accurately. In support of the modification of the MM3 constants, we report structural and mechanical properties for the binary oxides with the original MM3 constants (see ESI†) and modified constants (Table 1) in comparison with experimental values. Using the original MM3 constants results in large errors for second order elastic properties such as the elastic constants and bulk moduli.

Table 1 Comparison of structural and mechanical properties of metal oxides between experiment (exp.) and the VMOF forcefield with modified MM3 parameters (VMOF). Given are the elastic constants (C_ij), bulk moduli (B₀), unit cell parameters and metal–oxygen (M–O) bond lengths. All lengths are reported in Å and elastic constants and bulk moduli are reported in GPa. Percentage errors in the cell parameters relative to the experimental values are given in parenthesis under each value reported. DFT values are calculated with the PBEsol functional

	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)
ZrO2
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)
TiO2
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)
Al2O3
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.

3.1.2 Ligand parameters. Force constants for the ligand parameters are taken directly from the CHARMM library. During the derivation of forcefield parameters for the metal nodes, there was one 4-body interaction with the ligand that showed a particular propensity for producing phonon unstable structures if varied. The torsion across the head of the carboxylate groups did not always remain planar during optimisation, and enforcing planarity by increasing the force constant across this bonding connection often rendered structures phonon unstable. We therefore calculated the force constant across this interaction in an isolated BDC²⁻ ligand in the gas phase. The PBE0 functional³⁸ was used in the NWChem program³⁹ with the Dunning correlation consistent cc-pVTZ basis sets⁴⁰ to fully relax the ligand with a 0 and 90° torsion of the carboxylate heads in relation to the aromatic ring (depicted in the ESI†). We calculate the energy difference between the two configurations to be 0.530 eV. To maintain transferability of the forcefield we assume little variation of this energy would occur across different aromatic dicarboxylate ligands. Therefore, it is this value that the force constant is fitted to for the torsion between O_carb–C_carb–C_benz–C_benz for all structures with these atom types.

3.2 Structural properties

Following the derivation of the forcefield parameters, mechanical and vibrational properties have been calculated for eight different MOFs (Fig. 1), representing a range of ligand and metal node types: MOF-5, IRMOF-10, UIO-66, UIO-67, MOF-650, MIL-125 and NOTT-300.

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.⁴¹ Furthermore, we calculate the BPDC²⁻ 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.

Table 2 MOFs modelled with the VMOF forcefield in this work. Given are the metal cations with the corresponding formal oxidation state and the ligands comprising the MOFs, along with the reported space group number and name

Name	Metal	Ligand	Space group
MOF-5 (IRMOF-1)	Zn^2+	1,4-Benzene dicarboxylate	225 (Fm3̄m)
IRMOF-10	Zn^2+	4,4′-Biphenyl dicarboxylate	225 (Fm3̄m)
MOF-650	Zn^2+	2,6-Azulene dicarboxylate	225 (Fm3̄m)
UiO-66	Zr^4+	1,4-Benzene dicarboxylate	225 (Fm3̄m)
UiO-67	Zr^4+	4,4′-Biphenyl dicarboxylate	225 (Fm3̄m)
MIL-125	Ti^4+	1,4-Benzene dicarboxylate	139 (I4/mmm)
MOF-74	Zn^2+	2,5-Dihydroxyterephthalic acid	2 (P1̄)
NOTT-300	Al^3+	3,3′,5,5′-Biphenyltetracarboxylic acid	98 (I4122)

---

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.

Table 3 Comparison of unit cell parameters (a, b, c with all values in Å) from the VMOF forcefield with those calculated using DFT (PBEsol functional); percentage errors of FF values compared to DFT are given in brackets

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)

---

3.3 Mechanical properties

Bulk moduli have been calculated with DFT and VMOF to compare the mechanical strength of the materials predicted with the two approaches (Table 4).

Table 4 Comparison of bulk moduli (GPa) obtained with the VMOF forcefield with those calculated with DFT. DFT bulk moduli were calculated with the PBEsol functional with D3 dispersion correction, by fitting calculated energy/volume curves to a Birch–Murnaghan equation of state with ±3% sampling away from the equilibrium volume

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.⁴² 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²⁺ and organic ligands.

3.4 Phonon properties

3.4.1 IR spectra. The first approach to assessing the ability of a forcefield to reproduce accurate vibrational properties is to calculate the phonon density of states and the associated IR spectra (weighted by the mode intensities) to ensure the fingerprint of vibrational modes is similar between DFT and forcefield methods. Good agreement is observed between vibrational IR spectra and DOS between DFT (Fig. 3) and the forcefield (Fig. 4). The plots highlight the stability of the modelled MOFs with both methods, but also the small deviation between the two sets of calculated DOS and IR spectra. The biggest discrepancy between VMOF and DFT IR spectra is in the fingerprint (lower frequency region). This is due to the metal–oxygen bond stretching modes, and since IR activity ∝ charge × displacement, the discrepancy is primarily due to the use of formal charges for the metal ions. Importantly, we highlight that the DOS spectra between DFT and forcefield remain comparable, which suggests the use of a formal charge model has had little effect on the forces.

Fig. 3 Overlaid IR spectra and phonon density of states (DOS) calculated with DFT (top, black) and FF (bottom, purple) methods, plotted between 0–3200 cm⁻¹ for MOF-5 (a), IRMOF-10 (b), MOF-650 (c) and MOF-74 (d). All spectra are normalised to lie between 0 and 1 and area the under DOS is shaded (grey) for clarity.

Fig. 4 Overlaid IR spectra and phonon density of states (DOS) calculated with DFT (top, black) and FF (bottom, purple) methods, plotted between 0–3900 cm⁻¹ for UiO-66 (a), UiO-67 (b), MIL-125 (c) and NOTT-300 (d). All spectra are normalised to lie between 0 and 1 and area under the DOS is shaded (grey) for clarity.

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⁻¹). 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⁻¹, 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⁻¹ 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.

Fig. 5 Partial phonon density of states (PDOS) projected according to the elemental contribution calculated with DFT (top) and FF (bottom) methods, plotted between 0–3200 cm⁻¹ for MOF-5 (a), IRMOF-10 (b), MOF-650 (c) and MOF-74 (d). All spectra are normalised to lie between 0 and 1.

Fig. 6 Partial phonon density of states (PDOS) decomposed according to the elemental contribution calculated with DFT (top) and FF (bottom) methods, plotted between 0–3900 cm⁻¹ for UiO-66 (a), UiO-67 (b), MIL-125 (c) and NOTT-300 (d). All spectra are normalised to lie between 0 and 1.

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.

3.4.2 Helmholtz free energy. Reproducing accurate relative free energies of a system with temperature is crucial for the prediction of thermodynamic processes such as phase changes and reaction energies for the formation of defects within frameworks, for example in our previous work investigating the “missing linker phenomenon” in UiO-66.⁴³

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.

Fig. 7 Comparison of calculated Helmholtz free energy (black), vibrational entropy (blue) and constant volume heat capacity (C_v) (red) of (a) MOF-5, (b) IRMOF-10, (c) MOF-650, (d) MOF-74, (e) UiO-66, (f) UiO-67, (g) MIL-125 and (h) NOTT-300 between DFT (symbol) and FF (line) methods.

3.4.3 Quasi-harmonic approximation properties. Fitting forcefield models to ground state properties of the equilibrium structure does not account for the variation in volume and unit cell shape with temperature. As a method for extending the harmonic model of vibrations, the quasi-harmonic approximation (QHA) allows anharmonicity associated with volume change (thermal expansion) to be considered when calculating structural and thermodynamic properties.

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₄O(BDC-H₂). 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_BT. 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.

Fig. 8 Energy cost per ligand when rotating all C–C–C–O torsions in the ligands by successive angle increments away from the initial planar 0° structure. Relative energies are reported for both DFT and VMOF.

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⁻⁶ K⁻¹, the corresponding thermal expansion coefficient calculated with VMOF and DFT are 23 × 10⁻⁶ and 36 × 10⁻⁶ K⁻¹, respectively.⁴⁵

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.

Fig. 9 Comparison of structural and thermodynamic properties as a function of temperature calculated with the quasi harmonic-approximation for MOF-5, IRMOF-10, UiO-66 and NOTT-300 using DFT and the VMOF forcefield (FF). (a) Linear thermal expansion coefficient (b) Grüneisen parameter (c) bulk modulus (d) heat capacity at constant pressure.

Fig. 10 Structural and thermodynamic properties as a function of temperature calculated with the quasi-harmonic approximation for MOF-650, MOF-74, MIL-125 and UiO-67 using the VMOF forcefield (FF). (a) Linear thermal expansion coefficient (b) Grüneisen parameter (c) bulk modulus (d) heat capacity at constant pressure.

4 Conclusions

A new transferable forcefield for metal–organic frameworks named VMOF has been parameterised to reproduce accurate lattice dynamics and phonon properties. Such a forcefield contributes greatly to the current extensive field of MOF forcefields as it is unique in the number of thermodynamic properties that can be accurately determined in a rapid and transferable manner. For an initial training set of MOFs including MOF-5, IRMOF-10, UiO-66, UiO-67, NOTT-300, MIL-125, MOF-650 and MOF-74 we calculate numerous thermodynamic properties including bulk modulus, free energies and constant volume heat capacities. We further conduct quasi-harmonic calculations and find excellent agreement in thermal expansion, bulk moduli, Grüneisen parameter and heat capacity with temperature between DFT and the newly parameterised forcefield. This now opens the way for the future high-throughput computational screening of materials vibrational properties for a wide range of MOFs.

Data access statement

The VMOF forcefield is available as a library file for GULP from http://https://github.com/WMD-group/VMOF. A set of raw phonon data and program input/output files are available from http://https://researchdata.ands.org.au.

Acknowledgements

J. K. B. is funded by the EPSRC (Grant no. EP/G03768X/1). J. D. G. thanks the Australian Research Council for funding under the Discovery Programme, as well as the Pawsey Supercomputing Centre and NCI for the provision of computing resources. A. W. acknowledges support from the Royal Society University Research Fellowship scheme. K. L. S. is funded under ERC Starting Grant 277757 and J. M. S. is funded under EPSRC Grant no. EP/K004956/1. The work benefited from the high performance computing facility, Balena, at the University of Bath, and access to the ARCHER supercomputer through membership of the HPC Materials Chemistry Consortium (EP/L000202).

References

1. B. R. Brooks, C. L. Brooks, A. D. MacKerell, L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels and S. Boresch, J. Comput. Chem., 2009, 30, 1545–1614.
2. B. R. Brooks, R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan and M. Karplus, J. Comput. Chem., 1983, 4, 187–217.
3. N. L. Allinger, Y. H. Yuh and J. H. Lii, J. Am. Chem. Soc., 1989, 111, 8551–8566.
4. N. L. Allinger, F. Li and L. Yan, J. Comput. Chem., 1990, 11, 848–867.
5. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
6. A. K. Rappé, C. J. Casewit, K. Colwell, W. Goddard Iii and W. Skiff, J. Am. Chem. Soc., 1992, 114, 10024–10035.
7. M. A. Addicoat, N. Vankova, F. I. Akter and T. Heine, J. Chem. Theory Comput., 2014, 10, 880–891.
8. J. Borycz, D. Tiana, E. Haldoupis, J. C. Sung, J. I. Siepmann, O. K. Farha and L. Gagliardi, J. Phys. Chem. C, 2016, 120, 12819–12830.
9. S. Bureekaew, S. Amirjalayer, M. Tafipolsky, C. Spickermann, T. K. Roy and R. Schmid, Phys. Status Solidi B, 2013, 250, 1128–1141.
10. L. Vanduyfhuys, S. Vandenbrande, T. Verstraelen, R. Schmid, M. Waroquier and V. Van Speybroeck, J. Comput. Chem., 2015, 36, 1015–1027.
11. J. K. Bristow, D. Tiana and A. Walsh, J. Chem. Theory Comput., 2014, 10, 4644–4652.
12. C. E. Wilmer, M. Leaf, C. Y. Lee, O. K. Farha, B. G. Hauser, J. T. Hupp and R. Q. Snurr, Nat. Chem., 2012, 4, 83–89.
13. A. O. Yazaydin, R. Q. Snurr, T.-H. Park, K. Koh, J. Liu, M. D. LeVan, A. I. Benin, P. Jakubczak, M. Lanuza and D. B. Galloway, J. Am. Chem. Soc., 2009, 131, 18198–18199.
14. M. Haranczyk and R. Martin, J. Phys.: Conf. Ser., 2015, 012103.
15. M. V. Parkes, C. L. Staiger, J. J. Perry IV, M. D. Allendorf and J. A. Greathouse, Phys. Chem. Chem. Phys., 2013, 15, 9093–9106.
16. B. Huang, A. McGaughey and M. Kaviany, Int. J. Heat Mass Transfer, 2007, 50, 393–404.
17. B. Huang, Z. Ni, A. Millward, A. McGaughey, C. Uher, M. Kaviany and O. Yaghi, Int. J. Heat Mass Transfer, 2007, 50, 405–411.
18. H. Babaei and C. E. Wilmer, Phys. Rev. Lett., 2016, 116, 025902.
19. M. Arnold, P. Kortunov, D. J. Jones, Y. Nedellec, J. Kärger and J. Caro, Eur. J. Inorg. Chem., 2007, 60–64.
20. F. Jeremias, V. Lozan, S. K. Henninger and C. Janiak, Dalton Trans., 2013, 42, 15967–15973.
21. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 169.
22. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
23. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
24. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
25. J. D. Gale, J. Chem. Soc., Faraday Trans., 1997, 93, 629–637.
26. J. D. Gale and A. L. Rohl, Mol. Simul., 2003, 29, 291–341.
27. J. Gasteiger and M. Marsili, Tetrahedron, 1980, 36, 3219–3228.
28. J. Gasteiger and M. Marsili, Tetrahedron Lett., 1978, 19, 3181–3184.
29. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
30. A. Togo, L. Chaput, I. Tanaka and G. Hug, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 81, 174301.
31. F. I. Mopsik, J. Res. Natl. Inst. Stan. A-Phys. Chem., 1973, 407–409.
32. N. Allan, T. Barron and J. Bruno, J. Chem. Phys., 1996, 105, 8300–8303.
33. J. D. Gale, J. Phys. Chem. B, 1998, 102, 5423–5431.
34. R. Ahuja, L. Fast, O. Eriksson, J. Wills and B. Johansson, J. Appl. Phys., 1998, 83, 8065–8067.
35. G. Fadda, L. Colombo and G. Zanzotto, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 214102.
36. Z. Jian-Zhi, W. Guang-Tao and L. Yong-Cheng, Chin. Phys. Lett., 2008, 25, 4356.
37. J. Gladden, J. H. So, J. Maynard, P. Saxe and Y. Le Page, Appl. Phys. Lett., 2004, 85, 392–394.
38. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6170.
39. M. Valiev, E. J. Bylaska, N. Govind, K. Kowalski, T. P. Straatsma, H. J. Van Dam, D. Wang, J. Nieplocha, E. Apra and T. L. Windus, Comput. Phys. Commun., 2010, 181, 1477–1489.
40. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
41. A. Van Yperen-De Deyne, K. Hendrickx, L. Vanduyfhuys, G. Sastre, P. Van Der Voort, V. Van Speybroeck and K. Hemelsoet, Theor. Chem. Acc., 2016, 135, 1–14.
42. P. G. Yot, K. Yang, F. Ragon, V. Dmitriev, T. Devic, P. Horcajada, C. Serre and G. Maurin, Dalton Trans., 2016, 45, 4283–4288.
43. J. K. Bristow, K. L. Svane, D. Tiana, J. M. Skelton, J. D. Gale and A. Walsh, J. Phys. Chem. C, 2016, 120, 9276–9281.
44. D. E. Woon and T. H. Dunning Jr, J. Chem. Phys., 1993, 98, 1358–1371.
45. W. Zhou, H. Wu, T. Yildirim, J. R. Simpson and A. H. Walker, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 054114.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6cp05106e

---