Modeling the structural and thermal properties of loaded metal-organic frameworks.An interplay of quantum and anharmonic fluctuations

Modeling the structural and thermal properties of loaded metal-organic frameworks.
An interplay of quantum and anharmonic fluctuations

Venkat Kapil Laboratory of Computational Science and Modelling, Institute of Materials, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland    Jelle Wieme Center for Molecular Modeling, Ghent University, Tech Lane Ghent Science Park Campus A, Technologiepark 903, 9052 Zwijnaarde, Belgium    Steven Vandenbrande    Aran Lamaire    Veronique Van Speybroeck Center for Molecular Modeling, Ghent University, Tech Lane Ghent Science Park Campus A, Technologiepark 903, 9052 Zwijnaarde, Belgium    Michele Ceriotti Laboratory of Computational Science and Modelling, Institute of Materials, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
July 20, 2019

Metal-organic frameworks show both fundamental interest and great promise for applications in adsorption-based technologies, such as the separation and storage of gases. The flexibility and complexity of the molecular scaffold poses a considerable challenge to atomistic modeling, especially when also considering the presence of guest molecules. We investigate the role played by quantum and anharmonic fluctuations in the archetypical case of MOF-5, comparing the material at various levels of methane loading. Accurate path integral simulations of such effects are made affordable by the introduction of an accelerated simulation scheme and the use of an optimized force field based on first-principles reference calculations. We find that the level of statistical treatment that is required for predictive modeling depends significantly on the property of interest. The thermal properties of the lattice are generally well described by a quantum harmonic treatment, with the adsorbate behaving in a classical but strongly anharmonic manner. The heat capacity of the loaded framework – which plays an important role in the characterization of the framework and in determining its stability to thermal fluctuations during adsorption/desorption cycles – requires, however, a full quantum and anharmonic treatment, either by path integral methods or by a simple but approximate scheme. We also present molecular-level insight into the nanoscopic interactions contributing to the material’s properties and suggest design principles to optimize them.


si \altaffiliationThese authors contributed equally to this work \altaffiliationThese authors contributed equally to this work

1 Introduction

Tailor-made porous materials 1 like metal-organic frameworks (MOFs) 2 are at the core of emerging technologies due to their exceptional physical and chemical properties such as a tunable ultrahigh porosity and an associated enormous gas storage capacity. Therefore, they have been proposed for applications such as adsorbed natural gas (ANG) storage in vehicles 3, 4, adsorption-driven heat pumps 5, 6, and carbon capture and sequestration (CCS) 7, 8. While a lot of work still needs to be done to optimize the crucial adsorption and storage properties of these porous materials, 9 studies on other critical requirements such as heat management are gaining interest. 10 For instance, the heat capacity, i.e., the amount of energy required to increase the material’s temperature, is a fundamental thermodynamic property of interest in these applications which involve large thermal fluctuations as adsorption and desorption processes imply the release or consumption of energy. Moreover, the heat capacity of the MOF affects the energy penalty to regenerate the adsorbent in, for example, CCS. 11 To date, however, information on the heat capacity is lacking for most MOFs 12, 13 and the influence of adsorbed guest molecules on the heat capacity has not yet been investigated, in contrast to other thermal properties such as the thermal expansion behavior and the thermal conductivity. 14, 15, 16

Within this context, an efficient and accurate simulation protocol to tackle the structural and thermal properties of MOFs including all the relevant physical effects could facilitate a better understanding of the structure-property relations and suggest design principles for materials with improved properties. Due to the importance of finite-temperature effects, anharmonicity, and nuclear quantum effects (NQEs), the modeling of the thermophysics of MOFs is generally not a trivial exercise. The first two effects have already been the subject of many investigations and were included in our protocol to characterize the thermodynamics of MOFs. 17, 18 Furthermore, very recently, some of the present authors highlighted the necessity of an accurate theoretical framework for the design of thermoresponsive MOFs. 19 However, the impact of NQEs has so far received far less attention within the MOF community, despite the many light atoms contained in the crystal structure and present inside the pores. 20, 21

In this regard, path integral molecular dynamics (PIMD) 22 provides an ideal reference framework for the evaluation of thermodynamic averages, as it seamlessly captures both NQEs and the anharmonic motion of nuclei. The statistics of distinguishable quantum particles can be obtained through the equivalence between the thermodynamics of a quantum system and a classical ring polymer containing replicas of the system. 23 In the limit of large values, NQEs can then systematically be accounted for. The major downside of this technique is the associated high computational cost, i.e., times the cost of classical molecular dynamics (MD). However, several methodological advances 24, 25, 26, 27, 28, 29 that enable a reduction of the computational cost have made it a mainstream technique for material modeling. 30

Figure 1: The structure of the MOF-5 with a gas loading of 0, 50, 100, and 150 molecules of methane (left to right) in the conventional unit cell. (8(ZnO(CO)) CH) The oxygen, carbon, and zinc atoms are shown in red, silver, and blue respectively. For the sake of aesthetics the hydrogen atoms are not included. The methane molecules are represented by silver tetrahedra.

An additional difficulty arises from the fact that most experiments and practical applications are performed in isothermal-isobaric conditions, while the vast majority of atomistic simulations are performed with a fixed unit cell, corresponding to isochoric conditions. As most solid materials have a very small compressibility, the difference between the two ensembles is often negligible. For MOFs on the contrary – particularly when loaded with a gas – the behavior in isobaric and isochoric conditions can be very different. Some of us emphasized the importance of taking into account the variations of the cell shape to simulate properties of flexible MOFs. 17, 31 While algorithms for performing path integral simulations at constant pressure conditions exist 32, 33, an accurate evaluation of the thermophysical properties requires a very large number of replicas for convergence. In this article, we introduce a method, based on a recently developed implementation of high order path integral factorizations 26, to greatly accelerate the convergence of these simulations.

This method, in combination with a first-principles-based force field 34, 35, 18, makes it possible to characterize the structural and thermophysical properties of complex molecular systems such as guest-loaded MOFs. We investigate the archetypical case of the well-known MOF-5 36, 37 in the presence and absence of methane in its pores (see Figure 1). Evaluating and understanding the impact of methane adsorption on the properties of MOFs is especially important as they have been proposed as potential adsorbents for natural gas storage applications. 3, 9, 4 We demonstrate the crucial role of a complete statistical-mechanical description of the quantum and anharmonic fluctuations in MOFs for a correct description of structural properties and the heat capacity of guest-loaded MOFs. By meticulously disentangling anharmonic and nuclear quantum effects for both the lattice and the guest particles, we are able to propose an efficient empirical calculation scheme which may be used to screen MOFs with beneficial thermal properties on a larger scale.

2 Methodology

2.1 Materials

The materials that are considered in this theoretical work are pristine and methane-loaded MOF-5 scaffolds. 36 This framework consists of ZnO(CO) inorganic nodes connected through 1,4-benzenedicarboxylate (bdc) linkers. The unit cell is cubic and contains eight inorganic nodes, as shown in Figure 1. We consider three different loadings of , , and methane molecules in the conventional unit cell (8(ZnO(CO)) CH), which encompasses both the low- and high-adsorption regime. 38, 3 At 100 bar, for example, approximately 120 methane molecules are present per conventional unit cell, as measured by Mason et al. 3 (see SI Figure S2).

2.2 First-principles-derived force fields

The molecular simulations are performed using newly developed force fields for MOF-5 and methane. They are derived with QuickFF 34, 35, a software package developed to derive force fields for MOFs in an easy yet accurate way based on information obtained from first-principles input data. Isolated cluster models were used to generate the required first-principles input data, which includes the geometry and the Hessian in equilibrium together with the atomic charges. Within the QuickFF protocol, the quantum mechanical potential energy surface (PES) is approximated by a sum of analytical functions of the nuclear coordinates that describe the covalent and noncovalent interactions. The covalent interactions, which mimic the chemical bonds between the atoms, are approximated by different terms as a function of the internal coordinates (bonds, bends, out-of-plane distances, and dihedrals). The noncovalent interactions are composed of electrostatic and van der Waals interactions. The guest-host interactions between MOF-5 and methane only include noncovalent terms. A detailed discussion of the force field energy expression and derivation is provided in the Supporting Information (see SI Section S4).

2.3 Thermodynamic ensembles

In order to study the classical and quantum isobaric heat capacity of the guest-loaded MOF scaffold, the classical and quantum isothermal-isobaric thermodynamic ensembles – as defined in Ref. 17 – are respectively used. For a system with particles, subjected to an external mechanical stress , with the hydrostatic pressure and the deviatoric stress, the classical and quantum isothermal-isobaric ensembles at temperature and fixed normalized cell tensor are described by the partition functions


where are the corresponding canonical partition functions at volume , normalized cell tensor , and temperature . Similarly, the flexible partition functions, which allow the cell shape to change, are defined by the partition functions 39:


The classical ensembles are sampled by classical molecular dynamics 40, whereas the quantum ensembles are sampled using path integral molecular dynamics 22. The latter maps the quantum partition function of a system to that of a classical ring polymer Hamiltonian made of replicas of the system. 23 Formally, exact results can be obtained in the limit of . Unfortunately, the requisite value of to calculate structural properties rises rapidly with decreasing temperature 41 and for properties such as heat capacity 42 the increase in is even greater. This makes the standard scheme prohibitively expensive.

2.4 Accelerated simulation scheme

To remedy this problem, we present an accelerated scheme based on a constant pressure integrator for the high order path integral method 24. The quantum statistics of distinguishable particles arises from the non-commutative nature of the potential and kinetic energy operators. In standard (second order) path integral schemes, an approximate factorization of the high-temperature Boltzmann operator is introduced, that leads to an error that decreases as . High order techniques use an alternative splitting of the Boltzmann operator 43, 44, leading to an alternative ring polymer Hamiltonian with a faster, convergence to the exact quantum limit. This makes it possible to reduce the number of replicas and hence the computational cost. 24 While many high order schemes exist 43, 29, here we focus on the specific case 45 of a fourth-order Suzuki-Chin (SC) splitting 46, 44, which yields a so-called SC Hamiltonian . Considering for simplicity the case of a single particle in an external potential, the SC Hamiltonian takes the form:



is the ring polymer Hamiltonian of a free particle, subjected to cyclic boundary conditions (), and


Note that the odd and even replicas feel the physical potential scaled by factors of and respectively and that the high order term , that depends on the modulus of the force , only acts on the even replicas.

Figure 2: Fractional error in the standard (red) and Suzuki-Chin (blue) path integral estimators of the energy (circles) and the heat capacity (triangles) of empty MOF-5 modeled by the corresponding Debye crystal potential, as a function of the number of beads at 100 K (bottom) and 300 K (top). The values of the energy and the heat capacity were obtained analytically 47.

The improved efficiency of the high order scheme is demonstrated by studying the convergence of the total energy and its temperature derivative for a harmonic model of MOF-5 – obtained by computing the dynamical matrix associated with the first-principles-derived force field described in the previous section 2.2 – with respect to the number of replicas. The analysis is done analytically at two different temperatures: K and K, as described in Ref. 47. As shown in Figure 2, the high order scheme improves the convergence of both the energy and the heat capacity by a factor of two at room temperature. In the low temperature or high accuracy regime where quantum effects dominate, the efficiency of the high order scheme is even more significant.

The high accuracy afforded by the SC scheme makes it particularly useful to compute the heat capacity. This far, however, it has only been successfully applied to relatively simple systems 48, 49, 41 and heat capacities in particular have only been reported for small clusters of molecules and constant-volume conditions. 50 The difficulty in applying fourth-order schemes to complex materials can be understood by considering the fact that the force and virial contain derivatives of with respect to the atomic positions and the cell parameters,


Given that already contains first-order derivatives of the physical potential, the computation of the forces and the virial, required to sample the isothermal-isobaric ensemble by means of path integral dynamics, also demands the evaluation of higher-order derivatives of the potential, which is often cumbersome and computationally prohibitive. Much of the work on the practical implementation of high order path integrals has therefore focused on avoiding the calculation of these terms, by sampling the standard path integral Hamiltonian and introducing fourth-order statistics by re-weighting 24, 49. Unfortunately, re-weighting schemes have poor statistical performance for large systems 24, 51, so the application of the SC scheme has until now been limited to small systems and to constant-volume sampling.

To circumvent these limitations, we evaluate the virial using a finite-difference scheme that has recently been introduced by some of the present authors, to sample the SC canonical partition function directly. We discuss the essential ingredients of this scheme in Appendix A and report a detailed derivation in the SI. Direct access to the instantaneous force and virial allows us to develop an integration scheme that samples the quantum isothermal-isobaric ensemble. This scheme can be regarded as an extension of the second-order scheme introduced in the pioneering work of Martyna and co-workers 32. The rather cumbersome derivation and the equations of motion are given in the SI. It should be noted that our implementation in i-PI 52 is also fully compatible with multiple time stepping 53, 28 and stochastic thermostatting, extending the integrators introduced in Ref. 54 to the isothermal-isobaric ensemble.

By using the finite-difference expressions, only additional force evaluations are needed instead of Hessians. It is also useful to note that these components have a prefactor, which means that they become small and slowly varying for typical values of used in SC simulations. 26 This facilitates the use of a long time step for integrating the high order forces and virials, enabling us to sample the ensemble while evaluating the expensive terms rather infrequently.

To further reduce the computational effort, we combine this scheme with other accelerated PIMD methods that rely on the separation of the total potential into a cheap short-ranged term and an expensive long-ranged term, which is the case in our study, as discussed in subsection 2.2. We consider in particular ring polymer contraction (RPC) 27 and multiple time stepping (MTS) 53. Within the first method, the long-ranged components can be computed separately on a smaller ring polymer of beads, which are subsequently, without loss in accuracy, extrapolated to the case of beads. The second method reduces the frequency at which the long-ranged interactions are computed by the use of a longer time step for the integration. The two methods can also be used together to achieve substantial computational savings 28, 55 and can be seamlessly combined with our high order constant pressure scheme.

2.5 Calculation of thermodynamic observables

Within the SC scheme, the thermodynamic averages of all structural observables are estimated by an ensemble average over the odd replicas:


as first demonstrated by Jang and Voth 24. These estimators are commonly referred to as “operator” (OP) estimators, as opposed to the “thermodynamic” estimators obtained by derivation of the path integral partition function, that have often pathological statistical behavior. The simplicity of OP estimators makes the SC scheme very appealing in comparison to other high order schemes, in which ad hoc estimators need to be constructed for simple structural observables 49, 56. The simple OP estimators for the total energy and enthalpy are listed in Appendix B. Since the standard estimators for the heat capacity in path integral methods tend to be very complex and exhibit a large variance, we derive an OP double-virial estimator for the isobaric (and isochoric) heat capacity . The derivation is presented in the SI (Section S3.2) and the resulting expression is given in Appendix B, where it is also shown that this estimator has very good statistical properties and outperforms existing heat capacity estimators 57. However, in this study we computed thermophysical properties over a broad range of temperatures and found it more convenient to estimate by means of a finite difference approximation to the temperature derivative of the enthalpy:


A dedicated estimator will prove useful in simulations that are targeted at a single, specific temperature.

3 Computational Details

The required first-principles cluster data 58 for the determination of the covalent terms in the force field are generated with Gaussian 16 59 using the B3LYP 60, 61, 62 exchange-correlation functional. A 6-311G(,) basis set 63 is used for the C, O, and H atoms, together with the LanL2DZ basis set for Zn. 64 The atomic charges are derived with the Minimal Basis Iterative Stockholder (MBIS) partitioning scheme 65. The atomic charges of the MOF-5 clusters are obtained from the PBE 66 electron density computed with GPAW 67, 68 as implemented in Horton 69. For methane, the atomic charges are derived from the B3LYP all-electron density obtained with Gaussian 16. The parameters of the van der Waals interactions are taken from the MM3 force field. 70, 71. The van der Waals interactions are calculated up to a cutoff of 15 Å and a tail correction is added to the potential and its derivatives. 72 The initial configurations of the methane molecules are generated using RASPA 73 by inserting methane molecules at random positions, while ensuring that only realistic intermolecular distances are retained. Afterwards a canonical Monte Carlo algorithm was used to equilibrate the positions.

For MOF-5 with and without methane, we perform classical and path integral MD simulations at a mechanical pressure of 1 bar and at different temperatures in the range of 100 K to 500 K. The classical MD simulations of both loaded and pristine MOF-5 are performed using Yaff in the ensemble, i.e. with no constraints on the unit cell. 17 While the covalent interactions are calculated by Yaff, the expensive long-range interactions are computed by lammps 74 in a computationally efficient manner. The equations of motion are updated via a Verlet scheme, with a time step of 0.5 fs. The temperature is controlled via a single Nosé-Hoover chain consisting of three beads, with a relaxation time of 100 fs. 75, 76, 77 A Martyna-Tobias-Klein barostat with a relaxation time of 1000 fs is used to control the pressure. 78, 79, 17 We performed five independent runs of 500 ps, starting from a different random seed and from different methane positions. For the empty MOF-5, a single trajectory of 500 ps was used. An equilibration time of 100 ps was considered.

The PIMD simulations are performed with the universal force engine i-PI 33, 52 in the ensemble, where the cubic symmetry is kept fixed. The evaluation of the forces is carried out by Yaff and lammps, similar to the classical MD simulations, while the time evolution of the nuclei to sample the appropriate thermodynamic ensemble is done with i-PI. To control the temperature, a PILE-L thermostat 80 is applied to the system and a white noise Langevin thermostat 81 is applied to the cell. To control the pressure, a path-integral version of the Bussi-Zykova-Parinello (BZP) barostat 82, 33, adapted to the SC scheme, is used. The time constants for the thermostats and the barostats are same as the ones used in the classical simulations. A type 54 MTS scheme 53 (see SI Section S2.3) is used to integrate the equations of motion. The computationally cheap short-range terms of the force field are computed on 64 replicas and integrated with a time step of 0.25 fs. The remainder of the interactions, i.e. the expensive long-range interactions, are computed on replicas using RPC and integrated with a time step of 1 fs. As discussed above, a finite differences strategy is adopted to determine the heat capacity from the enthalpy with a temperature interval of 25 K. We performed thirty independent runs of 50 ps, starting from a different random seed and from different methane positions. For the empty MOF-5, five independent trajectories of 125 ps were used. An equilibration time of 25 ps was considered.

Figure 3: Panels (a) and (b) show the lattice parameter of MOF-5 with molecules of methane as a function of temperature (), obtained from classical MD and PIMD respectively. Panel (c) shows the linear thermal expansion coefficient () as a function of . The classical and quantum estimates are respectively shown with dashed and solid lines. Error bars indicate statistical uncertainity.

4 Results

Having developed an accelerated integration scheme for the quantum isothermal-isobaric ensemble, we now focus on the structural and thermal properties of methane-loaded MOF-5. Extensions towards other MOFs and adsorbates will be the topic of future studies. The importance of the inclusion of NQEs and anharmonicities in the modeling of the heat capacity is probed by comparing the results with other methods such as classical MD, which neglects NQEs, and the harmonic approximation, which neglects anharmonicity. We discuss the accuracy of these commonly-adopted approximations and provide empirical relations, which might resolve the general lack of knowledge on the heat capacity of this class of materials.

4.1 Structural properties

To unravel the influence of adsorbates on the framework and finally on the heat capacity, we start by investigating the structural response of MOF-5 for various loadings and temperatures. Here, one could expect that a proper inclusion of NQEs already becomes important as zero-point effects were recently found to substantially increase the volume of MOF-5 when comparing classical MD with PIMD. 83 Additionally, NQEs have previously been observed to change the volume of bulk alkanes by about . 84, 85, 86 A comparison of Figures 3 (a) and (b) indeed reveals that the inclusion of NQEs increases the volume by almost 1 % for all loadings and temperatures. Horizontally, this shift corresponds to a more substantial temperature reduction of about 100 K.

Interestingly, the qualitative ordering of the volume as a function of loading at the different temperatures does not change appreciably with or without the inclusion of NQEs. At low temperatures, the material slightly shrinks in the presence of methane. The observed adsorption-induced deformation can be understood by attractive van der Waals interactions between the framework and the adsorbed methane. 87, 88 At higher temperatures, by contrast, the empty framework has the lowest volume, as entropic and kinetic effects start dominating and the adsorbed molecules increase the internal pressure, which leads to a volumetric expansion when increasing the loading. The main effect of NQEs on the volume of the guest-loaded system is thus an upward volume shift, which is to a large extent independent of the number of guest molecules and the temperature, and thus originates from the zero-point fluctuations of the framework lattice. 83

Varying the concentration and the type of adsorbates in the framework was suggested by Calero and co-workers 15 as a way to control and tune the thermal expansion of a system based on classical MD simulations. We confirm that with methane, it is possible to go from the well-known negative thermal expansion behavior of the empty framework 89, 90, 91 towards positive thermal expansion. A proper inclusion of NQEs in our molecular dynamics simulations does not influence the predictions in this temperature window. The role of the quantum mechanical nature of the framework nuclei also remains limited when studying the thermal expansion coefficient, as both classical MD and PIMD simulations lie within the experimental range. 83 This conclusion does not change with methane in the pores, as shown in Figure 3 (c).

Figure 4: The methane distribution in the pores of MOF-5 at different temperatures as obtained from PIMD simulations. Orange spots indicate high probability adsorption sites. Other colors show the distribution of the low probability methane positions in the conventional unit cell and represent the probability representation (from very high (orange), to high (dark blue), to low (white) probability).

A more surprising picture emerges when looking at the distribution of methane inside MOF-5. Recent PIMD simulations 85 of bulk methane (at 110 K) have shown that NQEs lead to significant changes in the structure of methane at low temperature, corresponding to an overall softening of the structure and an increase in the intermolecular distance by about 0.1 Å. In contrast, in our study of methane confined in the pores of MOF-5, even at 100 K – where NQEs are expected to be the greatest – there is no appreciable difference between the shape of the classical and quantum distribution functions of methane, as shown in Figure 4. This can be understood by the fact that the change in the structure of bulk methane comes entirely from the isotropic expansion of the gas 92. In the case of methane molecules confined in the pores, the low compressibility of the framework makes the expansion as observed in bulk methane when including NQEs impossible.

This discussion shows that the structural response of MOF-5 to a varying number of methane molecules and temperature is largely unaffected by NQEs, except for the zero-point lattice fluctuations. Our observations also corroborate the common practice of ignoring NQEs when studying the loading of porous materials by Grand Canonical Monte Carlo simulations 73. Nevertheless, this conclusion cannot be generalized to other adsorbates, especially those possessing stronger intermolecular interactions such as hydrogen bonding.

Figure 5: Heat capacity of the empty MOF-5 as a function of temperature computed using classical MD (dashed), PIMD (solid) and the harmonic approximation (dotted). The right pointing arrow shows the Dulong-Petit limit. Different experimental results are shown in black using triangular 13, square 93 and diamond 12 markers. Error bars indicate statistical uncertainity.

4.2 Heat Capacity

MOF-5 has been the subject of a few experimental heat capacity studies 12, 94, 93, 13 which have shown that the material has a low specific (or molar) heat capacity, about 0.7 J/gK at room temperature, even when compared to other MOFs. Depending on the type of the application, a large (e.g. for ANG to limit temperature fluctuations) or a small (e.g. for CCS to limit the energy penalty) heat capacity is sought after. It is thus important to understand how this property changes at different levels of loading and temperature, and to determine the factors influencing the heat capacity, which is now possible for the first time using our high order PIMD scheme.

In the previous section, it has been shown that classical MD can – at least qualitatively – be used to model the structural response of MOF-5 in the presence of methane at various temperatures. This approach is however expected to fail for the description of the heat capacity since the heat capacity of many systems is dominated by NQEs at room temperature, as evidenced by experimental deviations from the classical Dulong-Petit law. The most common way of including NQEs for solids is the static harmonic approximation, using Einstein’s or Debye’s harmonic model for solids, which is able to reproduce the heat capacity of many solids and will therefore also be used for comparison.

We begin by presenting the estimates of the temperature dependence of the isobaric heat capacity of the empty MOF-5 framework. As shown in Figure 5, the classical MD estimates (dashed line) are in agreement with the Dulong-Petit law. The simulations yield an almost constant value of k per degree of freedom, which indeed results in large deviations from the experimental values. 12, 94, 93, 13 Upon inclusion of NQEs with our PIMD scheme (solid line), we find that the results follow the experimental measurements reasonably well up to almost 400 K. This agreement is remarkable as these measurements are typically carried out on the as-synthesized sample, which possibly includes solvents 12 and differs from the perfect crystal that we have simulated.

Figure 6: Panel (a) shows the comparison of the classical (dashed), quantum (solid), and harmonic estimates (dotted) of the isobaric heat capacity of MOF-5 with 100 molecules of methane, as a function of the temperature . Panel (b) shows the temperature dependence of the quantum isobaric heat capacity of MOF-5 with molecules of methane. Panel (c) shows the quantum isobaric heat capacity of the MOF with molecules of methane as a function of at 300 K. Error bars indicate statistical uncertainity.

Figure 5 also reveals that the results obtained using the simple and computationally inexpensive harmonic approximation (dotted line) are in good agreement with the exact values computed with PIMD. This implies that anharmonic quantum contributions to the heat capacity and the effect of an adequate anharmonic sampling are small for the empty MOF. Moreover, the harmonic approximation using the UFF4MOF force field 95 yields largely similar results (see SI Section S6.1), so that a less accurate inexpensive and generic model for the potential energy surface is capable of reproducing the heat capacity of the empty framework.

Another notable detail of our calculations is that the harmonic approximation was used to estimate the isochoric heat capacity instead of the isobaric one. As the isobaric and isochoric heat capacities are almost the same, the MOF behaves like a regular solid, despite its large negative thermal expansion coefficient. The harmonic approximation could therefore serve as an efficient procedure to accurately estimate the heat capacity of the empty framework in the increasing number of high-throughput MOF screenings. 96, 11, 9, 97, 98

In order to study the effect of adsorbates, we start by considering the case of a loading of 100 methane molecules per conventional unit cell (8(ZnO(CO)) 100 CH). Although a high-level PIMD strategy might not be required to estimate the heat capacity of the empty MOF host, PIMD proves to be crucial to capture the correct temperature dependence of the loaded system, as can be seen in the left panel of Figure 6. Here, anharmonic effects become important as we observe differences between the PIMD results and the harmonic approximation. The discrepancy in the qualitative behavior of the heat capacity between both techniques can be understood through the mobility of the guest molecules in the large pores of the framework, which cannot be adequately captured by a harmonic approximation. 99 These low-frequency anharmonic motions explain why we find at the classical MD level a similar low-temperature dependence as in PIMD, but only PIMD simulations include both anharmonic and nuclear quantum effects correctly. Interestingly, the combination of both effects yields a heat capacity that does not change monotonically, but exhibits a minimum at about 200 K.

Figure 7: Panels (a), (b) and (c) respectively show the decomposition of the specific heat capacity of the MOF and the adsorbate system into host, host-guest, and guest-guest contributions for different gas loadings . The curves were obtained by deriving a polynomial fit to the energy as a function of temperature.

Extending towards other loadings of methane in the middle panel of Figure 6, it becomes clear that the heat-capacity minimum as a function of temperature depends strongly on the number of guests and becomes more pronounced at higher loadings. Even when expressing the heat capacity normalized to the total mass of the system, one can see that at a fixed temperature increases almost linearly with the loading (see Figure 6 (c) at 300 K). For the volumetric heat capacity, i.e., the heat capacity per unit of volume of the system, similar results are obtained (see SI Section S6.3).

Figure 8: The decomposition of the total heat capacity of MOF-5 with 100 methane molecules per unit cell into covalent (COV), electrostatic (EI), and van der Waals (DISP) contributions. The curves were obtained by deriving a polynomial fit to the energy as a function of temperature.

To rationalize the origin of the non-monotonic temperature dependence of the heat capacity, we determine which interactions give the most substantial contribution to . To this end, the force-field energy contributions are decomposed in terms of the host, host-guest, and guest-guest interactions (see SI Section S6.2). Figure 7 displays the most important results of this analysis. The host and guest-guest contributions to the specific heat capacity are visualized in panels (a) and (c). The shape of the different host curves appears to be independent of the loading. In fact, when rescaled to the mass of the empty MOF, the curves coincide with one another and with the curve obtained within the harmonic approximation. This demonstrates that the degrees of freedom of the MOF-5 framework, which are more strongly quantized, are predominantly harmonic and do not significantly change due to the interaction with methane. Their contribution to the total heat capacity per unit mass, however, decreases with the loading due to a change in the mass balance. The guest-guest interactions, on the other hand, are relatively constant and only show a small increase when going from 100 K to 500 K, due to the activation of high-frequency vibrational modes. The most interesting contribution arises from the host-guest interactions, which explains the non-monotonic behavior of the specific heat capacity of the guest-loaded system. The contribution of these interactions decreases with a sharp temperature dependence when sufficient guest molecules are present inside the pores. The large heat capacity at temperatures lower than 100 K originates from the known first-order structural phase transition of methane in MOF-5 at 60 K, 100, 101 from which we observe the decreasing tail. Since the methane molecules are more localized at low temperatures, the attractive host-guest interactions allow to efficiently store thermal energy. At higher temperatures, from 250 K to 500 K, the host-guest contributions become negligible as the confined guests become more mobile and less bound to the framework, so that the increase in thermal energy can no longer be stored in the physical interactions between the methane guests and the MOF-5 host.

Another decomposition of the force-field energy in terms of the covalent, electrostatic, and van der Waals interactions shows that the short-range covalent interactions and thus the network of chemical bonds (Figure 8) dominate the contributions to the heat capacity. For the empty MOF-5 framework, the noncovalent interactions are negligble (see SI Section S6.2). This confirms that the heat capacity of empty MOF-5 can be approximated by considering only contributions from the separate molecular fragments of the material 12 and suggests why the harmonic approximation works well for this material. For the loaded framework, the noncovalent part starts to play a role, which is especially true for the host-guest interactions. Not surprisingly, in the case of nonpolar methane molecules, these interactions are dominated by the van der Waals terms in the force field (see SI Section S6.2). This suggests that the use of different, more polar, guests in which electrostatic interactions play a more prominent role (e.g. CO) could give rise to other interesting phenomena. However, care must be taken in interpreting these different terms as a separation is not unambiguously defined and might be force-field dependent.

4.3 The interplay of gas loading, anharmonicities, and quantum effects

Our analysis of the structural and thermal properties of methane-loaded MOF-5 shows that the total system does not always need a full treatment of anharmonicities and NQEs. This suggests that a full path integral sampling of the entire system may not be necessary, especially if qualitative trends are to be studied. Hence, inspired by our results, we propose an empirical formula for the volume and the heat capacity in which the most important effects, i.e., anharmonicities and/or NQEs, are captured and which might prove to be beneficial for future studies of guest-loaded MOFs.

As discussed above, the main difference between the volume with or without NQEs comes from zero-point effects in the lattice. The correct volume can therefore be estimated as follows:

where anh stands for the inclusion of anharmonicities with MD, and cl and qn denote the use of classical or path integral MD respectively. The left most panel of Figure 9 indicates that this approximate volume agrees very well with the exact results obtained from PIMD simulations. A more stringent test is the thermal expansion coefficient, which is – as shown in Figure 9 (b) – also in excellent agreement with the PIMD results. For systems where a first-principles treatment of the potential energy surface is required and PIMD simulations are too expensive, other approximate techniques such as the quasi-harmonic approximation or classical MD with a quantum thermostat 102 could be used to estimate the zero-point effects 83.

In contrast, we observed in the previous section that the heat capacity of the framework could be estimated with a harmonic approximation, while the guest-host interactions are dominated by anharmonicities. For that reason, we propose:

in which the high frequency modes of adsorbate and the MOF are treated in a harmonic fashion and the host-guest interactions are treated classically. For , the Dulong-Petit law can be used. As shown in the rightmost panel of Figure 9, our suggested approximation is even able to qualitatively reproduce the heat capacity minimum for a loading of methane molecules. Beyond 200 K, the agreement between the empirical expression and the exact PIMD becomes quantitatively correct. This method could thus be an inexpensive route to estimate the heat capacity of guest-loaded MOFs.

Figure 9: Panels (a), (b) and (c) respectively show the temperature dependence of the cell parameter () for MOF-5 with 50 and 150 molecules of methane, the linear thermal expansion coefficient of MOF-5 with molecules of methane as function of , and the isobaric heat capacity of MOF-5 with 100 molecules of methane, obtained with classical MD (dashed), PIMD (solid), and the approximation introduced in the work (dot-dashed). Error bars indicate statistical uncertainity.

5 Conclusions

To summarize, we developed an efficient and accurate methodology to calculate the isobaric thermophysical properties of materials, that is generally applicable and therefore ideally suited to the study of guest-loaded MOFs. For this purpose, we derived and implemented the necessary algorithms in i-PI 52 to perform simulations in the quantum isothermal-isobaric ensemble using the Suzuki-Chin path integral molecular dynamics framework. The method is rigorous and can be seamlessly combined with other cost-reduction techniques, which facilitates a huge reduction of the computational cost compared to standard techniques.

We demonstrated the applicability of our approach by investigating the heat capacity of the prototypical MOF-5 loaded with different numbers of methane molecules. We observed that the level of statistical sampling that is needed to achieve quantitative accuracy depends on the property of interest. For all the cases we considered, we found the framework to behave in a strongly quantized manner, but to be largely amenable to a harmonic treatment. The adsorbates, on the other hand, show only mild quantum effects in their intermolecular interactions, but require a full anharmonic description. The heat capacity shows a particularly subtle interplay of quantum and anharmonic fluctuations, that results in a non-monotonic temperature dependence of the heat capacity with a minimum.

Through a decomposition of the heat capacity into molecular interactions, we find that the host-guest interactions are responsible for this behavior, as their contribution to the total heat capacity decreases with temperature. By comparing the behavior of different classes of framework materials and guest molecules, this may reveal new design rules to optimize the thermal behavior of a storage material over a broad range of temperatures and levels of loading. Our approach provides an affordable route to perform benchmark studies and approximation strategies to carry out the high-throughput studies that are needed to obtain a complete understanding of the interplay between framework, adsorbate, and quantum mechanical and anharmonic fluctuations that determine the thermophysical properties of MOFs.

6 Acknowledgements

We acknowledge financial support by the Swiss National Science Foundation (project ID 200021-159896), the Fund for Scientific Research Flanders (FWO), the Research Board of Ghent University and the European Union’s Horizon 2020 research and innovation programme (consolidator ERC grant agreement No. 647755-DYNPOR (2015-2020)). The computational resources and services used in this work were provided by VSC (Flemish Supercomputer Center), funded by Ghent University, FWO, and the Flemish Government department EWI. We also want to acknowledge R. Semino, S.M.J. Rogge and L. Vanduyfhuys for their valuable input and discussions.

7 Supporting Information

Derivation of the Suzuki-Chin partition functions, forces, virials and total stress; the equations of motion and integration schemes for sampling the isothermal-isobaric scheme ensemble; derivation of the heat capacity estimators; details of the generation and validation of the first principles force field; temperature dependence of the adsorption isotherms for methane; the formula of the harmonic heat capacity, the decomposition of the total heat capacity into contributions from different physical interactions and the temperature dependence of the volumetric heat capacity.

8 Author Contributions

VK and JW contributed equally to this work.

Appendix A Suzuki-Chin forces and virials

Recently, some of the present authors showed that the high-order force can be estimated without computing the Hessian in a finite difference fashion 26 (see SI Section S1.5):


where and is a normalization factor, so that represents the root mean square displacement applied to each atom. This avoids the explicit calculation of the Hessian and allows for the direct sampling of the Suzuki-Chin canonical ensemble. Following a similar strategy (see SI Section S1.6), we show that the high-order component of the virial can be estimated as:


Appendix B Suzuki-Chin Estimators

Figure 10: The top panel shows the convergence of the isobaric heat capacity of liquid water at 300 K, modeled by the q-TIP4P/f potential, as a function of number of replicas. The blue curve was obtained from standard PIMD using the Yamamoto estimator while the red and green curves were obtained from SC PIMD using the Yamamoto (TD) and OP estimator respectively. The dashed black line represents the experimental result. Error bars indicate statistical uncertainity. The bottom panel shows that instantaneous values of the the computationally expensive term when computed with the OP and the Yamamoto (TD) estimator.

The operator (OP) and thermodynamic (TD) estimators of the energy for Suzuki-Chin isothermal-isobaric PIMD simulations of a single particle system, described in Section 2.4, are given by


where is the position of the centroid, the total Suzuki-Chin force, and denotes a thermodynamic average over the quantum isothermal-isobaric ensemble. The generalization to a many particle system is straightforward. The OP or TD estimators of the enthalpy are obtained through the simple expression:


The OP or TD (that we will refer to as the Yamamoto estimator 57) double virial estimators of the heat capacity take the following form:


where for TD is given in Ref. 50 and


The double virial force , which depends on the Hessian of the physical potential, can be estimated in a computationally efficient way by finite differences:


We test this estimator on the well known isobaric heat capacity of liquid water which is 1 or 9 per molecule at 300 K. The water is modeled with the q-TIP4P/f force field 103 and the heat capacity is estimated with SC PIMD using the aforementioned DCV estimators. As shown in Figure 10, both estimators converge to the same value and are in excellent agreement with the experiments. More importantly, however, the variance of the computationally expensive term with the OP method is almost two orders of magnitude smaller than with the Yamamoto estimator. The OP estimator also requires one fourth of the number of force evaluations than its Yamamoto counterpart and should therefore always be preferred.


  • Slater and Cooper 2015 Slater, A. G.; Cooper, A. I. Function-led design of new porous materials. Science 2015, 348, aaa8075.
  • Furukawa et al. 2013 Furukawa, H.; Cordova, K. E.; O’Keeffe, M.; Yaghi, O. M. The chemistry and applications of metal-organic frameworks. Science 2013, 341, 1230444.
  • Mason et al. 2014 Mason, J. A.; Veenstra, M.; Long, J. R. Evaluating metal-organic frameworks for natural gas storage. Chem. Sci. 2014, 5, 32–51.
  • He et al. 2017 He, Y.; Chen, F.; Li, B.; Qian, G.; Zhou, W.; Chen, B. Porous metal-organic frameworks for fuel storage. Coord. Chem. Rev. 2017, 373, 167–198.
  • de Lange et al. 2015 de Lange, M. F.; Verouden, K. J. F. M.; Vlugt, T. J. H.; Gascon, J.; Kapteijn, F. Adsorption-Driven Heat Pumps: The Potential of Metal-Organic Frameworks. Chem. Rev. 2015, 115, 12205–12250.
  • Wang et al. 2018 Wang, S.; Lee, J. S.; Wahiduzzaman, M.; Park, J.; Muschi, M.; Martineau-Corcos, C.; Tissot, A.; Cho, K. H.; Marrot, J.; Shepard, W.; Maurin, G.; Chang, J.-S.; Serre, C. A robust large-pore zirconium carboxylate metal-organic framework for energy-efficient water-sorption-driven refrigeration. Nat. Energy 2018, 3, 985–993.
  • Sumida et al. 2012 Sumida, K.; Rogow, D. L.; Mason, J. A.; McDonald, T. M.; Bloch, E. D.; Herm, Z. R.; Bae, T.-H.; Long, J. R. Carbon Dioxide Capture in Metal-Organic Frameworks. Chem. Rev. 2012, 112, 724–781.
  • Trickett et al. 2017 Trickett, C. A.; Helal, A.; Al-Maythalony, B. A.; Zamani, Z. H.; Cordova, K. E.; Yaghi, O. M. The chemistry of metal-organic frameworks for CO capture, regeneration and conversion. Nat. Rev. Mater. 2017, 2, 17045.
  • Simon et al. 2015 Simon, C. M.; Kim, J.; Gomez-Gualdron, D. A.; Camp, J. S.; Chung, Y. G.; Martin, R. L.; Mercado, R.; Deem, M. W.; Gunter, D.; Haranczyk, M.; Sholl, D. S.; Snurr, R. Q.; Smit, B. The materials genome in action: identifying the performance limits for methane storage. Energy Environ. Sci. 2015, 8, 1190–1199.
  • Mason et al. 2015 Mason, J. A.; Oktawiec, J.; Taylor, M. K.; Hudson, M. R.; Rodriguez, J.; Bachman, J. E.; Gonzalez, M. I.; Cervellino, A.; Guagliardi, A.; Brown, C. M.; Llewellyn, P. L.; Masciocchi, N.; Long, J. R. Methane storage in flexible metal-organic frameworks with intrinsic thermal management. Nature 2015, 527, 357–361.
  • Huck et al. 2014 Huck, J. M.; Lin, L.-C.; Berger, A. H.; Shahrak, M. N.; Martin, R. L.; Bhown, A. S.; Haranczyk, M.; Reuter, K.; Smit, B. Evaluating different classes of porous materials for carbon capture. Energy Environ. Sci. 2014, 7, 4132–4146.
  • Mu and Walton 2011 Mu, B.; Walton, K. S. Thermal Analysis and Heat Capacity Study of Metal-Organic Frameworks. J. Phys. Chem. C 2011, 115, 22748–22754.
  • Kloutse et al. 2015 Kloutse, F. A.; Zacharia, R.; Cossement, D.; Chahine, R. Specific heat capacities of MOF-5, Cu-BTC, Fe-BTC, MOF-177 and MIL-53(Al) over wide temperature ranges: Measurements and application of empirical group contribution method. Microporous Mesoporous Mater. 2015, 217, 1–5.
  • Babaei and Wilmer 2016 Babaei, H.; Wilmer, C. E. Mechanisms of Heat Transfer in Porous Crystals Containing Adsorbed Gases: Applications to Metal-Organic Frameworks. Phys. Rev. Lett. 2016, 116, 025902.
  • Balestra et al. 2016 Balestra, S. R. G.; Bueno-Perez, R.; Hamad, S.; Dubbeldam, D.; Ruiz-Salvador, A. R.; Calero, S. Controlling Thermal Expansion: A Metal-Organic Frameworks Route. Chem. Mater. 2016, 28, 8296–8304.
  • Auckett et al. 2018 Auckett, J. E.; Barkhordarian, A. A.; Ogilvie, S. H.; Duyker, S. G.; Chevreau, H.; Peterson, V. K.; Kepert, C. J. Continuous negative-to-positive tuning of thermal expansion achieved by controlled gas sorption in porous coordination frameworks. Nat. Commun. 2018, 9, 4873.
  • Rogge et al. 2015 Rogge, S. M. J.; Vanduyfhuys, L.; Ghysels, A.; Waroquier, M.; Verstraelen, T.; Maurin, G.; Van Speybroeck, V. A Comparison of Barostats for the Mechanical Characterization of Metal-Organic Frameworks. J. Chem. Theory Comput. 2015, 11, 5583–5597.
  • Vanduyfhuys et al. 2018 Vanduyfhuys, L.; Rogge, S. M. J.; Wieme, J.; Vandenbrande, S.; Maurin, G.; Waroquier, M.; Van Speybroeck, V. Thermodynamic insight into stimuli-responsive behaviour of soft porous crystals. Nat. Commun. 2018, 9, 204.
  • Wieme et al. 2018 Wieme, J.; Lejaeghere, K.; Kresse, G.; Van Speybroeck, V. Tuning the balance between dispersion and entropy to design temperature-responsive flexible metal-organic frameworks. Nat. Commun. 2018, 9, 4899.
  • Paesani 2012 Paesani, F. Water in metal-organic frameworks: structure and diffusion of HO in MIL-53(Cr) from quantum simulations. Mol. Simul. 2012, 38, 631–641.
  • Borges et al. 2017 Borges, D. D.; Semino, R.; Devautour-Vinot, S.; Jobic, H.; Paesani, F.; Maurin, G. Computational Exploration of the Water Concentration Dependence of the Proton Transport in the Porous UiO-66(Zr)-(COH) Metal-Organic Framework. Chem. Mater. 2017, 29, 1569–1576.
  • Parrinello and Rahman 1984 Parrinello, M.; Rahman, A. Study of an F center in molten KCl. J. Chem. Phys. 1984, 80, 860.
  • Chandler and Wolynes 1981 Chandler, D.; Wolynes, P. G. Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids. J. Chem. Phys. 1981, 74, 4078–4095.
  • Jang et al. 2001 Jang, S.; Jang, S.; Voth, G. A. Applications of higher order composite factorization schemes in imaginary time path integral simulations. J. Chem. Phys. 2001, 115, 7832–7842.
  • Ceriotti and Manolopoulos 2012 Ceriotti, M.; Manolopoulos, D. E. Efficient First-Principles Calculation of the Quantum Kinetic Energy and Momentum Distribution of Nuclei. Phys. Rev. Lett. 2012, 109, 100604.
  • Kapil et al. 2016 Kapil, V.; Behler, J.; Ceriotti, M. High order path integrals made easy. J. Chem. Phys. 2016, 145, 234103.
  • Markland and Manolopoulos 2008 Markland, T. E.; Manolopoulos, D. E. An efficient ring polymer contraction scheme for imaginary time path integral simulations. J. Chem. Phys. 2008, 129, 024105.
  • Kapil et al. 2016 Kapil, V.; VandeVondele, J.; Ceriotti, M. Accurate molecular dynamics and nuclear quantum effects at low cost by multiple steps in real and imaginary time: Using density functional theory to accelerate wavefunction methods. J. Chem. Phys. 2016, 144, 054111.
  • Poltavsky and Tkatchenko 2016 Poltavsky, I.; Tkatchenko, A. Modeling quantum nuclei with perturbed path integral molecular dynamics. Chem. Sci. 2016, 7, 1368–1372.
  • Markland and Ceriotti 2018 Markland, T. E.; Ceriotti, M. Nuclear quantum effects enter the mainstream. Nat. Rev. Chem. 2018, 2, 0109.
  • Rogge et al. 2018 Rogge, S. M. J.; Caroes, S.; Demuynck, R.; Waroquier, M.; Speybroeck, V. V.; Ghysels, A. The Importance of Cell Shape Sampling To Accurately Predict Flexibility in Metal–Organic Frameworks. J. Chem. Theory Comput. 2018, 14, 1186–1197.
  • Martyna et al. 1999 Martyna, G. J.; Hughes, A.; Tuckerman, M. E. Molecular dynamics algorithms for path integrals at constant pressure. J. Chem. Phys. 1999, 110, 3275.
  • Ceriotti et al. 2014 Ceriotti, M.; More, J.; Manolopoulos, D. E. i-PI: A Python interface for ab initio path integral molecular dynamics simulations. Comp. Phys. Commun. 2014, 185, 1019–1026.
  • Vanduyfhuys et al. 2015 Vanduyfhuys, L.; Vandenbrande, S.; Verstraelen, T.; Schmid, R.; Waroquier, M.; Van Speybroeck, V. QuickFF: A program for a quick and easy derivation of force fields for metal-organic frameworks from ab initio input. J. Comput. Chem. 2015, 36, 1015–1027.
  • Vanduyfhuys et al. 2018 Vanduyfhuys, L.; Vandenbrande, S.; Wieme, J.; Waroquier, M.; Verstraelen, T.; Van Speybroeck, V. Extension of the QuickFF force field protocol for an improved accuracy of structural, vibrational, mechanical and thermal properties of metal-organic frameworks. J. Comput. Chem. 2018, 39, 999–1011.
  • Li et al. 1999 Li, H.; Eddaoudi, M.; O’Keeffe, M.; Yaghi, O. M. Design and synthesis of an exceptionally stable and highly porous metal-organic framework. Nature 1999, 402, 276–279.
  • Eddaoudi et al. 2002 Eddaoudi, M.; Kim, J.; Rosi, N.; Vodak, D.; Wachter, J.; O’Keeffe, M.; Yaghi, O. M. Systematic Design of Pore Size and Functionality in Isoreticular MOFs and Their Application in Methane Storage. Science 2002, 295, 469–472.
  • Zhou et al. 2007 Zhou, W.; Wu, H.; Hartman, M. R.; Yildirim, T. Hydrogen and Methane Adsorption in Metal-Organic Frameworks: A High-Pressure Volumetric Study. J. Phys. Chem. C 2007, 111, 16131–16137.
  • Martyna et al. 1994 Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant pressure molecular dynamics algorithms. J. Chem. Phys. 1994, 101, 4177.
  • Rahman 1964 Rahman, A. Correlations in the Motion of Atoms in Liquid Argon. Physical Review 1964, 136, A405–A411.
  • Uhl et al. 2016 Uhl, F.; Marx, D.; Ceriotti, M. Accelerated path integral methods for atomistic simulations at ultra-low temperatures. J. Chem. Phys. 2016, 145, 054101.
  • Shiga and Shinoda 2005 Shiga, M.; Shinoda, W. Calculation of heat capacities of light and heavy water by path-integral molecular dynamics. J. Chem. Phys. 2005, 123, 134502.
  • Takahashi and Imada 1984 Takahashi, M.; Imada, M. Monte Carlo calculation of quantum systems. II. Higher order correction. J. Phys. Soc. Jap. 1984, 53, 3765–3769.
  • Chin 1997 Chin, S. A. Symplectic integrators from composite operator factorizations. Phys. Lett. A 1997, 226, 344–348.
  • 45 The standard fourth order Suzuki-Chin Splitting contains a parameter which can be used to share the high order term between odd and even beads. In this case we have set it to zero so that the high order term is evaluated only on the even beads reducing the computational cost by a factor of two.
  • Suzuki 1995 Suzuki, M. Hybrid exponential product formulas for unbounded operators with possible applications to Monte Carlo simulations. Phys. Lett. A 1995, 201, 425–428.
  • Brain 2011 Brain, G. Higher Order Propagators in Path Integral Molecular Dynamics. Ph.D. thesis, Part II Chemistry Thesis, Oxford University, 2011.
  • Jang et al. 2014 Jang, S.; Sinitskiy, A. V.; Voth, G. a. Can the ring polymer molecular dynamics method be interpreted as real time quantum dynamics? J. Chem. Phys. 2014, 140, 154103.
  • Pérez and Tuckerman 2011 Pérez, A.; Tuckerman, M. E. Improving the convergence of closed and open path integral molecular dynamics via higher order Trotter factorization schemes. J. Chem. Phys. 2011, 135, 064104.
  • Yamamoto 2005 Yamamoto, T. M. Path-integral virial estimator based on the scaling of fluctuation coordinates: Application to quantum clusters with fourth-order propagators. J. Chem. Phys. 2005, 123, 104101.
  • Ceriotti et al. 2011 Ceriotti, M.; Brain, G. A. R.; Riordan, O.; Manolopoulos, D. E. The inefficiency of re-weighted sampling and the curse of system size in high order path integration. Proc. Royal Soc. A 2011, 468, 2–17.
  • Kapil et al. 2018 Kapil, V. et al. i-PI 2.0: A universal force engine for advanced molecular simulations. Comput. Phys. Commun. 2018, in press.
  • Tuckerman et al. 1992 Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible multiple time scale molecular dynamics. J. Chem. Phys. 1992, 97, 1990.
  • Leimkuhler and Matthews 2013 Leimkuhler, B.; Matthews, C. Robust and efficient configurational molecular sampling via Langevin dynamics. J. Chem. Phys. 2013, 138.
  • Kapil et al. 2018 Kapil, V.; Cuzzocrea, A.; Ceriotti, M. Anisotropy of the Proton Momentum Distribution in Water. J. Phys. Chem. B 2018, 122, 6048–6054.
  • Poltavsky et al. 2018 Poltavsky, I.; Jr., R. A. D.; Tkatchenko, A. Perturbed path integrals in imaginary time: Efficiently modeling nuclear quantum effects in molecules and materials. J. Chem. Phys. 2018, 148, 102325.
  • Yamamoto 2005 Yamamoto, T. M. Path-integral virial estimator based on the scaling of fluctuation coordinates: Application to quantum clusters with fourth-order propagators. J. Chem. Phys. 2005, 123, 104101.
  • Tafipolsky et al. 2007 Tafipolsky, M.; Amirjalayer, S.; Schmid, R. Ab initio parametrized MM3 force field for the metal-organic framework MOF-5. J. Comput. Chem. 2007, 28, 1169–1176.
  • Frisch et al. 2016 Frisch, M. J. et al. Gaussian16 Revision B.01. 2016; Gaussian Inc. Wallingford CT.
  • Becke 1988 Becke, A. D. Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior. Phys. Rev. A 1988, 38, 3098–3100.
  • Becke 1993 Becke, A. D. Density-Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648–5652.
  • Lee et al. 1988 Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-Energy Formula into a Functional of the Electron Density. Phys. Rev. B 1988, 37, 785–789.
  • Krishnan et al. 1980 Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. 20. Basis Set for Correlated wave-Functions. J. Chem. Phys. 1980, 72, 650–654.
  • Hay and Wadt 1985 Hay, P. J.; Wadt, W. R. Ab initio effective core potentials for molecular calculations. Potentials for the transition metal atoms Sc to Hg. J. Chem. Phys. 1985, 82, 270.
  • Verstraelen et al. 2016 Verstraelen, T.; Vandenbrande, S.; Heidar-Zadeh, F.; Vanduyfhuys, L.; Van Speybroeck, V.; Waroquier, M.; Ayers, P. W. Minimal Basis Iterative Stockholder: Atoms in Molecules for Force-Field Development. J. Chem. Theory Comput. 2016, 12, 3894–3912.
  • Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865.
  • Mortensen et al. 2005 Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-space grid implementation of the projector augmented wave method. Phys. Rev. B 2005, 71, 035109.
  • Enkovaara et al. 2010 Enkovaara, J. et al. Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. J. Phys.: Condens. Matter 2010, 22, 253202.
  • Verstraelen et al. 2017 Verstraelen, T.; Tecmer, P.; Heidar-Zadeh, F.; González-Espinoza, C. E.; Kim, T. D.; Boguslawski, K.; Fias, S.; Vandenbrande, S.; Berrocal, D.; Ayers, P. W. Horton 2.1.0. 2017;
  • Lii and Allinger 1989 Lii, J. H.; Allinger, N. L. Molecular Mechanics. The MM3 Force Field for Hydrocarbons. 3. The van der Waals Potentials and Crystal Data for Aliphatic and Aromatic Hydrocarbons. J. Am. Chem. Soc. 1989, 111, 8576–8582.
  • Allinger et al. 1994 Allinger, N. L.; Zhou, X.; Bergsma, J. Molecular Mechanics Parameters. J. Mol. Struc.-THEOCHEM 1994, 312, 69–83.
  • Sun 1998 Sun, H. COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase Applications - Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338–7364.
  • Dubbeldam et al. 2016 Dubbeldam, D.; Calero, S.; Ellis, D. E.; Snurr, R. Q. RASPA: molecular simulation software for adsorption and diffusion in flexible nanoporous materials. Mol. Simul. 2016, 42, 81–101.
  • Plimpton 1995 Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1–19.
  • Nosé 1984 Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. Mol. Phys. 1984, 52, 255–268.
  • Hoover 1985 Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695–1697.
  • Martyna et al. 1992 Martyna, G. J.; Klein, M. L.; Tuckerman, M. Nose-Hoover Chains - the Canonical Ensemble Via Continuous Dynamics. J. Chem. Phys 1992, 97, 2635–2643.
  • Martyna et al. 1994 Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177–4189.
  • Martyna et al. 1996 Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit reversible integrators for extended systems dynamics. Mol. Phys. 1996, 87, 1117–1157.
  • Ceriotti et al. 2010 Ceriotti, M.; Parrinello, M.; Markland, T. E.; Manolopoulos, D. E. Efficient stochastic thermostatting of path integral molecular dynamics. J. Chem. Phys. 2010, 133, 124104.
  • Bussi and Parrinello 2007 Bussi, G.; Parrinello, M. Accurate sampling using Langevin dynamics. Phys. Rev. E 2007, 75, 56707.
  • Bussi et al. 2009 Bussi, G.; Zykova-Timan, T.; Parrinello, M. Isothermal-isobaric molecular dynamics using stochastic velocity rescaling. J. Chem. Phys. 2009, 130, 074101.
  • Lamaire et al. submitted Lamaire, A.; Wieme, J.; Rogge, S. M. J.; Waroquier, M.; Van Speybroeck, V. On the importance of anharmonicities and nuclear quantum effects in modelling the structural properties and thermal expansion in MOF-5. submitted,
  • Balog et al. 2000 Balog, E.; Hughes, A. L.; Martyna, G. J. Constant pressure path integral molecular dynamics studies of quantum effects in the liquid state properties of n-alkanes. J. Chem. Phys. 2000, 112, 870–880.
  • Pereyaslavets et al. 2018 Pereyaslavets, L.; Kurnikov, I.; Kamath, G.; Butin, O.; Illarionov, A.; Leontyev, I.; Olevanov, M.; Levitt, M.; Kornberg, R. D.; Fain, B. On the importance of accounting for nuclear quantum effects in ab initio calibrated force fields in biological simulations. Proc. Natl. Acad. Sci. U.S.A. 2018, 115, 8878–8882.
  • Veit et al. 2018 Veit, M.; Jain, S. K.; Bonakala, S.; Rudra, I.; Hohl, D.; Csányi, G. Equation of state of fluid methane from first principles with machine learning potentials. arXiv preprint arXiv:1810.10475v1 2018,
  • Ravikovitch and Neimark 2006 Ravikovitch, P. I.; Neimark, A. V. Density Functional Theory Model of Adsorption Deformation. Langmuir 2006, 22, 10864–10868.
  • Joo et al. 2013 Joo, J.; Kim, H.; Han, S. S. Volume shrinkage of a metal-organic framework host induced by the dispersive attraction of guest gas molecules. Phys. Chem. Chem. Phys. 2013, 15, 18822–18826.
  • Zhou et al. 2008 Zhou, W.; Wu, H.; Yildirim, T.; Simpson, J. R.; Hight Walker, A. R. Origin of the exceptional negative thermal expansion in metal-organic framework-5 ZnO(1,4-benzenedicarboxylate). Phys. Rev. B 2008, 78, 054114.
  • Lock et al. 2010 Lock, N.; Wu, Y.; Christensen, M.; Cameron, L. J.; Peterson, V. K.; Bridgeman, A. J.; Kepert, C. J.; Iversen, B. B. Elucidating Negative Thermal Expansion in MOF-5. J. Phys. Chem. C 2010, 114, 16181–16186.
  • Lock et al. 2013 Lock, N.; Christensen, M.; Wu, Y.; Peterson, V. K.; Thomsen, M. K.; Plitz, R. O.; Ramirez-Cuesta, A. J.; McIntyre, G. J.; Norén, K.; Kutteh, R.; Kepert, C. J.; Kearley, G. J.; Iversen, B. B. Scrutinizing negative thermal expansion in MOF-5 by scattering techniques and ab initio calculations. Dalton Trans. 2013, 42, 1996–2007.
  • 92 A 10 increase in volume due to NQEs, assuming the molecules to be spherical is associated with a times increase in the effective radius. The classical ”radius” is close to 4 Å. And thus the increase in the inter-molecular distance due to NQEs coming from the isotropic expansion of the gas is Å.
  • Ming et al. 2014 Ming, Y.; Purewal, J.; Liu, D.; Sudik, A.; Xu, C.; Yang, J.; Veenstra, M.; Rhodes, K.; Soltis, R.; Warner, J.; Gaab, M.; Müller, U.; Siegel, D. J. Thermophysical properties of MOF-5 powders. Microporous Mesoporous Mater. 2014, 185, 235–244.
  • Liu et al. 2012 Liu, D.; Purewal, J. J.; Yang, J.; Sudik, A.; Maurer, S.; Müller, U.; Ni, J.; Siegel, D. J. MOF-5 composites exhibiting improved thermal conductivity. Int. J. Hydrogen Ener. 2012, 37, 6109–6117.
  • Addicoat et al. 2014 Addicoat, M. A.; Vankova, N.; Akter, I. F.; Heine, T. Extension of the Universal Force Field to Metal-Organic Frameworks. J. Chem. Theory Comput. 2014, 10, 880–891.
  • Wilmer et al. 2012 Wilmer, C. E.; Leaf, M.; Lee, C. Y.; Farha, O. K.; Hauser, B. G.; Hupp, J. T.; Snurr, R. Q. Large-scale screening of hypothetical metal-organic frameworks. Nat. Chem. 2012, 4, 83–89.
  • Moghadam et al. 2018 Moghadam, P. Z.; Islamoglu, T.; Goswami, S.; Exley, J.; Fantham, M.; Kaminski, C. F.; Snurr, R. Q.; Farha, O. K.; Fairen-Jimenez, D. Computer-aided discovery of a metal-organic framework with superior oxygen uptake. Nat. Commun. 2018, 9, 1378.
  • Tabor et al. 2018 Tabor, D. P.; Roch, L. M.; Saikin, S. K.; Kreisbeck, C.; Sheberla, D.; Montoya, J. H.; Dwaraknath, S.; Aykol, M.; Ortiz, C.; Tribukait, H.; Amador-Bedolla, C.; Brabec, C. J.; Maruyama, B.; Persson, K. A.; A., A.-G. Accelerating the discovery of materials for clean energy in the era of smart automation. Nat. Rev. Mater. 2018, 3, 5–20.
  • Tsivion and Head-Gordon 2017 Tsivion, E.; Head-Gordon, M. Methane Storage: Molecular Mechanisms Underlying Room-Temperature Adsorption in ZnO(BDC) (MOF-5). J. Phys. Chem. C 2017, 121, 12091–12100.
  • Wu et al. 2009 Wu, H.; Zhou, W.; Yildirim, T. Methane Sorption in Nanoporous Metal-Organic Frameworks and First-Order Phase Transition of Confined Methane. J. Phys. Chem. C 2009, 113, 3029–3035.
  • Kuchta et al. 2017 Kuchta, B.; Dundar, E.; Formalik, F.; Llewellyn, P. L.; Firley, L. Adsorption-Induced Structural Phase Transformation in Nanopores. Angew. Chem. Int. Ed. 2017, 56, 16243–16246.
  • Ceriotti et al. 2009 Ceriotti, M.; Bussi, G.; Parrinello, M. Langevin Equation with Colored Noise for Constant-Temperature Molecular Dynamics Simulations. Phys. Rev. Lett. 2009, 102, 020601.
  • Habershon et al. 2009 Habershon, S.; Markland, T. E.; Manolopoulos, D. E. Competing quantum effects in the dynamics of a flexible water model. The Journal of Chemical Physics 2009, 131, 024501.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description