Transferable Force Field for Metal-Organic Frameworks from First-Principles: BTW-FF.

We present an ab-initio derived force field to describe the structural and mechanical properties of metal-organic frameworks (or coordination polymers). The aim is a transferable interatomic potential that can be applied to MOFs regardless of metal or ligand identity. The initial parametrization set includes MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1. The force field describes the periodic crystal and considers effective atomic charges based on topological analysis of the Bloch states of the extended materials. Transferable potentials were developed for the four organic ligands comprising the test set and for the associated Cu, Zn, and Zr metal nodes. The predicted materials properties, including bulk moduli and vibrational frequencies, are in agreement with explicit density functional theory calculations. The modal heat capacity and lattice thermal expansion are also predicted.

Introduction

Metal–organic frameworks (MOFs), or coordination polymers, are hybrid materials that display a range of unique properties that can be used for a variety of applications ranging from gas storage and catalysis to photovoltaics and drug delivery systems.[1−3] MOFs are formed via the coordination of metals and organic ligands in a self-assembled manner to create an extended crystalline material. The number of possible combinations of building blocks is therefore vast, and the associated crystal structures and properties are often unpredictable. To reduce efforts spent on refining the optimum conditions for the synthesis of a particular MOF, and to identify compositions of particular interest, accurate property predictions from computer simulations would be highly beneficial.

The emerging field of “materials design” has largely been based around the application of modern electronic structure techniques, such as density functional theory (DFT), to predict the structures and properties of new materials.[4−8] Such approaches are appropriate for high-symmetry close-packed inorganic materials but are challenging for the large and complex crystal structures associated with MOFs. For example, a “complex” quaternary semiconductor such as Cu2ZnSnS4 can be described using only 8 atoms in its primitive unit cell,[9] while a “simple” MOF can require several hundred atoms. For example, the popular framework MOF-5 has 106 atoms in its primitive cell and 424 atoms in its conventional crystallographic cell. A high-quality calculation of a single MOF poses a heavy computational burden; a large-scale screening is prohibitively expensive.

An alternative to a direct quantum mechanical treatment, which usually involves a self-consistent numerical solution of the Kohn–Sham or Hartree–Fock equations, is the use of an analytical interatomic potential or force field that can describe a range of properties of interest.[10,11] An accurate and transferable force field for MOFs would provide a means of rapidly screening material compositions and properties for particular applications.

A number of initial MOF force fields have been recently reported.[12−16] These have been mainly fitted for specific MOF structures that were highlighted experimentally as possessing functional properties. An example by Vaduyfhuys et al. is a force field for the Al containing MIL-53.[17] One common approximation for MOF potentials, to reduce the complexity of parametrization, is to fix the atom positions within the unit cell. This approach is generally used for probing gas adsorption using Grand Canonical Monte Carlo (GCMC) or Molecular Dynamic (MD) simulations. The advantage is that the intraframework interactions need not be considered, facilitating fast screening,[2,18−23] but materials response functions are not described. As an extension to this approach one can conduct hybrid GCMC and MD calculations to model structural transitions of flexible MOFs.[24−26]

A second approach is to use a generic force field derived for molecules such as proteins, hydrocarbons, and common gases; OPLS-aa (Optimized Potentials for Liquid Simulations—all atom), DREIDING, and MM3 (Molecular Mechanics 3) are examples of these.[27−29] The application of generic force fields to MOFs, which consist of both organic and inorganic components, is highly convenient but approximate. While large-scale screenings can be again performed,[30,31] complex geometries and interactions are often not modeled with quantitative accuracy. Recent progress includes UFF4MOF, an extension of the universal force field (UFF)[30] to describe some common metal–organic framework motifs.[32]

A third common approach to MOF force fields is to construct representative finite clusters (fragments) of the full MOF crystals.[33] The advantage is faster derivation with a reduced number of interaction parameters (and degrees of freedom) in the model. The disadvantage is that the neglect of periodicity and long-range interactions is unphysical and standard mechanical, thermal, and dielectric properties cannot be described in the absence of a sophisticated embedding procedure.

In this paper, we report a force field to describe metal–organic frameworks parametrized by first-principles crystal structures and electron density. The aim is a transferable potential form suitable to describe the majority of ligand and metal combinations for MOFs, including predicting properties of novel compositions. In contrast to the universal force field approach in which general parameters not fitted for MOFs and fixed generic charges are employed, we develop a force field that has been fitted explicitly to the periodic frameworks. Some initial parameters have been refined from existing force fields (MM3 and MOF-FF) and the functional format of MM3 is preserved.[34−36] Ligand and metal interatomic potentials are parametrized and validated across multiple MOF structures and atomic charges are determined based on topological analysis using Bader’s atom in molecules theory (AIM)[37] of the equilibrium electron density from solid-state DFT calculations.

Metal–Organic Frameworks

We consider six representative MOFs covering three of the most common isoreticular frameworks: [MOF-5 (Zn4O(BDC)3), IRMOF-10 (Zn4O(BPDC)3), IRMOF-14 (Zn4O(PDC)3), UiO-66 (Zr6O4(OH)4(BDC)6), UiO-67 (Zr6O4(OH)4(BPDC)6), and HKUST-1 (Cu3(BTC)2)], as illustrated in Figures 1 and 2.

Figure 1

Comparative illustrations of the repeating units of Zn-containing MOF-5, IRMOF-10 and IRMOF-14 with the organic ligands shown underneath each structure. Torsion angle (atom types 170-913-902-912) labeled as x has been highlighted.

Figure 2

Comparative illustrations of the repeating units of Cu-containing HKUST-1 and Zr-containing UiO-66 and UiO-67 with the organic ligands shown underneath each structure. Torsion angle (atom types 912-903-903-912) labeled as y has been highlighted.

MOF-5 was first reported by the group of Yaghi.[38] It is formed of Zn4O clusters and 1,4-benzenedicarboxylate (BDC) organic ligands. The Zn metal centers are in tetrahedral coordination with respect to oxygen (Figure 1), which is similar to the bulk metal oxide. The crystal structure can be described by its primitive rhombohedral unit cell of dimensions a = 18.289 Å, α = 60.0°. The cubic unit cell of MOF-5 has the associated dimension of a = 25.832 Å and space group Fm3̅m.[39] The ligands form the edges of the cubic structure and the metal clusters the corners. Each unit cell contains half the ligand molecules orientated face-on and half rotated 90° into the plane. The structures of IRMOF-10 and IRMOF-14 are analogous to that of MOF-5, differing only in the ligand identity. The organic ligand comprising IRMOF-10 is 4,4′-biphenyldicarboxylate (BPDC), whereas for IRMOF-14 it is pyrene-2,7-dicarboxylate (PDC). The cell parameter for the cubic unit cell of IRMOF-10 is a = 34.281 Å, while IRMOF-14 is a = 34.381 Å.[39]

Beyond the isoreticular Zn-MOFs, UiO-66, UiO-67, and HKUST-1 were also modeled (see Figure 2). UiO-66 is formed of Zr6O4(OH)4 clusters 12-coordinated to BDC ligands, such that each Zr ion is in octahedral coordination with capping η3-OH and η3-O anions.[40,42,43] The cubic cell dimension for UiO-66, with the space group F4̅3m, is a = 20.978 Å.[40] The crystal structure of UiO-67 is similar but composed of the biphenyl ligand (BPDC); the extension is analogous to that of MOF-5 to IRMOF-10. The behavior of the BPDC ligand differs for UiO-67 due to a ligand twist (labeled y in Figure 2). First-principles calculations predict the twist across the central carbon atoms connecting the two aromatic rings to be approximately 31°. UiO-67 also has the space group F4̅3m, with cell dimension a = 27.094 Å.[40,44] The UiO series is of increasing interest for MOF applications due to their high thermal stability up to 813 K and resistivity to water decomposition.[45]

HKUST-1 differs in relation to the structures previously discussed. Here, the ligand is 1,3,5-benzenetricarboxylic acid (BTC), as shown in Figure 2. This MOF, first reported by Chui et al.,[41] is of interest not only for its high gas storage capacities but also for its unique electronic structure originating from Cu–Cu interactions. The room temperature experimental crystal structure infers metal–metal separation of 2.63 Å. The metal nodes in HKUST-1 are Cu, which are in a square planar coordination resulting in a 3D network with a paddlewheel geometry and three voids of diameters 5, 11, and 13.5 Å, respectively. Each metal center has 4 BTC ligands, with additional water coordination in vertical alignment with the Cu–Cu interaction, resulting in a pseudo-octahedral geometry for hydrated crystals. The water can be removed and/or substituted to further expose the metals.[46,47] In this paper, we will consider only the dehydrated HKUST-1 structure. HKUST-1 has a cubic crystal symmetry (space group Fm3̅m) with a lattice dimension of a = 26.343 Å. Table 1 gives a summary of unit cell parameters of each structure described.

Table 1

Experimentally Determined Crystal Structure Parameters of Six Metal–Organic Frameworksa

MOF	metal	oxidation	ligand	space group	a (Å)	N
MOF-5[39]	Zn	II	BDC	Fm3̅m	25.832	424
IRMOF-10[39]	Zn	II	BPDC	Fm3̅m	34.281	664
IRMOF-14[39]	Zn	II	PDC	Fm3̅m	34.381	760
UiO-66[40]	Zr	IV	BDC	F4̅3m	20.978	456
UiO-67[40]	Zr	IV	BPDC	F4̅3m	27.094	684
HKUST-1[41]	Cu	II	BTC	Fm3̅m	26.343	576

Oxidation refers to the formal metal oxidation state, while N refers to the number of atoms in the unit cell described.

Theoretical Approach

Reference Solid-State Electronic Structure Calculations

In order to provide systematic benchmark data, electronic structure calculations for the periodic crystals were performed using DFT as implemented in the VASP (Vienna Ab-Initio Simulation Package) code.[48,49] All calculations were performed using the PBEsol functional,[50] which is a semilocal functional that usually predicts equilibrium structures in very good agreement with experiment; its success for MOFs has been well demonstrated.[3,51] Comparison of different exchange-correlation functionals for a range of materials shows that the structures and frequencies of PBEsol are among the most accurate currently available for solids.[52,53]

The computational setup differed between structures due to variations in the unit cell size. A 2 × 2 × 2 k-point grid was used for the optimization of UiO-66. Due to the larger unit cells of MOF-5, IRMOF-10, IRMOF-14, UiO-67, and HKUST-1, only Γ-point sampling was performed. The quasi-Newtonian relaxation employed for structural optimization was converged to forces of 0.005 eV/Å or lower. A kinetic energy cutoff of 500 eV was employed for the plane-wave basis set. Scalar-relativistic projector-augmented wave (PAW) potentials were used to model the interactions between the core and valence electrons on each atom, with the 3d electrons treated explicitly for Zn. Effective atomic charges for each atom type were derived using the AIM theory on the total electron density (i.e., the sum of the PAW and valence density) for the optimized structure.[54] Vibrational frequencies were calculated with Γ-point sampling of the Brillouin zone using the finite displacement method.

Force Field Parametrization and Testing

The functional form of MM3 is maintained in our parametrization. Thus, the overall energy expression iswhere str = stretch, tor = torsion, opb = out of plane bend, Coul = Coulomb, and vdW = van der Waals interactions and where the usual bond stretching and bending modes are described, in addition to longer range van der Waals (dispersion) and Coulombic interactions. Nonbonding interactions were calculated using the Buckingham equation:vdW radii and ε values are given in Table 4. Default values for the A, B, and C constants were used (184000.0, 12.0, 2.25, respectively) as defined in TINKER.

Table 4

van der Waals Radii and ε Values Used for the Transition Metals within the MOF Structuresa

elements	vdW radii (Å)	ε (kcal mol–1)
Zn	2.290	0.276
Cu	2.290	0.276
Zr	3.520	0.367

Epsilon refers to the polarizability of the atoms, which is an energy term within the van der Waals function in the MM3 format (eq 1).

The MM3 force field has been shown to recreate organic systems accurately and, more recently, has been applied to MOF structures.[33−36,55] The TINKER package[56,57] contains the full set of the MM3 force field parameters and has a range of capabilities for crystal structure and property calculations. The XTALMIN program inside of TINKER was used to apply periodic boundary conditions allowing the full crystal to be modeled and optimized. The internal positions are described using connectivity (in the absence of space group symmetry). When considering the interaction between organic and inorganic building blocks both “bonding” and “nonbonding” contributions were included. Atom types were assigned based on the element and the environment of the crystal. The reparameterization of the MM3 force field include the terms describing the carboxylic head and interaction between metal node and ligand. New parameters were also derived for the metal node, particularly for the metal–inorganic oxygen interaction. TINKER default values of the Buckingham potentials were used in association with Coulomb long-range terms to describe nonbonding interactions.

Vibrational frequencies were calculated at the Γ-point using the VIBRATE program, and bulk moduli with the DYNAMIC program in the TINKER package through molecular dynamic simulations. A series of Canonical (NVT) ensemble calculations were employed to fix the unit cell volume at 1 K with negligible kinetic energy contribution. Unit cell volumes were sampled every 0.01 Å for ±1% from the equilibrium volume (V0). Simulations were run for 250 000 dynamic steps with a time step of 1 fs at 1 atm with an Ewald cutoff of 11 Å. Constant temperature of the system was maintained with a Nosé–Hoover thermostat. The Velocity Verlet algorithm was used to integrate the Newton equations. Average values of the potential energy were taken from the final 50 ps (500 dynamic steps). The bulk moduli (B0) of each structure was calculated using the isothermal Birch–Murnaghan equation of state (eq 3). Data processing was implemented in Octave using the Asturfit program, which performs a least-squares fit to the Birch–Murnaghan equation of state.[58,59]

Linear (α) (eq 4) and volumetric (β) (eq 5) thermal expansion coefficients were calculated from a series of anisotropic isothermal–isobaric (NPT) ensemble calculations. The Berendensen bath coupling method was used as a barostat to maintain defined pressures. The temperatures were sampled at 1 K and between 80–500 at 50 K increments with 500 K as an additional temperature point. The program DYNAMIC was used to calculate the average lattice parameters over 50 ps following equilibration for 200 ps.

For isotropic expansion 3 × α = β; this relationship was used when calculating the volumetric thermal coefficients from the linear thermal coefficients. The value of the lattice parameter a when calculating the linear thermal expansion coefficient was taken as the average of the 298 K MD simulation and ∂a/∂T is the average slope over the temperature range 80–500 K.

Finally, isochoric heat capacities (CV) were calculated within the harmonic approximation from the vibrational frequencies of each structure for temperatures ranging between 80 and 500 at 50 K intervals using the standard (Einstein) phonon summation:

Molecular versus Periodic Charges

One distinction of our MOF potential is the choice of effective atomic charges. When deriving effective charges the core atomic density is often not considered due to the use of a pseudopotential or effective core potential. This leads to a systematic error in the calculation of the atomic surfaces, making bonds appear more ionic in nature.[60] With the PAW method, we can reconstruct the total charge density of the system as a sum of the valence and PAW density. The effective charges derived from topological analysis of charge density for different structural representations of MOF-5, including those of BTW-FF (with and without the core density) are listed in Table 2.

Table 2

Comparative Bader Charges Derived for a Small Cluster (SC) (79 Atom, Zn4O13C42H30) and Large Cluster (LC) (328 Atom, Zn32O104C120H72) of MOF-5 Compared with Those Derived for the Periodic System Used in the BTW-FF and UFF Modelsa

		Bader charges (au)
atom type	element	BTW-FFinc.core	BTW-FFexc.core	SC	LC	UFF
172	Zn	1.281	1.408	1.291	1.291	1.308
913	C (acid)	1.497	2.683	1.558	1.536	1.912
912	C (benz)	–0.053	–0.011	–0.058	–0.055	1.912
902	C (C–Cacid)	–0.008	–0.007	0.010	–0.003	1.912
170	O (acid)	–1.151	–1.768	–1.195	–1.168	–2.300
171	O (inorganic)	–1.115	–1.336	–1.171	–1.207	–2.300
915	H (H–C)	0.126	0.083	0.090	0.123	0.712

Also given are the periodic charges with and without the inclusion of core density. Charges are given in atomic units.

Analysis of data first confirms the importance of including core charges for MOFs. The differences are large, especially for the carboxylic head of the ligand. Second, comparing charges derived from the Bloch states versus finite molecular orbitals, we show that similar results can be obtained when a large cluster is used. In contrast, the standard molecular UFF parametrization has a more ionic description with charges larger in magnitude for all atoms.

Results

To begin, topological analysis of charge density of the equilibrium MOF structures from the electronic structure calculations was performed to obtain the effective atomic charges to be used in the force field. The values are presented in Table 3. Comparable charges were derived for MOF-5, IRMOF-10, and IRMOF-14, and also for UiO-66 and UiO-67. The higher charge of Zr is consistent with its higher oxidation state. Similar derived charges for each atom type allowed average charges to be used, resulting in a fully transferable model, analogous to that of an UFF implementation.

Table 3

Atomic Charges Derived Via Topological Analysis for Each Atom Type in Each MOFa

		effective atomic charges (au)
element	atom type	MOF-5	IRMOF-10	IRMOF-14	UiO-66	UiO-67	HKUST-1	avg. charges
C (benz)	912	–0.054	–0.046	–0.044	–0.058	–0.060	–0.023	–0.050
C (acid)	913	1.497	1.538	1.538	1.576	1.548	1.540	1.539
C (Cbenz–Cacid)	902	0.008	–0.028	0.051	–0.056	0.006	–0.011	–0.008
C (Cbenz–Cbenz′)	903		–0.012			–0.035		–0.024
O (acid)	170	–1.151	–1.163	–1.156	–1.181	–1.182	–1.091	–1.154
O (inorganic)	171	–1.115	–1.214	–1.224	–1.189	–1.881		–1.186
O (O–H)	75				–1.242	–1.244		–1.243
H (H–C)	915	0.126	0.105	0.092	0.129	0.096	0.158	0.118
H (H–O)	21				0.622	0.623		0.622
metal		1.281	1.295	1.297	2.601	2.610	1.036	1.291 (Zn), 2.605 (Zr)

Atomic charges are given in atomic units and total average charges of all structures are given.

The values used for the van der Waals radii and polarizability (ε) parameters describing the nonbonding interactions in each structure are given in Table 4. The derived parameters of Zn and Cu are similar to those derived by Schmid et al.[34] Zr is heavier and its higher oxidation state in UiO-66 results in greater polarizability; hence, its values are larger. The Zr parameters are similar to those already used within inorganic crystals.[61]

The final set of force field parameters to describe all MOFs reported here are included as Supporting Information (SI) and an online repository, which will be updated with future potentials.[62]

IRMOF Series

The lattice parameter agreement between the equilibrium BTW-FF and DFT (PBEsol) structures is very good (errors of less than 1%). The full internal structural parameter comparison for MOF-5, IRMOF-10, and IRMOF-14 is provided as SI. Atom type assignment is shown in Figure 3.

Figure 3

Atom type definitions used for MOF-5/IRMOF-10/IRMOF-14.

An interesting observation was made when deriving the potentials for both MOF-5 and IRMOF-10. Initial parametrization resulted in slightly distorted ligand structures with a nonplanar torsion angle occurring between the O-C(carb)-C(phen)-C(phen) (170-3-902-2) (labeled as x in Figure 1). For comparison a DFT(PBEsol) calculation was run starting from a 2.73° angle, where the equilibrium structure gave a distortion of 1.04°. Our data suggests that the ground-state of the ligands in MOF-5, IRMOF-10, and IRMOF-14 is slightly distorted, with a space group of lower symmetry than the average one identified at room temperature. This distortion is likely to fluctuate at higher temperatures, with an average structure that is planar (double well potential). The force field parameters presented contain a torsion potential to suppress this distortion and recreate experimental data; however, it can be removed without disrupting the remainder of the framework.

UiO-Series

Compared to the IRMOF series, the structural fit of UiO-66 and UiO-67 proved more challenging to converge. The equilibrium lattice constants are within 0.51% of the DFT(PBEsol) values. The ligand potentials used to model the IRMOF structures are transferable to UiO-66 and UiO-67 with small errors in bond lengths and angles (see SI). Atom type assignments are shown in Figure 4. The results and the error associated with the Zr–inorganic oxygen bond (2.23 and 2.67% in UiO-66 and UiO-67, respectively) suggest a more robust nonbonding interaction for these bonds may be required to improve the overall accuracy of the models (e.g., higher-order multipoles in the electrostatic summations).

Figure 4

Definition of atom types for UiO-66 and UiO-67.

To stress the close relationship between MOF-5 and UiO-66, the ligand potentials for the BDC ligand remained unchanged between the structures. This is also the case for IRMOF-10 and UiO-67. The one exception is that in order to recreate the 31° twist across the central bond of the biphenyl ligand within UiO-67, the applied k-force for the torsion (912-903-903-912 atom types) was reduced. The resultant twist of UiO-67 was 32.2°. The Zr potential is fully transferable between UiO-66 and UiO-67.

HKUST-1

The equilibrium lattice constants of HKUST-1 show errors of less than 0.1%, with the internal structure parameters also well described (see SI). See Figure 5 for atom type assignment. The nature of the Cu–Cu bonding is somewhat in debate; chemically each Cu ion is formally divalent (d9) and the unpaired spin form open-shell and closed-shell singlets, as well as a (ferromagnetic) triplet state.[63,64] Here, it was necessary to model Cu–Cu as 5 coordinate, that is, bonded to 4 carboxylic acid oxygens in the equatorial positions and 1 Cu in the axial position. Reasonable k-force was required to model this metal–metal interaction. Future work could extend these HKUST-1 potentials to include water coordination as described in the introduction of this paper to form a pseudo Cu octahedral environment.

Figure 5

Definition of atom types for HKUST-1.

Property Calculations

To validate the accuracy and transferability of the potential model beyond equilibrium structures, a series of materials property calculations have been performed. We have determined the bulk moduli and vibrational frequencies for each structure using our force field and compared these to available reference data. These properties are extremely computationally expensive to calculate on the DFT potential energy landscape but are straightforward using our potential model.

The calculated bulk moduli are given in Table 5, with the energy/volume curves provided in Figure 6. An excellent agreement with reference values for all structures besides UiO-66 is shown. The trends in the obtained bulk moduli of the MOFs is also consistent with available experimental and DFT calculated values. UiO-66 appears to be an unique case due to the large increase in mechanical strength for UiO-66 when comparing this value with the structurally similar UiO-67. To our knowledge, experimental values for the bulk modulus of UiO-66 are not available currently to provide a definitive reference. The large difference between the bulk moduli of MOF-5 and UiO-66 is due to the coordination and oxidation state differences of the metal centers in either structure. In MOF-5 the metals are 4-coordinate and formally Zn2+ and in UiO-66 the metals are 6-coordinate with 6 further capping ligands to create a total outer coordination of 12. In addition, the formal oxidation state of the metals in UiO-66 is Zr4+. UiO-66 can therefore be considered as being formed of stronger interactions, which increase the mechanical strength of the material.

Table 5

Bulk Moduli of MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1a

MOF	Bulk moduli (GPa)
	BTW-FF	Reference
MOF-5	11.95	18.20[65]
IRMOF-10	8.25	6.00[66]
IRMOF-14	8.40	5.90[66]
UiO-66	27.15	41.01[65]
UiO-67	19.15	17.15[65]
HKUST-1	25.05	24.53[65]

Values reported are those using BTW-FF and available reference data from DFT calculations. Reference calculations used Density Functional based Tight Binding (Kuc et al.)[66] and PBE functional (Wu et al.).[65] Note that the bulk modulus is related to the second derivatives of the energy with respect to volume and hence is sensitive to the theoretical approach.

Figure 6

Energy/volume curves for MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67 and HKUST-1, from which the bulk moduli were derived via an equation of state.

It should be noted that the bulk modulus of a material is temperature dependent. Direct comparison with experimental data is therefore not often appropriate when discussing accuracy to athermal calculations.[67] Furthermore, where values from other DFT studies are compared, it should be remembered that these values are sensitive to the choice of functional, where the bulk modulus is proportional to the equilibrium cell volume, that iswhere standard local and semilocal functionals can produce large errors.

Vibrational frequencies were calculated and compared with those from DFT (PBEsol) (see Figure 7). A generally good agreement is found for all structures. There are, however, frequencies in the 1900 cm–1 to 3000 cm–1 range for UiO-67 and IRMOF-10, which are not calculated with DFT. We can assign these anomalous modes to stretches involving the central C–C bond (atom types 903-903). The force constants involving the connection of this central bond were strengthened to prevent the previously discussed carboxylate head twisting distortion occurring. Further additional modes for UiO-67 and also for UiO-66 were predicted by BTW-FF in the same frequency region. Analysis of the associated eigenvectors confirms these to be due to stretching modes within the Zr–O(O–H) (atom types Zr-75) bonds. These bonds have the largest error (2.23% and 2.67% for UiO-66 and UiO-67, respectively, Supporting Information Table S2) when compared to the equilibrated DFT structure. HKUST-1 was not included in the vibrational frequency analysis as the DFT proved to be too expensive to compute due to the larger crystal structure. It should be noted that soft modes associated with the carboxylate torsion of the ligands were present in all cases. This torsion has previously been described and labeled as x in Figure 1.

Figure 7

Γ-point vibrational frequencies between 500–3500 cm–1 for MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1. DFT calculated (left) and BTW-FF calculated (right). Note: DFT (PBEsol) frequencies for HKUST-1 could not be computed owing to the computational expense.

As a further analysis, we compare the vibrational frequencies of the IRMOF and UiO-series with DFT (B3-LYP) data previously reported and analyzed by Civalleri et al. and Valenzano et al. for the respective series of structures (Table 6 and Table 7).[44,68] Assignment of the vibrational frequencies for prominent stretching and bending modes shows an excellent agreement in our reported results using BTW-FF and those from DFT (B3-LYP). This agreement further supports the accuracy of our model.

Table 6

Assigned Vibrational Frequencies (ω) for MOF-5, IRMOF-10, and IRMOF-14a

	ω (cm–1)
description	DFT (MOF-5)	MOF-5	IRMOF-10	IRMOF-14
Oacid–Zn–Oacid bend	114	111	123	136
Zn–Oinorganic–Zn bend	136	174	181	178–179
Zn–Oacid–Cacid bend	263	283	302–305	284
Zn–Oinorganic asymmetric stretch	512	497	497–498	493
Zn–Oacid symmetric stretch	579	558	544–575	545–565
Zn–Oacid asymmetric stretch	606	563–568	590–602	575–590
Cacid–Oacid symmetric stretch	1421	1394	1471–1473	1355–1382

Reference DFT calculations B3-LYP level of theory in the CRYSTAL code (Civalleri et al.). Reported DFT values are for MOF-5.[68]

Table 7

Assigned Vibrational Frequencies (ω) for UiO-66 and UiO-67a

	ω (cm–1)
description	DFT (UiO-66)	UiO-66	UiO-67
Zr–Oacid asymmetric stretch	556	558	593
μ3–O stretch	673	671–672	667–671
O–H bend	771, 814	778–779, 810–812	872
Cacid–Oacid symmetric stretch	1408	1380, 1408	1363

Reference DFT calculations B3-LYP level of theory in the CRYSTAL code (Valenzano et al). Reported DFT values are for UiO-66.[44]

Occupation of the lattice phonons at finite temperatures leads to changes in the crystal structure parameters. Linear and volumetric thermal expansion coefficients have been calculated to determine the change in unit cell size with increasing temperature. Many MOFs are known to contract with increasing thermal energy due to transverse vibrational modes of the ligands. Details describing the causes of negative thermal expansion in solid materials are detailed in a referenced review.[69] The phenomenon is not to be confused with structural changes that occur following the evacuation of internal solvent molecules.

Negative thermal expansion was calculated for each structure as shown by the thermal expansion coefficients (Table 8). The extent of unit cell contraction with temperature is shown in Figure 8. IRMOF-10 experiences the greatest thermal contraction due to bending modes of the biphenyl ligand. IRMOF-14 and MOF-5, although containing the same metal node as IRMOF-10, have more rigid ligands and therefore contract less with temperature. The UiO series contract the least due to the high charge on the Zr metal node resulting in a rigid extended structure. The influence of the charge on the metal is highlighted by the difference in behavior of UiO-67 and IRMOF-10. These structures contain the same ligand, but due to the higher ionicity of UiO-67, soft bending modes of the biphenyl ligand are no longer present. The final structure, HKUST-1 also exhibits very weak negative thermal expansion at high temperatures. A low charge was calculated on the Cu metal centers (Table 3) suggesting that HKUST-1 exhibits harder structural properties due to the rigidity of the tricarboxylate ligand and not increased ionicity. Calculated values are consistent with those previously found from experiment and MD simulations.[70−72]

Table 8

Calculated Linear (α) and Volumetric (β) Thermal Expansion Coefficients of MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1

MOF	α (×10–6) (K–1)	β (×10–6) (K–1)
MOF-5	–5.27	–15.80
IRMOF-10	–8.11	–24.32
IRMOF-14	–4.95	–14.86
UiO-66	–1.04	–3.11
UiO-67	–2.22	–6.66
HKUST-1	–3.18	–9.53

Figure 8

Thermal expansion profiles for MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1. In contrast to standard quasi-harmonic approaches, molecular dynamic simulations include high order anharmonicity from phonon–phonon interactions. Shown is the difference in a parameter at a given temperature, with respect to the a parameter at 1 K. Error bars indicate the variation in the lattice constants at each temperature due to thermal fluctuations.

Finally, volumetric heat capacities were calculated at constant volume from the vibrational frequencies (eq 6) for each structure under the harmonic approximation. The heat capacity describes the energy required to increase the temperature of a material by a given quantity and is determined from the changes in vibrational occupancy with increasing temperature. The values plateau at relatively low temperatures (200–300 K) suggesting low Debye temperatures (Figure 9). Unfortunately, little experimental data is available on the heat capacity of MOFs to date for comparison.

Figure 9

Temperature dependence of the predicted volumetric heat capacities for MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1.

Conclusions

We have parametrized a new interatomic potential to describe the crystal structures and properties of metal–organic frameworks. The force field is shown to accurately reproduce the structural parameters of MOF-5, IRMOF-10, IRMOF-14, UiO-66, UiO-67, and HKUST-1. The ligand parameters are transferable between the Zr and Zn frameworks and are expected to be valid for other systems of interest. Bulk moduli and vibrational frequencies have been calculated and are in agreement with calculations on the DFT (PBEsol) potential energy surface. Finally, we highlighted the importance of periodic boundaries when deriving empirical parameters for MOFs. Future work will concern extending this model to other systems and extending the range of materials response functions that can be calculated for hybrid frameworks.