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

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$\mathrm{\ddot{u}}$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.


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 parametrising 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 parametrisation 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 nonbonding 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][13][14][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 pardepenameters.
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 finitetemperature properties. 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 deter-mined 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 additonal temperature dependent structural properties.

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) 16 code was used to perform Kohn-Sham density functional theory (DFT) calculations using the PBEsol exchange-correlation functional. 17 The projector-augmented wave method 18 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 sets 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 19 was included and found to be necessary to remove phonon instabilities for some MOF structures. Reference calculations for the binary oxides: ZnO, Al 2 O 3 , TiO 2 and ZrO 2 , 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 SI along with the optimised unit cell parameters.

Forcefield calculations
VMOF was derived using the General Utility Lattice Program (GULP) code, which has extensive capabilities suited to both inorganic and organic materials. 20,21 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 (Equation 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 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; ε i j = √ ε ii ε j j and d i j = d ii d j j , to reduce the number of parameters that require fitting, thus increasing the transferability. The long-range cutoff 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 ε 0 is the vacuum permittivity. Note that the harmonic bonding terms in GULP possess a multiplication factor of 1 2 , and so CHARMM force constants were appropriately scaled.
The charges of the ligands are derived within GULP using the charge equilibration model of Gasteiger 22,23 , 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. 23 Initial charge parametrisation 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 (TPDC) and 2,7pyrene-dicarboxylate (PDC). Once derived for a specific atom type, the same charges are used for all the structures modelled.

Property calculations
2.3.1 Mechanical properties. Bulk moduli (B 0 ) 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.

Vibrational properties.
IR frequencies and intensities, as well as densities of states, were post-processed from both forcefield and first principles methods to apply a consistent broadening factor of 10 cm −1 with a frequency sampling resolution of 1 cm −1 . 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.

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. 24,25 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 SI.
The quasi-harmonic approximation 26 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 27 (Zero Static Internal Stress Approximation) approximation. While this can be routinely achieved using some forcefield implementations 28 , the requirement to compute third-order derivatives makes this particular demanding with quantum mechanical approaches.
Throughout the calculation of thermodynamic properties a consistent Brillioun-zone phonon sampling (q-point) mesh of 32 × 32 × 32, 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.

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 SI) 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.
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.

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 2− ligand in the gas phase. The PBE0 functional 33 was used in the NWChem program 34 with the Dunning correlation consistent cc-pVTZ basis sets 35 to fully relax the ligand with a 0 and 90 • torsion of the carboxylate heads in relation to the aromatic ring (depicted in the SI). 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.
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 ( Figure 6). 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. 36 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 SI).
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/primitive unit cell more stable, with no ex-  Table 2. ternal 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.
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.

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).
The forcefield predicts smaller bulk moduli than DFT. 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. 37 The trend in mechanical strength for each framework is the same between the two methods, with increasing in 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 Znisorecticular MOFs, which possess large internal voids and weaker van der Waals interactions between Zn 2+ and organic ligands.

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 ( Figure 2) and the forcefield (Figure 3). Comparison of the plots highlights 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 of 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.
The decomposition of the DOS into elemental contributions is shown in the partial DOS plots (Figures 4 and 5). 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 parametrising 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.

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. 38 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 (Figure 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 c d reproduce the vibrational properties of the subset of MOFs studied.

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 sys-tems 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. 34,35,39 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 SI). Rotating the ligand in this manner in MOF-5 is shown to cost little energy with both DFT and VMOF up to 10 • ( Figure  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 method with first principles calculations, we report a comparison of properties for DFT and VMOF for only four structures; MOF-5, IRMOF-10, UiO-66 and NOTT-300 (Figure 9), while quasi-harmonic calculations are also performed for MIL-125, UiO-67, MOF-650 and MOF-74 with the forcefield (Figure 10).
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.
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 highpressure 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 (Figures 9 and 10). Firstly, a comparison of thermal expansion coefficients show the vast difference in structural changes with temperature. The "softer" Zn-isorecticular 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 (Figure 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 constant pressure heat capacity at constant pressure are reproduced accurately and appear unaffected by the discrepency 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.

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 quasiharmonic 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.

Acknowledgments
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).

Data Access Statement
The VMOF forcefield is available as a library file for GULP from https://github.com/WMD-group/VMOF. A set of raw phonon data and program input/output files are available from https://researchdata.ands.org.au. 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 i j ), bulk moduli (B 0 ), unit cell parameters and metal-oxygen (M-O) bond lengths. Percentage errors in the cell parameters are given in parenthesis.  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.