On the importance of anharmonicities and nuclear quantum effects in modelling the structural properties and thermal expansion of MOF-5

In this article, we investigate the influence of anharmonicities and nuclear quantum effects (NQEs) in modelling the structural properties and thermal expansion of the empty MOF-5 metal-organic framework. To introduce NQEs in classical molecular dynamics simulations, two different methodologies are considered, comparing the approximate, but computationally cheap, method of generalised Langevin equation thermostatting to the more advanced, computationally demanding path integral molecular dynamics technique. For both methodologies, similar results were obtained for all the properties under investigation. The structural properties of MOF-5, probed by means of radial distribution functions (RDFs), show some distinct differences with respect to a classical description. Besides a broadening of the RDF peaks under the influence of quantum fluctuations, a different temperature dependence is also observed due to a dominant zero-point energy (ZPE) contribution. For the thermal expansion of MOF-5, by contrast, NQEs appear to be only of secondary importance with respect to an adequate modelling of the anharmonicities of the potential energy surface (PES), as demonstrated by the use of two differently parametrised force fields. Despite the small effect in the temperature dependence of the volume of MOF-5, NQEs do however significantly affect the absolute volume of MOF-5, in which the ZPE resulting from the intertwining of NQEs and anharmonicities plays a crucial role. A sufficiently accurate description of the PES is therefore prerequisite when modelling NQEs. © 2019 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5085649


I. INTRODUCTION
Within the past few decades, the field of atomistic computer simulations has established itself as an important and complementary tool to experimental research in the characterisation of various types of materials. This computational approach does not only provide information about the fundamental molecular features underlying the material's properties but also allows for a large scale screening of potentially interesting materials, where new material designs can be assessed without the burdening requirement of synthesis that precedes experimental characterisation. 1 For metal-organic frameworks (MOFs), [2][3][4] in particular, a screening by means of atomistic simulations is recommended as these crystalline materials are built up out of inorganic building blocks that are connected to one another through organic linkers, giving rise to a virtually unlimited number of possible frameworks that can be constructed out of the vast amount of building blocks. [5][6][7] The designer aspect of MOFs along with the versatility in properties of the different MOF structures makes this class of nanoporous materials most attractive for industrial applications, 8 ranging from gas storage and gas separation 9,10 over heterogeneous catalysis 11 to shock absorption for the more flexible framework materials. 12 Computationally, material properties can be obtained from molecular dynamics (MD) simulations as time averaged quantities by integrating the equations of motion throughout time. By relying on the Born-Oppenheimer approximation, as is usually done in MD simulations, one can separate the dynamics of the nuclei from the dynamics of the electronic structure of the system. Depending on the required length of the simulation and the size of the molecular system, the potential energy surface (PES) is then described ab initio, typically by means of density functional theory (DFT), or analytically, using force fields. However, irrespective of the description of the PES, a vast majority of the MD simulations treats the atomic nuclei as classical particles, in spite of their quantum mechanical nature. For sufficiently heavy atoms at sufficiently high temperatures, subjected to moderate interatomic interactions, this approximation is expected to be reasonable. Lighter atoms at lower temperatures or atoms subjected to strong interatomic interactions on the contrary cannot be described as classical particles, but have to be treated quantum mechanically as nuclear quantum effects (NQEs) such as zero-point energy (ZPE) and tunnelling can induce significant deviations from the classical behaviour. For water, a prototypical case study within the context of NQEs, the importance of NQEs has been demonstrated for ample of aspects, 13,14 including the occurrence of transient autoprotolysis events in liquid water 15 and a proper description of the heat capacity of both liquid water and ice. 16 Given the large ZPE of about 21 kJ/mol per O− −H stretch, 13 which is in sharp contrast with the classical contribution to the kinetic energy of only 3.7 kJ/mol at 300 K (as predicted by the equipartition law), the importance of NQEs in water is abundantly clear.
Within the context of MOFs, NQEs are hardly ever taken into account, despite the presence of many light framework atoms, such as hydrogen and carbon. Only when studying water molecules confined in the pores of a framework, NQEs have been included to some extent. 17,18 In this work, the importance of NQEs in the structural and thermal properties of the empty framework of MOF-5 19 is assessed by means of MD simulations. In this way, the NQEs inherent to the MOF itself are investigated rather than the NQEs associated with water present in its pores, using MOF-5 as a prototypical rigid framework (Fig. 1). The structural properties of MOF-5 are probed by means of radial distribution functions (RDFs), which yield information about the interatomic equilibrium distances and the thermal activation of bond stretches of different atom pairs. The thermal property under review is the well-known negative thermal expansion of MOF-5, 20-25 for which a delicate interplay between the modelling of anharmonicities and NQEs is demonstrated.
For each of these properties, the importance of NQEs is assessed by comparing ordinary, classical MD simulations to MD simulations including NQEs, where the quantum mechanical description of the nuclei is pursued using two different methodologies. The first method is the so-called generalised Langevin equation (GLE) quantum thermostat, 26 which allows us to include harmonic NQEs at virtually no additional computational cost. The second method that will be considered is the more advanced technique of path integral MD (PIMD), 27,28 which entails a more rigorous description of both harmonic and anharmonic NQEs, at a sometimes considerable additional computational expense. By combining both methods, the computational cost of PIMD simulations can however be reduced by a factor of about four 29,30 and therefore, this synergistic approach will also be considered, alongside the GLE quantum thermostat and regular PIMD, as explained in Sec. II.

II. METHODOLOGY: MODELLING NUCLEAR QUANTUM EFFECTS
A. GLE quantum thermostat A first, but approximate technique to include NQEs in MD simulations is the so-called generalised Langevin equation (GLE) quantum thermostat. 26 Just as any other thermostat, the GLE thermostat imposes a sampling of the canonical ensemble instead of the microcanonical ensemble by adapting Hamilton's equations of motion, to allow for an exchange of heat with the system's surroundingṡ with q the generalised coordinate describing the position of a particle with a unit mass, p the conjugate momentum, V(q) the potential energy, K(t) the memory kernel of the GLE, and ζ (t) a random Gaussian distributed force. To impose a canonical sampling, the time correlation function H(t) = ζ (t)ζ (0) of the random force ζ (t) has to be related to the memory kernel K(t) by means of the fluctuation-dissipation theorem 31,32 However, if one wants to turn the GLE thermostat into a GLE quantum thermostat, which can be used as a tool to introduce NQEs in MD simulations, the constraint imposed by the fluctuation-dissipation theorem has to be relaxed to provide the necessary freedom to include harmonic NQEs. To that end, an effective temperature T , given by is introduced in order to match the probability distribution ρ of the displacement q of a classical harmonic oscillator with frequency ω and mass m, to the quantum mechanical probability distribution of the displacement (see Sec. S1), which is proportional to In this way, the classical canonical sampling imposed by the Boltzmann distribution [Eq. (5)] can be transformed into a quantum mechanical sampling that is in accordance with the probability distribution of the quantum harmonic oscillator [Eq. (6)] and as a consequence introduces a sampling of the harmonic NQEs. Therefore, one should picture the GLE quantum thermostat as a thermostat that excites every vibrational mode present in a molecular system to an effective temperature T , which is the temperature the mode is supposed to have according to the probability distribution of the quantum harmonic oscillator.
The appropriate canonical sampling is then once again obtained by imposing the fluctuation-dissipation theorem, where the effective, frequency-dependent temperature T now takes the role of the physical temperature T so that However, since the GLE quantum thermostat is a so-called non-equilibrium thermostat that violates the fluctuation-dissipation theorem (i.e., the fluctuationdissipation theorem for the physical temperature T), one now also needs to explicitly impose the equipartition theorem, as this cannot longer be taken for granted for a non-equilibrium thermostat. In view of an adequate modelling of harmonic NQEs, one therefore also has to satisfy the equipartition theorem for the quantum harmonic oscillator, which requires that The only difference between an ordinary GLE thermostat and the GLE quantum thermostat is hence the way in which the memory kernel K(t) and the random force ζ (t) are related to one another so that the GLE thermostat allows us to include harmonic NQEs in MD simulations without almost any additional computational cost. For a more detailed discussion of the GLE quantum thermostat, the reader is referred to Refs. 26, 33, and 34

B. Second order path integral MD
A more rigorous method to model NQEs in MD simulations is given by the path integral MD formalism. 35 In this formalism, which can be traced back to Feynman's path integral formulation of quantum mechanics, 36 NQEs are evoked by replacing every atom of the molecular structure by a so-called ring polymer that consists of P atom replicas, called beads, that interact with one another via a harmonic nearest-neighbour interaction, as depicted schematically in Fig. 2.
For an N-atom molecular system characterised by the Hamiltonian with p i and r i , respectively, the momentum and position of the ith particle and V(r 1 , . . ., r N ) the potential energy, this more rigorous introduction of NQEs can be readily obtained by means of a Trotter factorisation of the quantum mechanical canonical partition function Z = Tr[e − βĤ ], expanding the canonical density matrix e − βĤ as follows: with β = (k B T) −1 the inverse temperature and P the number of beads. The resulting expression for the quantum mechanical canonical partition function Z then becomes (see Sec. S2) where the Hamiltonian of the N-atom system is replaced by the Hamiltonian of an extended molecular system given by with and with r i , respectively, the position and momentum of the kth bead of the ith particle, ω P = √ P β the angular frequency of the harmonic interaction between the neighbouring beads of the ring polymer, and r The higher the number of beads P, the smaller the error in the Trotter factorisation of Eq. (10), yielding quantum mechanically exact results in the limit of an infinite number of beads. PIMD is therefore not only more accurate than the GLE quantum thermostat, as it takes both harmonic and anharmonic NQEs into account, but it also allows us to systematically increase the accuracy with which NQEs are modelled by increasing the number of beads, a desirable feature that is lacking for the GLE quantum thermostat.
The major downside of PIMD when compared to the GLE quantum thermostat is the considerable computational overhead that is induced by substituting every atom in the system by multiple bead replicas. Given the efficiency of the GLE quantum thermostat in modelling harmonic NQEs, it is however possible to develop a hybrid PI + GLE scheme, 29,30 combining PIMD with a GLE quantum thermostat in order to reduce the number of beads required to converge PIMD simulations. In the so-called PIGLET thermostatting scheme, 30 the potential and kinetic energy of the system are enforced to match the harmonic quantum limit by subjecting the path centroid, defined as to classical thermostatting, while the individual beads are subjected to a GLE quantum thermostat. As a consequence, the required number of beads to converge the energy is reduced by a factor of about four since the GLE quantum thermostat already takes care of the harmonic NQEs so that higher bead multiples are only necessary to converge the remaining anharmonic NQEs. 13 Besides PIGLET thermostatting, there also exist other cost reducing techniques for PIMD, 37 such as multiple time stepping (MTS), 38 ring polymer contraction (RPC), 28,39 and perturbed path integrals (PPI). 40,41 In the MTS algorithm, an algorithm first devised for classical MD, a more efficient force evaluation is obtained by separating the slowly varying longrange interactions from the rapidly varying short-range interactions. Given their different characteristic time scales, the long-range interactions only need to be evaluated at fixed time multiples of the short-range interactions, which are calculated at every time step. A similar splitting of the interaction is used in RPC, where the long-range part of the potential is evaluated on a contracted ring polymer containing fewer than P beads. In the limit of contraction to a single bead, the classical long-range force evaluation is retrieved. A different approach is used in PPI, where the theory of Feynman path integrals is combined with quantum mechanical perturbation theory to reduce the required number of beads by means of an a posteriori correction to the canonical partition function [Eq. (11)].
Finally, the convergence of PIMD simulations can also be accelerated by so-called high-order PI schemes, which rely on a more accurate decomposition of the canonical density matrix than the second-order Trotter factorisation [Eq. (10)]. A more accurate decomposition scheme such as the Suzuki-Chin (SC) scheme 42,43 then reduces the error to O(P −4 ), which also reduces the required number of beads to reach convergence. The force evaluation in SC PIMD is, however, more expensive than in Trotter PIMD, 44 as the equivalent of Eq. (14) not only depends on the physical potential energy, but also on the derivative of the physical potential energy. Nonetheless, SC PIMD becomes very advantageous at low temperatures or when a high accuracy is required. Therefore, this technique is used to validate our results, as outlined in Sec. IV.

III. COMPUTATIONAL DETAILS
As already mentioned in the Introduction, the importance of NQEs in MOF-5 will be assessed within this work by comparing ordinary, classical MD simulations to MD simulations including NQEs, using either a GLE quantum thermostat or PIMD. Given the rather large conventional unit cell of MOF-5, containing 424 atoms, an ab initio description of the PES is to be viewed as computationally too expensive, especially in combination with PIMD, so that all MD simulations within this work are performed using an in-house developed force field, generated from ab initio data using the in-house developed software package QuickFF. 45,46 In view of recent publications exemplifying a delicate interplay between a correct description of anharmonic force field contributions and NQEs, 14,47-49 the bond and bend contributions in the force field of MOF-5 are parametrised both by harmonic potentials, resulting in a "harmonic" force field, and by higher order polynomials, resulting in an "anharmonic" force field (see Sec. S3). For both types of force fields, the simulation results are discussed in Sec. IV.
All MD simulations are furthermore conducted in the NPT ensemble at a pressure of 1 bar, using a cut-off radius of 15 Å for the long-range interactions. For the classical and GLE quantum thermostatted MD simulations, the in-house MD code Yaff 50 is used in combination with LAMMPS, 51 which allows us to speed up the evaluation of the non-covalent force field contributions. The PIMD simulations on the other hand are performed with the iPI software package, 52,53 delegating the force evaluation to Yaff.
In the classical MD simulations, temperature and pressure are, respectively, controlled by a three bead Nosé-Hoover  52 The relaxation times of the different thermostats and barostats are set to, respectively, 0.1 ps and 1 ps. 59 For all simulations, except for the ordinary PIMD simulations using a white noise Langevin thermostat, an MD integration step of 0.5 fs is used. The ordinary PIMD simulations require a smaller integration step of 0.25 fs due to the presence of artificial high-frequency modes introduced by the harmonic nearest neighbour coupling. 28 The total simulation time of all PIMD simulations is 1.25 ns, whereas the classical and GLE quantum thermostatted simulations have a total simulation time of 3 ns. The equilibration time is set to 150 ps. Finally, to fix the number of replicas in the PIMD simulations, the convergence of the volume was tested as function of the number of beads (Fig. S1), reaching convergence at, respectively, 8 beads for PIGLET thermostatted PIMD simulations and 32 beads for ordinary PIMD simulations.

A. Structural properties: RDFs
To evaluate the influence of NQEs upon the structural properties of MOF-5, a series of radial distribution functions (RDFs) is constructed, focussing on two different atom pairs. On the one hand, we consider C− −H pairs of atoms, as the C− −H bonds in MOF-5 have the highest characteristic frequency, i.e., about 3195 cm −1 , as obtained from a normal mode analysis performed with TAMkin. 60 On the other hand, we also consider Zn− −O pairs of atoms, as the Zn− −O bond in MOF-5 is characterised by a much lower characteristic frequency of about 450 cm −1 and the low-frequency modes involving ZnO 4 tetrahedra are believed to contribute to the negative thermal expansion of MOF-5. 20 Based on their characteristic frequency, these bonds can also be attributed a so-called characteristic temperature or Debye temperature which can be used as a measure to quantify the importance of NQEs, as it represents the temperature at which the quantum mechanical energy quantum ω equals the thermal energy quantum k B T and NQEs are only expected to be of importance for temperatures T Θ D . For the C− −H and Zn− −O bonds in MOF-5, the characteristic temperatures are, respectively, given by 4600 K and 650 K. For both types of atom pairs, the RDFs are shown in Fig. 3, comparing classical MD simulations to MD simulations including NQEs at different temperatures, using the anharmonic force field of MOF-5. The first peak of each RDF is always located at the equilibrium bond distance of the atom pair in MOF-5. The atom pairs contributing to the other RDF peaks are indicated on Fig. S2. For the C− −H pairs of atoms, the equilibrium bond length is observed to increase by an amount of about 0.02 Å when including NQEs (see Table S2 and Fig. S6), shifting the first RDF peak towards slightly higher equilibrium distances. This shift is also present for the more distant RDF peaks, indicating that the molecular structure as a whole is affected by the presence of NQEs. Furthermore, the RDF peaks are also observed to be much broader when modelled quantum mechanically instead of classically, yielding widths that are more than twice as large (see Fig. S3). In addition, the quantum mechanical widths appear to be more or less temperature independent so that, in contrast to the classical description, the C− −H stretch is not yet thermally active at temperatures below 500 K. Quantum mechanically, the behaviour of the C− −H stretch is therefore dominated by the ZPE it possesses, as it is significantly larger than the thermal energy up to temperatures as high as Θ D /2 = 2300 K (within a harmonic approximation).
For the Zn− −O pairs of atoms, the differences between the RDFs with and without NQEs are less pronounced, but nevertheless still present. Given the much lower characteristic temperature of the Zn− −O stretch in comparison with the C− −H stretch, the ZPE associated with the Zn− −O bond only dominates the thermal energy up to a temperature of about 300 K. Contrary to the C− −H pairs of atoms, the first RDF peak of the Zn− −O pairs of atoms shows no appreciable shift when including NQEs (see Table S3 and Fig. S6), as it is comprised of two overlapping RDF peaks corresponding to the Zn− −O 1 and Zn− −O 2 bonds [see Fig. 4(a)], limiting the peak's resolution. If the probability distribution of the Zn− −O 1 bond is however isolated, as done in Fig. 4(b), one does observe a minor shift in the equilibrium bond length, just as for the C− −H atom pairs. The different RDF peaks are furthermore also observed to broaden under the influence of NQEs, albeit by less than half their classical width (see Fig. S4).
When comparing the RDFs obtained using different methods to model NQEs, as shown for GLE quantum thermostatted MD simulations and PIGLET thermostatted PIMD simulations in the right column of Fig. 3, the different methods are found to agree well with one another. For the Zn− −O pairs of atoms, there is almost no observable difference between the RDFs, whereas for the C− −H pairs of atoms only minor differences are present, which are limited to the first RDF peak. The GLE quantum thermostat therefore performs remarkably well in comparison with PIMD, given its approximate nature and low computational cost. In order to benchmark the accuracy of the PIGLET thermostatted PIMD simulations, a 16 bead SC PIMD simulation was performed at 300 K, 44 confirming the validity of the results obtained with PIGLET thermostatting, as shown in Fig. S5.
Performing the above mentioned simulations with a harmonic force field of MOF-5 instead of an anharmonic force field leads to similar results (see Sec. S5). The inclusion of anharmonic contributions in the force field mainly causes the RDF peaks to shift towards larger interatomic distances, as the equilibrium bond lengths are larger for anharmonic potentials with the same force constant, as illustrated by Fig. 5. When including NQEs, this effect becomes even more pronounced for the bonds with a high characteristic frequency.  (Table S3). For the C− −H RDFs, by contrast, the difference in peak positions is about 0.02 Å for the first peak at 100 K, a value that further increases both with the peak number and the temperature, reaching a value of 0.06 Å for the sixth peak at 500 K (Table S2), a value that is two times larger than in classical simulations. This correlated effect between the inclusion of anharmonicities in the bond description and an adequate modelling of NQEs can be inferred from the quantum mechanical description of an anharmonic potential such as the well-known Morse potential. In comparison with a harmonic potential, the anharmonic potential systematically yields a larger equilibrium bond distance, even at the lowest vibrational energy level, where the bond only possesses a ZPE (Fig. 5). As higher vibrational energy levels of the anharmonic potential are occupied, the equilibrium bond distance increases up to a point where the occupation of additional vibrational energy levels no longer affects the equilibrium bond distance, explaining the observed temperature dependence of the RDF peak shifts. The more distinct differences for the C− −H RDFs can be related to the higher characteristic frequency of the C− −H bond, resulting in a larger ZPE and a larger separation between the vibrational energy levels. The higher the peak number, the larger the peak shift due to the accumulated effect of increased equilibrium bond lengths for the intermediate atoms. This coupling between the description of anharmonicities and NQEs will also prove to be important in an accurate description of the negative thermal expansion of MOF-5, as explained in Sec. IV B. Apart from RDFs, the structural properties of MOF-5 can also be characterised by other means such as probability distributions of internal coordinates including (dihedral) angles, X-ray diffraction (XRD) patterns, or pore volume sizes. An example of an interesting internal coordinate is the C 1 − −C 2 − −C 2 angle [ Fig. 4(d)], which can be viewed as a measure for the transverse translational motion of the benzene ring, an important motion in the negative thermal expansion of MOF-5. 20 Just as for atomic bonds, the probability distribution of the of Chemical Physics angle is observed to have a different mean value and to become noticeably broader when including NQEs, especially at lower temperatures. The XRD pattern by contrast shows no appreciable differences (Fig. S7), as its resolution is less adequate to identify smaller structural changes. Similarly, also the pore

B. Thermal expansion
The second property of MOF-5 for which the importance of NQEs is examined is the negative thermal expansion of MOF-5. 20-25 Quantitatively, the thermal expansion can be characterised by the volumetric thermal expansion coefficient α V , defined as the relative change in volume with respect to temperature at constant pressure, that is, Within first order, the volumetric thermal expansion coefficient α V is then obtained as the slope of a linear fit to the logarithm of the volume as a function of temperature. For periodic framework materials such as MOF-5, one can also define a lattice thermal expansion coefficient α a , which is related The Journal of Chemical Physics ARTICLE scitation.org/journal/jcp to the volumetric thermal expansion coefficient by α V = 3α a , assuming an isotropic thermal expansion. In Fig. 6, the temperature dependence of the volume of MOF-5 is shown for (PI)MD simulations performed with both the harmonic and anharmonic force fields, for which the corresponding lattice thermal expansion coefficients are listed in Table I. As all the thermal expansion coefficients calculated with the anharmonic force field lie within the experimental range, including the value derived from classical MD simulations, NQEs appear to be only of secondary importance in modelling the thermal expansion of MOF-5. A proper description of the anharmonicities of the underlying PES is observed to be far more important, as evidenced by the simulation results obtained with the harmonic force field of MOF-5. Also in the harmonic case, NQEs only yield small corrections to the classically obtained results, which does however partly correct for the overestimation of the lattice thermal expansion coefficient in comparison with experimental values (see Table I and Fig. S11). As a consequence, more accurate MD simulations accounting for NQEs are only useful if the PES is modelled sufficiently accurately, as the most decisive contribution to the thermal expansion stems from an adequate modelling of the anharmonicities, reducing the thermal expansion coefficient by about 20%-30% for MOF-5.
The absolute volume on the other hand is significantly affected by the presence of NQEs in an anharmonic description of the atomic bonds, yielding differences of the order of 150 Å 3 with respect to the classically predicted volume, an effect which is completely absent within a harmonic description of the atomic bonds (Fig. 6). To explain this interplay between NQEs and anharmonicities, the equilibrium volume at 0 K is determined by means of a series of fixed volume geometry optimisations, as shown in Fig. 7. NQEs are then only accounted for within the harmonic limit, yielding a ZPE contribution of k ω k 2 , with ω k the normal mode  frequencies as determined with TAMkin. The resulting equilibrium volumes, with and without NQEs, are indicated in Fig. 6 and are observed to agree well with the extrapolated (PI)MD simulation data so that the explanation for the difference in volumes is to be found within the phenomenon of ZPE.
In the harmonic limit, the ZPE only depends on the normal mode frequencies of the structure so that it is instructive to determine how those frequencies depend on the volume of the structure. For a harmonic potential, the normal mode frequencies as given by the second order positional derivative of the potential are volume independent, yielding a constant ZPE contribution as a function of volume (see Sec. S6). An anharmonic potential by contrast has a second order positional derivative that depends on the position, due to the presence of higher order polynomial contributions to the potential, and is hence characterised by a set of volume dependent normal mode frequencies. As a consequence,   (Fig. S8), resulting in an increase of the corresponding ZPE contributions. As the total increase in ZPE is about one order of magnitude larger than the corresponding increase in thermal energy (i.e., 41.2 kJ/mol versus 4.0 kJ/mol), the difference in volume due to the presence of NQEs is perceived within a wide temperature range. Only at very high temperatures, classical and quantum mechanical volumes can be observed to converge to one another (Fig. S10). The volume dependence of the ZPE of MOF-5 in the harmonic force field (Fig. 7) is to be attributed to other energy contributions in the force field that have an anharmonic character, such as the non-covalent electrostatic and van der Waals contributions. With decreasing volume, or equivalently increasing temperature, the normal mode frequencies are now observed to decrease (Fig. S8) so that also the ZPE decreases with temperature and the thermal energy becomes the dominant energy contribution. The difference in volume due to the presence of NQEs is therefore only noticeable in a relatively small low temperature range, in contrast with the anharmonic case.
To conclude, we also evaluate the different methodologies used to model NQEs in the thermal expansion of MOF-5. Just as for the RDFs, the GLE quantum thermostat performs well when compared to the more demanding PIMD simulations, yielding similar absolute volumes and a matching volumetric temperature dependence. Only at the lowest temperatures, below 100 K, the GLE quantum thermostat predicts a deviating behaviour from the PIMD results, as for the GLE thermostat the volume is observed to saturate. Given the approximate nature of the GLE quantum thermostat, modelling only harmonic NQEs, this saturation of the crystal volume can be understood as an artefact of the constant harmonic ZPE, which causes the bond length fluctuations to become constant and thus saturates the crystal volume.

V. CONCLUSIONS
In this article, we examined the influence of anharmonicities and NQEs upon an accurate modelling of the structural and thermal properties of MOF-5. The structural properties, probed by means of RDFs, were shown to be significantly affected by the presence of NQEs. Under the influence of quantum fluctuations, the RDF peaks tend to broaden and show furthermore a small shift towards larger interatomic distances, due to a coupling between NQEs and the anharmonicities of the PES. Also for the temperature dependence of the RDFs and hence, the structure itself, a different picture emerges. At lower temperatures, the ZPE becomes the most dominant energy contribution and as a consequence, reduces the importance of the thermal energy. For the highest characteristic frequencies in MOF-5, this effect is perceivable up to very high temperatures, making the C− −H bond virtually temperature independent below 500 K. For the negative thermal expansion of MOF-5, by contrast, the inclusion of NQEs appears to be only of secondary importance with respect to a proper modelling of the anharmonicities of the PES, as it yields only small corrections to the thermal expansion coefficient. The interplay between NQEs and anharmonicities does however affect the absolute volume, as volume dependent characteristic frequencies lead to a volume dependent ZPE, resulting in a different equilibrium volume. A sufficiently accurate description of the PES is therefore prerequisite when modelling NQEs.
A comparison of different methodologies capable of modelling NQEs furthermore showed an excellent agreement between the approximate technique of GLE thermostatting and the more rigorous framework of PIMD, at least for the structural properties and the negative thermal expansion of MOF-5. Therefore, the inclusion of NQEs in MD simulations of MOF-5 can be achieved most efficiently by means of a GLE thermostat, which produces almost no computational overhead.

SUPPLEMENTARY MATERIAL
In the supplementary material, a more elaborate discussion of the GLE quantum thermostat and PIMD can be found, as well as additional information concerning the RDFs and thermal expansion of MOF-5. Also the convergence of the PIMD simulations and the force field parametrisation are addressed, providing a force field parameter input file for MOF-5.