Sorption-Induced Breathing in the Flexible Metal Organic Framework

Dec 17, 2008 - Density functional theory (DFT) and force-field-based calculations have been carried out on the breathing metal−organic framework MIL...
1 downloads 0 Views 1MB Size
Download PDF
544

J. Phys. Chem. C 2009, 113, 544–552

Sorption-Induced Breathing in the Flexible Metal Organic Framework CrMIL-53: Force-Field Simulations and Electronic Structure Analysis David S. Coombes, Furio Cora`, Caroline Mellot-Draznieks, and Robert G. Bell* Department of Chemistry, UniVersity College London, 20 Gordon Street, London WC1H 0AJ, United Kingdom ReceiVed: July 28, 2008

Density functional theory (DFT) and force-field-based calculations have been carried out on the breathing metal-organic framework MIL-53(Cr) in both its large- and narrow-pore forms. In its sorbate-free form, the large-pore structure appears to be the global minimum. We develop a hybrid force field combining ionicmodel potentials, used for modeling inorganic solids, with molecular mechanics terms for the organic part. This gives an energy difference of close to 30 kJ mol-1 between the large- and narrow-pore forms. Calculations in which water molecules are introduced into the structures illustrate how the energetics of physisorption are able to drive the pore breathing process, with the water molecules, being more strongly stabilized in the narrow-pore form, favoring pore closure at loadings of more than one molecule per unit cell. 1. Introduction Recent years have seen an effusion of work on nanoporous metal-organic frameworks (MOFs), also often known as coordination networks. The defining feature of these crystalline materials is that they possess metal centers linked through multidentate organic ligands, forming bonded three-dimensional networks. The inorganic part of the network may consist of isolated metal atoms or of metal-oxide moieties, themselves forming either one-dimensional clusters or two-dimensional chains; the ligands commonly bond via oxygen (as in the carboxylate function), or nitrogen donor atoms. The huge numbers of potential ligands and network topologies lead to the possibility of tailor-made voids of molecular size. Furthermore, the “spacer” ligands may be functionalized, leading to a preordained chemical reactivity. Many applications have thus been proposed for MOF materials, including gas storage and separation, catalysis, and sensing. A number of reviews have appeared in which the structure and chemistry of MOFs are discussed.1-3 An intriguing property of many MOFs, generally not found in purely inorganic nanoporous materials such as zeolites, is that of framework flexibility. This can lead, for example, to a “gate-opening” phenomenon4 where a guest molecule is able to pass (into a larger cavity) through a pore that would ostensibly seem too narrow when measured from the static crystal structure. The incorporation of flexible ligands in order to introduce deliberately the property of framework flexibility has also been reported.5 Another manifestation of framework flexibility is that of pore “breathing”, where the whole structure can expand to admit guest molecules, with an analogous shrinking of the overall pore structure under desorption conditions. A unique example of a MOF material exhibiting such pore breathing is that of MIL-53, particularly the chromiumcontaining version6,7 with which we shall be concerned here. This structure is part of a series of isotypical MOFs synthesized by the group of Fe´rey.6-10 The structure of the chromiumcontaining MIL-53 [chemical formula Cr(OH)(O2C-C6H4-CO2] consists of infinite chains of corner-sharing CrO4(OH)2 octa* To whom correspondence should be addressed: e-mail r.g.bell@ ucl.ac.uk.

hedral units interconnected by p-benzenedicarboxylate ligands. The result is a three-dimensional metal-organic framework containing monodirectional diamond-shaped channels. This structure has been shown to exhibit unusual framework flexibility (breathing) upon adsorption/desorption of various sorbates, such as water11 and CO2,12 and the coadsorption of water and CO2.13 The more open, large-pore form (to which we refer as MIL-53lp, consistent with the nomenclature of Ramsahye et al.14) is obtained upon calcination of the as-synthesized compound, while the adsorption of water molecules at room temperature leads to the closed, narrow-pore form (MIL-53np). Upon heating, the water adsorption is reversed, leading to recovery of the open form, MIL-53lp. The difference between the open and closed forms is striking: although the bonded topology of the framework is unchanged, some atomic positions move by more than 5 Å, and the unit cell volume decreases by over 30% upon going from open to closed forms. Clearly a subtle balance of forces is at play in the MIL-53 system, which causes the pore breathing to be triggered by the adsorption or desorption of a relatively small number of guest molecules. A previous study15 has explored the complex nature of the interactions between the framework and guest species, as well as illuminating many structural details concerned with the breathing process. Molecular simulations clearly have a role to play in the study of flexible MOFs, considering their structural complexity and the wide variety of sorption response behavior which is observed. Although it does not exhibit pore breathing to the same degree as MIL-53, MOF-5 has been widely studied. Greathouse and Allendorf16 were able to simulate, using molecular dynamics, the variation of cell parameters, and eventual structural collapse, upon adsorption of water. More recently they have reproduced the negative thermal expansion of this material as well as its response to the adsorption of various organic guest molecules.17 The force field used in both cases, was based on cVff18 bonded terms for the organic linker with nonbonded Lennard-Jones potentials representing the inorganic Zn-O interactions. Schmid and co-workers,19,20 using an MM3-based force field, demonstrated the importance of a flexible framework model when simulating the diffusion of guest molecules such as benzene in MOF-5.

10.1021/jp809408x CCC: $40.75  2009 American Chemical Society Published on Web 12/17/2008

Sorption-Induced Breathing in Metal-Organic Framework In order to more fully understand and describe the pore breathing effect in MIL-53 under a variety of sorption conditions, one approach would be to carry out fully flexible molecular simulations on the system. In a parallel simulation study,21 this technique has been used to examine pore breathing as a response to CO2 adsorption as well as its thermal activation. In addition, computational models employing static calculations could conceivably tell us a great deal about the role of the different inter- and intramolecular forces in the crystal chemistry of MIL-53, and specifically in the breathing process. We investigate this topic with density functional theory calculations on the sorbate-free crystal, in an attempt to sample the energy surface of the pore breathing process. We then develop a force field that is able to describe the structure and thermodynamics of the Cr-MIL-53 system with reasonable accuracy. This force field is further applied to investigate framework breathing upon adsorption of water, demonstrating how pore breathing may be thermodynamically driven by adsorption. 2. Methodology The main goal of our work is to understand the forces that enable large structural relaxations in metal-organic framework materials and provide an atomic-level description of the process in the MIL-53 class of materials (exemplified here by the Cr3+ version). The availability of a force field that would describe these materials is an essential part of this process. Our strategy involved the combination of components from two force fieldssone for the inorganic subunit and one for the organic subunitswhich are both used extensively for modeling inorganic and organic materials, respectively. The strategy differs from others that use charges derived from electronic structure calculations, in that we use a combination of formal charges on the inorganic subunit and partial charges, assigned from bond/ valence considerations, on the organic subunits. This choice has the main advantage that the force field should be more easily transferable to other MOF materials. DFT calculations have been performed to characterize the bonding in Cr-MIL-53 and justify the choice of force field. The quantum mechanical results have also been used to sample the potential energy surface of the material, against which the force-field results can be compared. In the following sections we describe first the DFT calculations, followed by our approach to developing the force field. 2.1. Periodic DFT Calculations. The Cr-containing MIL53 material discussed in the paper was investigated with periodic density functional theory (DFT) calculations, employing the PBE22 and hybrid-exchange PBE023 functionals and localized orbitals, as implemented in the CRYSTAL0624 code. The use of orbital-dependent functionals, obtained here via hybrid exchange, is of particular interest in the field of MOFs as they share the main electronic feature of strongly correlated materials such as NiO and superconducting cuprates in that they contain transition metal ions with localized d electrons. These compounds are often poorly described in DFT by standard (LDA or GGA) (local) exchange and correlation functionals, and they require orbital-dependent terms to yield the correct ground electronic state. The basis sets chosen are expressed in a series of Gaussiantype functions and correspond to a 86411-41G contraction for Cr25 and a split valence plus polarization 6-311G** for H, C, and O.26 All elements are treated at the all-electron level. We used a Pack-Monkhorst grid of 4 × 4 × 4 and truncation thresholds of (7 7 7 7 14) for the Coulomb and exchange integrals, while the SCF convergence threshold was set to 1 × 10-7 hartree for both eigenvalues and total energy. These tolerances ensure sufficient numerical accuracy in the calculations.

J. Phys. Chem. C, Vol. 113, No. 2, 2009 545 Geometry optimizations were checked against the root-meansquare (rms) and absolute value of the largest component for both gradient and nuclear displacements, using the values of 1.2 × 10-4 atomic unit for the maximum gradient and 1.8 × 10-4 atomic unit for the maximum displacement. Optimizations were iterated until the lattice parameters showed no further change upon restart. The starting coordinates for our geometry optimizations were taken from the experimentally derived geometries for the large- and narrow-pore phases. Both structures have been fully optimized in the experimentally determined space groups Imma for the open form and C2/c for the closed form. In addition, we also performed a number of partial optimizations at fixed lattice parameters taken both from experiment and force-field optimizations, to provide a more detailed sampling of the potential energy surface of the solid. Finally, full geometry optimizations were also performed after each constrained fixed-lattice optimization, to ensure that the structure obtained computationally does not depend on the starting geometry, as the probability of locating local minima is high in such a complex potential energy surface as that of the flexible MIL-53 framework. Due to the need to use symmetry at this level of theory for reasons of computational efficiency, it was not possible to introduce sorbate molecules into the DFT calculations. The required drop in symmetry to P1 would make such calculations impractical. 2.2. Derivation of Hybrid Force Field. Development of a new force field for modeling MIL-53 required the ability to describe the organic ligand and its possible flexibility together with an accurate representation of the metal coordination. To describe intramolecular (bond stretch, angle bend, and torsion) and intermolecular interactions within the organic ligand, we adopted the parametrization of the well-known cVff force field,18 for its robustness and its availability with a large variety of functionalized organic molecules. Furthermore, the cVff force field is well-suited to be combined with the type of ionic model often used for inorganic oxide materials and has been successfully used in the past by a number of groups 27-29 investigating adsorption phenomena in zeolitic and layered materials. Bond lengths and bond angles are represented through harmonic bond stretching and angle bending potentials, plus a torsional term (functional forms are given in Table 1). For the Cr(III) coordination sphere, an ionic model based on that of Cormack et al.30 was adopted. The overall charge on chromium was +3.0. For each Cr(III) center, two distinct types of Cr-O linkage are present in MIL-53: four “equatorial” bonds to oxygens of the carboxylate ligands (denoted here as O1), and two along the chain direction (i.e. axial) to the bridging hydroxyl oxygen (O2). For the Cr-O2 potential, Buckingham parameters were modified slightly from those of Cormack et al.30 by fitting to Cr-O2 axial bond lengths in order to reproduce the observed cell parameter along the Cr-O2-Cr chain direction. The overall charge on the bridging hydroxyl group was -1.0. The Cr-O1 (carboxylate) interaction was also modeled as a Buckingham potential, but the parameters were refitted to a greater extent to account for the fact that lower cVff partial charges were assigned to the organic ligand, the charge on O1 being -0.57. All parameters are given in Table 1. On the chromium ions, a core-shell model was used to account for possible polarization of this transition metal caused by ligand field effects:30 the hydroxyl group was modeled by use of a Morse term and the parametrization of Saul et al.31 The other interactions between the inorganic subnetwork and the organic ligands were modeled as follows: a four-body term

546 J. Phys. Chem. C, Vol. 113, No. 2, 2009

Coombes et al.

TABLE 1: Force-Field Parameters Inorganic Subunit

()

EBuckingham ) Aij exp

Cij rij - 6 Fij rij

1 Ecore-shell ) kc-src-s2 2 EMorse ) De{(1 - exp[-R(r - r0)])2 - 1} Buckingham Buckingham shell model Morse

A ) 1150 eV A ) 22 764 eV Y ) 0.97 e De ) 7.0525 eV

O2-Cr O2-O2 Cr H2-O2

F ) 0.307 035 Å F ) 0.149 Å kc-s ) 67 eV Å-2 R ) 2.1986 Å-1

C ) 0 eV Å6 C ) 27.88eV Å6 re ) 0.9485 Å

fitted from ref 30 ref 35 ref 30 ref 31

Organic Liganda

EijLennard-Jones )

A B r12 r6

1 1 Eintramolecular ) kij(rij - r0)2 + kijk(θijk - θ0)2 + kφ[1 + s cos (nφ - φ0)] 2 2 Inorganic-Organic Interactions Buckingham Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones Lennard-Jones three body torsion a

O1-Cr C1-Cr C1-O2 H2-C1 C2-Cr C2-O2 H2-C2 H1-Cr H1-O2 H1-H2 O1-O2 H2-O1 C2-O1-Cr C1-C2-O1-Cr

A ) 2811.073 32 eV A ) 82 612.877 29 eV Å12 A ) 39 031.73351 eV Å12 A ) 0 eV Å12 A ) 82 612.877 29 eV Å12 A ) 39 031.733 51 eV Å12 A ) 0 eV Å12 A ) 4042.486 782 eV Å12 A ) 1909.935 496 eV Å12 A ) 0 eV Å12 A ) 11 833.9137 eV Å12 A ) 0 eV Å12 kijk ) 1.2 eV rad-2 kΦ ) 0.23 eV

F ) 0.238 721 Å B ) 136.951 275 5 eV Å6 B ) 35.26589306 eV Å6 B ) 0 eV Å6 B ) 136.951 275 5 eV Å6 B ) 35.265 893 06 eV Å6 B ) 0 eV Å6 B ) 21.564 727 28 eV Å6 B ) 5.553 065 229 eV Å6 B ) 0 eV Å6 B ) 21.6336 eV Å6 B ) 0 eV Å6 θ0 ) 120° n)2

C ) 0 eV Å6

Φ0 ) 180°

fitted from ref 30 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 cVff 18 this work this work

18

Harmonic, three body, torsion, and Lennard-Jones parameters were obtained by cVff.

was introduced to describe the C1-C2-O1-Cr torsion, as well as a three-body term for the C2-O1-Cr angle. The three-body potential was adapted from the cVff term acting about a protonated carboxyl oxygen. The torsional term, which influences how the carboxyl group, acting as a hinge about its oxygen atoms, can pivot to allow opening and closing of the pore, was fitted so as to give the closest match to the experimental CrO6 geometry of the open MIL53ht form. The implications of the value of this term are discussed later. All other nonbonded interactions between the inorganic subunit and organic ligand were represented by Lennard-Jones potentials: The force-field atom typing and charges used are shown in Figure 1, and the corresponding potential parameters are shown in Table 1. For water, the standard cVff three point charge model was used (charges on O and H were -0.82 and +0.41, respectively), with Lennard-Jones terms acting between water oxygen and the framework atoms. No Lennard-Jones term was used between water and chromium, since water is not anticipated to enter the inner coordination sphere of Cr. The electrostatic interactions

were calculated with an Ewald summation, with the cutoff parameters for the real and reciprocal space terms chosen to minimize the total number of terms in the two summations. 3. Results and Discussion 3.1. Periodic Density Functional Theory Calculations. Full geometry optimizations were performed on the sorbate-free framework of CrMIL-53, in the space groups determined experimentally7 for the closed low-temperature (lt) (C2/c) and open high-temperature (ht) (Imma) framework structures. Initial calculations assumed a ferromagnetically ordered structure, as it retains the full space-group symmetry and minimizes the computational cost. Results obtained with either PBE or PBE0 functional indicate that only one stable structure (one minimum in the potential energy surface) exists for the sorbate-free CrMIL-53 framework. This corresponds to the open structure labeled as lp. Geometry optimizations starting from the closed structure (np) relax gradually until the open framework structure is recovered (although described in the different C2/c space group). It is,

Sorption-Induced Breathing in Metal-Organic Framework

J. Phys. Chem. C, Vol. 113, No. 2, 2009 547

Figure 1. Force-field parametrization for MIL-53: (a) atom types; (b) charges. (Atom types refer to those used in Table 1.)

TABLE 2: Comparison of Experimental7 and Optimized Structures of CrMIL-53lpa exptl a, Å 16.733 b, Å 13.038 c, Å 6.812 volume, Å3 1486.1 R, deg 90 β, deg 90 γ, deg 90 Cr-O1, Å 1.977 Cr-O2, Å 1.896 C2-O1, Å 1.363 C2-O1-Cr, deg 151.6 a

calcd GULP/ calcd DFT/ calcd DFT/ potentials PBE PBE0 16.281 14.073 6.794 1556.6 90 90 90 1.936 1.894 1.249 138.6

16.893 13.196 6.965 1552.7 90 90 90 1.989 1.955 1.275 134.0

16.770 13.141 6.932 1527.6 90 90 90 1.981 1.943 1.258

Atom types are as in Table 1 and Figure 1.

however, important to note that the potential energy surface is very flat, and all calculations performed at fixed lattice parameters fall within 100 kJ/mol of the global minimum, despite sampling a change of more than 50% in the unit cell volume. The geometry optimized with periodic DFT calculations is reported in Table 2, where it is compared to both experimental and force-field results. Lattice parameters are slightly overestimated at the PBE level, by approximately 1% as is usually the case for PBE calculations on crystalline solids. The agreement with experiment is improved when PBE0 results are considered. Internal coordinates differ considerably from the published structure; in particular, the benzene rings of the organic ligands are bent in the experimental structural refinement, while these aromatic rings are planar in the calculated geometry. We shall discuss the implications of this structural feature later in this section, when we examine magnetic order. In addition to the PBE and PBE0 calculations detailed above, we attempted to study the CrMIL-53 structure with hybrid functionals containing a higher percentage (above 50%) of exact Hartree-Fock (HF) exchange. Such a choice has been shown to improve the reproduction of structural and electronic properties of inorganic transition metal compounds.32 We found, however, that increasing the HF exchange has a detrimental effect in the description of the aromatic organic ligands. Metal-organic frameworks therefore have different require-

ments from inorganic solids in terms of the best-performing functional, and a fraction of HF around 25% (as in the PBE0 functional) appears optimal for frameworks with aromatic bridging ligands. Upon analysis of the electronic density, we find that bonding between the transition metal ion and the organic ligands is mostly ionic in nature, with very little charge transfer due to Cr-carboxylate covalence. A suitable measure of the electronic localization is provided by the spin density on the Cr3+ ions, which is calculated as 2.918 |e| with PBE and 2.999 |e| with PBE0 by a Mulliken analysis. The spin transferred (by covalent backbonding) to the ligands is less than 0.04 |e| for both the benzenedicarboxylate and hydroxide units at the PBE level. The choice of using metal and organic subunits described with formal charges in the force field is therefore entirely justified by our electronic structure calculations. Let us now consider magnetic order in the structure. The MIL53 framework contains only one type of magnetic ion, Cr3+ in our case, and two different types of linkages between adjacent magnetic ions provided by the benzenedicarboxylate and hydroxide ligands. Each Cr3+ ion is surrounded by six nearestneighbor Cr3+ ions, of which four are linked via benzenedicarboxylate units and two via hydroxide ions. Broken-symmetry calculations have been performed on the geometry-optimized structures by use of both PBE and PBE0 functionals. The ground-state magnetic order was found to be antiferromagnetic across both benzenedicarboxylate and hydroxide ligands and at both levels of theory. We found, however, that the coupling constants between adjacent Cr ions are much higher for interactions across the benzenedicarboxylate ligands than across the hydroxide ones. With the PBE0 functional, the difference in energy (proportional to the magnetic coupling constants) between ferromagnetic and antiferromagnetic ordered structures is 2.4 kJ mol–1 across the benzenedicarboxylate ligand and only 0.07 kJ mol–1 across the hydroxide ones. The corresponding values calculated with the PBE functional are respectively 7.0 and 0.13 kJ mol–1; as expected, these values are higher than when the hybrid exchange functional PBE0 is used because the covalent Cr-ligand interactions responsible for superexchange are overestimated by local DFT. The magnetic interaction between Cr3+ ions bridged by a benzenedicarboxylate ligand is of superexchange type; it is

548 J. Phys. Chem. C, Vol. 113, No. 2, 2009 TABLE 3: Comparison of Experimental7 and Optimized Structures of CrMIL-53npa a, Å b, Å c, Å volume, Å3 R, deg β, deg γ, deg Cr-O1, Å Cr-O2, Å C2-O1, Å C2-O1-Cr, deg a

exptl

calcd (GULP/force field)

19.685 7.849 6.782 1012.6 90 104.9 90 2.077, 1.971 1.913 1.239, 1.304 135.1, 141.8

19.975 6.730 6.753 881.6 90.03 103.8 89.97 1.956, 1.942 1.903 1.250, 1.247 128.5, 132.6

Atom types are as in Table 1 and Figure 1.

mediated by the aromatic π-electron system of the organic ligand. Planarity of the benzenedicarboxylate molecule is essential for this superexchange interaction to achieve maximum efficiency. The open form of CrMIL-53 has been found experimentally to be an antiferromagnetic insulator, with a Neel temperature of 65 K;7 this value is consistent with the results of our calculations and validates the geometry optimized structure with planar benzenedicarboxylate ligands. Full geometry optimizations performed with the ground-state antiferromagnetic order show that the CrMIL-53 framework has negligible magnetostriction; the largest variation of the lattice parameters upon a change of magnetic order is less than 0.01 Å, indicating that magnetic interactions are decoupled from the large structural relaxation responsible for framework breathing. This result has two practical consequences: first, at the quantum chemical level, calculations with the less computationally demanding ferromagnetic structure will be able to reproduce the driving force for framework breathing; and second, calculations based on force fields and interatomic potentials that neglect magnetism altogether can also be employed reliably to study the large framework relaxations induced by molecular adsorption processes. As we shall see below, the energy difference between magnetic structures, being at most 7.0 kJ mol–1 at the PBE level, is smaller than the energy difference between open and closed forms of the framework, and thus magnetism is unlikely to drive any structural phase transition in CrMIL-53. It may, however, be possible that other framework topologies and compositions designed to maximize magnetic coupling may be able to couple magnetism with breathing in multifunctional applications. A final comment should be made here on the choice of functional to be used within DFT for calculations on CrMIL53. We have found that PBE and PBE0 functionals yield different quantitative values of magnetic coupling, in line with results observed earlier for purely inorganic solids; however, qualitative results are consistent at both levels of theory. Calculations employing local functionals are therefore able to reproduce correctly, at least qualitatively, all features of the electronic ground state of CrMIL-53. All results discussed following this section refer therefore only to PBE results. Having ascertained that the sorbate free framework displays only one stable structure, corresponding to the open (lp) form, we employed the periodic PBE calculations and ferromagnetic order to sample the relative energy of structures at different lattice parameters. These included the values obtained experimentally for the open and closed forms (reported in Tables 2 and 3), and the closed form optimized with our newly derived force field (Table 3). Constrained optimizations of the fractional parameters at fixed lattice parameters were performed in these

Coombes et al. cases. The energies relative to the fully optimized open form are 27.02, 70.04, and 100.15 kJ mol-1 per unit cell [i.e., 6.75, 17.51, and 25.04 kJ mol-1 per formula unit, Cr(OH)(O2C-C6H4CO2)] for the open form at the experimental lattice parameters, the closed form at the experimental lattice parameters (in the experiment, water is additionally present), and the closed form at the force-field-optimized lattice parameters (see below). These values sample the potential energy surface of CrMIL-53 around the two experimentally observed phases and provide a reference against which the force-field results can be validated. The calculations indicate that the large-pore form is the energy minimum, with the narrow-pore form being metastable. Large changes in equilibrium volume are observed (the phases examined range from 881.6 to 1552.7 Å3). Despite this, however, the two values of 70.04 and 100.15 kJ mol-1 per unit cell between the open and closed forms (experimental and forcefield-minimized, respectively) seem quite high when one considers the ease with which the structures may be interconverted upon the adsorption of a small number of guest molecules. Because we could find no minimum for the narrowpore form and merely sampled two sets of lattice parameters, the value of 70 kJ mol-1 should therefore be taken very much as an upper bound for the energy difference. We note that the difference between the structures with the experimental lp and np lattice parameters is 42 kJ mol-1, which intuitively seems a more realistic figure. 3.2. Force-Field Validation. The GULP program33 was used to obtain the minimized structures starting from both the MIL53ht and MIL-53lt experimental structures, without the presence of water or any sorbate molecules. A key feature of the force field is that, under constant pressure conditions in which the cell parameters are allowed to relax in addition to the atomic coordinates, two minima are found, depending on whether the open or closed form is chosen as a starting configuration. The overall minimum appears to be an open lp structure, which lies lower in energy than the dehydrated “closed” form by 29.9 kJ mol-1 per unit cell. This is in line with the experimental observation7 that the structure opens upon dehydration, suggesting that the open form is the thermodynamically preferred form when the pores are empty (albeit at elevated temperatures). Moreover, the small energy difference is consistent with the fact that the physisorption of only one or two water molecules, or CO2 molecules, per unit cell causes pore closure. However, the presence of two minima (presumably one global and one local) is at odds with the periodic DFT calculations described above, where only one minimum (large-pore) was found at constant pressure, with no evidence of any local minima corresponding to the np phase (though only two sets of lattice parameters were sampled). The force field could be amended so as to give only one minimum. In particular, the C1-C2-O1-Cr torsional term could be increased, reinforcing the planarity of the chromium/carboxylate moiety. However, although this has the effect of greatly increasing the energy difference between closed (e.g., calculated at constant volume) and open forms, more in line with the DFT energies, it also introduced other distortions into the structure. This approach was thus deemed unsatisfactory. In fact it is unlikely that the complex potential energy surface of the system is very accurately described by the potentials or even by DFT. The energy barrier between the two minima is almost certainly small, and the detection of the higher minimum could be dependent on a number of computational parameters, for example, those concerning the minimization algorithm. Upon careful reflection, therefore, we decided, rather than “overfitting” the force field, to persist with one that

Sorption-Induced Breathing in Metal-Organic Framework

J. Phys. Chem. C, Vol. 113, No. 2, 2009 549

Figure 2. Comparison of (a) experimental, (b) force-field-optimized, and (c) DFT-optimized crystal structures for CrMIL-53lp.

gave two minima, and a difference of around 30 kJ mol-1 per unit cell between lp and np. This had the added advantage that minimization calculations could be carried out in which sorbate molecules could be introduced into open and narrow-pore versions separately, structures could be minimized, and comparisons could be made (see Section 3.4). 3.3. Comparison of Optimized Structures. Cell parameters and structural details of the large- and narrow-pore minimized structures are given in Tables 2 and 3, respectively. For the large-pore system, we compare the experimental structure (MIL53ht7) with both force-field- and DFT-optimized periodic structures. The structures are also shown in Figure 2. Looking at the cell parameters, we see that the c parameter, along the direction of the inorganic chain, is very close for all three structures. In the DFT structure, the a and b parameters are slightly larger than in the experimental structure, on the order of 1%, with the cell volume correspondingly greater. We note that the experimental structure determination was carried out at 523 K, whereas simulations are at 0 K, and that framework structures are known to contract as well as expand with temperature. The force-field-minimized structure has a somewhat smaller and b larger, though the cell volume is very close to the DFT structure. a and b are the parameters that vary during pore breathingsas one contracts, the other expandssand it is clear that quite small energies would be required for minor cooperative distortions of the cell of this type. Furthermore, in the DFT structure there is a slight kink of the carboxylate groups, such that they are not coplanar with the aromatic ring. This is not reproduced by the force field. We also see that some inaccuracies in the experimental structure are corrected in the simulations: as mentioned earlier, the benzene rings of the benzenedicarboxylate ligands become planar and the carboxylate C-O bond length is more realistic (the experimental value of 1.36 Å seems too long). The Cr-O-C angle (139°) is much

closer to that from the DFT calculations (134°) than to the experimental value (152°). Again, it is clear that the experimentally refined coordinates of the organic linker, derived from powder diffraction, are not entirely accurate, and so this bond angle, like the carboxylate C-O bond length, should not be regarded as reliable. For comparison, we note that a very recent structure solution of the isostructural MIL-53(Al), obtained by neutron diffraction by Liu et al.,34 has Al-O-C angles of 137° in the large-pore form at 77 K, with the carboxylate C-O bond constrained at around 1.28 Å. Finally we note that, in general, bonds are somewhat longer in the DFT structure than the forcefield one, which can at least in part be attributed to the tendency of DFT to overestimate bond lengths. The narrow-pore force-field-minimized structure is compared to the experimental MIL-53lt structure7 in Table 3 and Figure 3. We do not compare with DFT since a narrow-pore minimum could not be obtained under constant pressure conditions. Again, the c parameter is closely reproduced. However, compared to the experimental structure, the simulated a parameter is larger and b smaller, leading to a considerably smaller cell volume. This is not very surprising considering that the experimental structure also contains four water molecules per unit cell (not shown in the figures). In the structures of MIL-53(Al) reported by Liu et al.,34 the dehydrated narrow-pore form had a ) 20.824 Å and b ) 6.871 Å at 77 K. At 295 K the values were a ) 20.756 Å and b ) 7.055 Å, emphasizing the temperature dependence of a/b even in the absence of large-scale pore breathing. Bond lengths and angle are otherwise closely reproduced, with the carboxylate C-O distances and C-O-Cr angles being again slightly lower than in experiment. Also reproduced is the alternating canting of the CrO6 octahedra going along c, and the concomitant rotation of the phenyl rings in the organic group, with respect to the carboxylate functions.

550 J. Phys. Chem. C, Vol. 113, No. 2, 2009

Coombes et al.

Figure 3. Comparison of (a) experimental and (b) force-field-optimized crystal structures for CrMIL-53np.

Overall, the force field reproduces the structures to a high degree of accuracy and captures all the main structural differences observed between the large- and narrow-pore versions. Furthermore, it gives the large-pore version as the energy minimum, with an energy of transformation from the narrowpore structure of -29.9 kJ mol-1. 3.4. Water Adsorption. Further minimization calculations were carried out in which increasing numbers of water molecules were introduced into the frameworks. For both the large- and narrow-pore structures, 1, 2, 3, 4, and 6 water molecules were placed at locations within the pores of the material. For low loadings, the most favorable sites were found to be adjacent to the hydroxyl group in the inorganic chain, with hydrogen bonding between the water and the OH group. The energies of adsorption of water were calculated according to

Eads )

Ecalc(hydrated MOF) - Ecalc(dehydrated MOF) n

where Ecalc(hydrated MOF) is the minimized energy of the particular hydrated material, Ecalc(dehydrated MOF) is the minimized energy of the large-pore dehydrated MIL-53, and n is the number of water molecules. The energy of isolated water molecules is zero in this force field. All structures were fully minimized at constant pressure, including all framework atoms. Under constant pressure conditions, all narrow-pore configurations minimized as narrow-pore, and all large-pore as largepore (with the exception of six water molecules in the largepore system for which no minimum could be found); that is, there was no pore opening or closing observed as a result of minimization. However, at each loading, one form is clearly lower in energy than the other. We do not necessarily expect a standard minimization technique to find the global minimum. The energies are given in Table 4: it is important to note that the values of Eads for both large- and narrow-pore forms are

TABLE 4: Calculated Energies of Adsorption of n Water Molecules in the Large- and Narrow-Pore Forms of CrMIL-53a dehydrated large pore narrow pore

0 29.9

n)1

n)2

n)3

n)4

n)6

-57.5 -58.2 -57.5 -56.9 -38.6 -64.7 -80.6 -84.4 -67.4 -68.5 -79.6 -90.6 -91.7 -72.4

a Energies in roman type are given in kJ mol–1 per water molecule relative to the dehydrated large-pore structure (except for dehydrated narrow pore, which is simply given per unit cell). For comparison, values for the narrow-pore form relative to the empty CrMIL-53np are given in italic type.

calculated with respect to the global minimum (that is, the dehydrated large-pore structure). This enables direct comparisons to be made between all the structures (np values relative to dehydrated np are given in italic type). Upon introduction of one water molecule, the energy of adsorption in the large-pore structure is calculated to be -57.5 kJ mol-1. By contrast, in the narrow-pore form the energy compared to the empty narrow-pore was -68.5 kJ mol-1 since a water molecule is able to interact with opposite sides of the pore. Crucially, however, with one water molecule adsorbed the closed form is still 18.9 kJ mol-1 unstable with respect to the open form; that is, for the adsorption of one water per unit cell, there is no thermodynamic force that would drive the open form to close. For two molecules, a different situation prevails. Here, for the large-pore version of the material, the overall adsorption energy remains closely the same at -58.2 kJ mol-1 per water molecule. In the narrow-pore, the adsorption energy is now -64.7 kJ mol-1. In other words is it now thermodynamically favorable for the structure to adopt the closed form. The energy difference between the two widens as three and then four molecules are introduced, and the energy of adsorption of water continues to increase up to -84.4 kJ mol-1 per water

Sorption-Induced Breathing in Metal-Organic Framework

J. Phys. Chem. C, Vol. 113, No. 2, 2009 551 with both sides of the pore and from ordering of the water molecules, and the pore closes. This is entirely consistent with the observation that the dehydrated large-pore CrMIL-53 is unstable in air and collapses to the narrow-pore form. However, upon thermal dehydration, the large-pore structure is recovered. 4. Conclusions

Figure 4. Force-field-optimized crystal structures for (a) CrMIL-53lp and (b) CrMIL-53np (fragment), with four water molecules per unit cell adsorbed in each.

molecule at a loading of four water molecules; that is, the same loading as found experimentally in the narrow-pore MIL53lt. The energy decreases in magnitude to -67.4 kJ mol-1 upon going to 6 molecules/unit cell., as steric constraints take effect in the small pores. The relatively constant value of the energy per water in the large-pore form, up to 4 molecules/unit cell, is due to the fact that there are four equivalent sites per unit cell, adjacent to the OH groups that project into the pore (see Figure 4a). These sites are so distributed that there is little selfinteraction between water molecules when there is one at each site, with O-O distances of 6.8 Å (along pore) and 8.5 Å (across pore). In the case of the narrow-pore version however, cooperative effects are seen to start at a loading of 2 molecules/unit cell as the energy per water increases, even with respect to the guest-free np form (italic type in Table 4). In Figure 4b we show a view along a channel of MIL53np at a loading of 4 molecules/unit cell in which the water molecules are clearly ordered, with each interacting both with the pendant OH groups (H · · · Ow distance of 1.5 Å) and with the carboxylate ligands (both water protons pointing toward a carboxylate O at 1.9 Å each). Both of these interactions could be classified as hydrogen bonding. Futhermore the water molecules, separated along the tunnel direction by 3.4 Å, alternate in their orientation, with their dipoles aligned successively up and down with respect to each other when looking along the channel direction (note that no symmetry was imposed during force-field minimizations). Overall, we see that the closure of the open pore upon the adsorption of a small quantity of water is explained entirely by the thermodynamics of physisorption: as water is increasingly adsorbed onto the hydroxyl sites, the narrow-pore from becomes more energetically favorable as energy is gained from interaction

In this work we have studied the structure and energetics of the chromium version of MIL-53 in both large- and narrowpore forms. The use of torsional and three-body terms involving chromium and the carboxylate function are needed to provide the correct energetic balance between the narrow- and largepore systems. By use of the force field, the GULP program is able to find local minima for hydrated and dehydrated versions of the large- and narrow-pore structures, which enable the thermodynamics of the breathing process to be further illuminated. The experimental observation that the open-pore form is favored on dehydration is rationalized by the fact that this is the global minimum. Upon introduction of water molecules, however, the narrow-pore form becomes favored at a certain concentration, leading to pore closure. Therefore it seems quite clear that, in the case of MIL-53, the phenomenon of pore breathing is driven overwhelmingly by the thermodynamics of physisorption. Moreover, this is only possible because there exist two forms of the material, close to each other in energy and easily interconvertible via a cooperative series of bond rotations. In trying to design another such material that might exhibit pore breathing, these are the factors that should be borne in mind when the inorganic units and interconnecting organic ligands are selected. Future extensions will investigate the phenomenon further by molecular dynamics simulations and by also considering carbon dioxide as a sorbate gas. Acknowledgment. This work was funded by the EU via the FP6 STREP project “DeSANNS” (SES6-CT-2005-020133). References and Notes (1) Kitagawa, S.; Kitaura, R.; Noro, S. Angew. Chem., Int. Ed. 2004, 43, 2334. (2) Rosseinsky, M. J. Microporous Mesoporous Mater. 2004, 73, 15. (3) Rowsell, J. L. C.; Yaghi, O. M. Microporous Mesoporous Mater. 2004, 73, 3. (4) Uemura, K.; Kitagawa, S.; Fukui, K.; Saito, K. J. Am. Chem. Soc. 2004, 126, 3817. (5) Hawxwell, S. M.; Espallargas, G. M.; Bradshaw, D.; Rosseinsky, M. J.; Prior, T. J.; Florence, A. J.; van de Streek, J.; Brammer, L. Chem. Commun. 2007, 1532. (6) Millange, F.; Serre, C.; Ferey, G. Chem. Commun. 2002, 822. (7) Serre, C.; Millange, F.; Thouvenot, C.; Nogues, M.; Marsolier, G.; Louer, D.; Ferey, G. J. Am. Chem. Soc. 2002, 124, 13519. (8) Loiseau, T.; Mellot-Draznieks, C.; Muguerra, H.; Ferey, G.; Haouas, M.; Taulelle, F. C. R. Chim. 2005, 8, 765. (9) Surble, S.; Serre, C.; Mellot-Draznieks, C.; Millange, F.; Ferey, G. Chem. Commun. 2006, 284. (10) Mellot-Draznieks, C.; Serre, C.; Surble, S.; Audebrand, N.; Ferey, G. J. Am. Chem. Soc. 2005, 127, 16273. (11) Loiseau, T.; Serre, C.; Huguenard, C.; Fink, G.; Taulelle, F.; Henry, M.; Bataille, T.; Ferey, G. Chem.sEur. J. 2004, 10, 1373. (12) Bourrelly, S.; Llewellyn, P. L.; Serre, C.; Millange, F.; Loiseau, T.; Ferey, G. J. Am. Chem. Soc. 2005, 127, 13519. (13) Llewellyn, P. L.; Bourrelly, S.; Serre, C.; Filinchuk, Y.; Ferey, G. Angew. Chem., Int. Ed. 2006, 45, 7751. (14) Ramsahye, N. A.; Maurin, G.; Bourrelly, S.; Llewellyn, P.; Loiseau, T.; Ferey, G. Phys. Chem. Chem. Phys. 2007, 9, 1059. (15) Serre, C.; Bourrelly, S.; Vimont, A.; Ramsahye, N. A.; Maurin, G.; Llewellyn, P. L.; Daturi, M.; Filinchuk, Y.; Leynaud, O.; Barnes, P.; Ferey, G. AdV. Mater. 2007, 19, 2246. (16) Greathouse, J. A.; Allendorf, M. D. J. Am. Chem. Soc. 2006, 128, 10678. (17) Greathouse, J. A.; Allendorf, M. D. J. Phys. Chem. C 2008, 112, 5795.

552 J. Phys. Chem. C, Vol. 113, No. 2, 2009 (18) Dauber-Osguthorpe, P.; Roberts, V. A.; Osguthorpe, D. J.; Wolff, J.; Genest, M.; Hagler, A. T. Proteins: Struct., Funct., Genet. 1988, 4, 31. (19) Amirjalayer, S.; Tafipolsky, M.; Schmid, R. Angew. Chem., Int. Ed. 2007, 46, 463. (20) Tafipolsky, M.; Amirjalayer, S.; Schmid, R. J. Comput. Chem. 2007, 28, 1169. (21) Salles, F.; Ghoufi, A.; Maurin, G.; Bell, R. G.; Mellot-Draznieks, C.; Fe´rey, G. Angew. Chem., Int. Ed. 2008, 47, 8487. (22) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (23) Perdew, J. P.; Ernzerhof, M.; Burke, K. J. Chem. Phys. 1996, 105, 9982. (24) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL06 User’s Manual; Università di Torino, Torino, Italy, 2006. (25) Catti, M.; Sandrone, G.; Valerio, G.; Dovesi, R. J. Phys. Chem. Solids 1996, 57, 1735.

Coombes et al. (26) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650. (27) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255. (28) Plant, D. F.; Maurin, G.; Bell, R. G. J. Phys. Chem. B 2007, 111, 2836. (29) Shubin, A. A.; Catlow, C. R. A.; Thomas, J. M.; Zamaraev, K. I. Proc. R. Soc. London, Ser. A 1994, 446, 411. (30) Cormack, A. N.; Lewis, G. V.; Parker, S. C.; Catlow, C. R. A. J. Phys. Chem. Solids 1988, 49, 53. (31) Saul, P.; Catlow, C. R. A.; Kendrick, J. Philos. Mag. B 1985, 51, 107. (32) Corà, F. Mol. Phys. 2005, 103, 2483. (33) Gale, J. D.; Rohl, A. L. Mol. Simul. 2003, 29, 291. (34) Liu, Y.; Her, J. H.; Dailly, A.; Ramirez-Cuesta, A. J.; Neumann, D. A.; Brown, C. M. J. Am. Chem. Soc. 2008, 130, 11813. (35) Sanders, M. J.; Leslie, M.; Catlow, C. R. A. J. Chem. Soc., Chem. Commun. 1984, 1271.

JP809408X