Hybrid Monte Carlo Simulations Combined with a Phase Mixture

Mar 23, 2010 - (1) These solids built up from inorganic clusters and organic linkers have ..... was represented by the rigid model developed by Harris...
1 downloads 0 Views 1MB Size
Download PDF
6496

J. Phys. Chem. C 2010, 114, 6496–6502

Hybrid Monte Carlo Simulations Combined with a Phase Mixture Model to Predict the Structural Transitions of a Porous Metal-Organic Framework Material upon Adsorption of Guest Molecules A. Ghoufi†,‡ and G. Maurin*,† Institut Charles Gerhardt Montpellier, UMR 5253 CNRS UM2, ENSCM, UniVersite´ Montpellier 2, Pl. E. Bataillon, 34095 Montpellier cedex 05, France, and Institut de Physique de Rennes, UMR 6251 CNRS, UniVersite´ Rennes 1, 263 aVenue du Ge´ne´ral Leclerc, 35042 Rennes, France ReceiVed: December 3, 2009; ReVised Manuscript ReceiVed: February 5, 2010

Anisotropic isobaric/isothermal molecular dynamics (MD) and grand canonical Monte Carlo (GCMC) techniques are combined in a hybrid scheme to get an osmotic Monte Carlo approach able to deal with a guest-assisted structural transition of a metal organic framework (MOF) porous solid corresponding to a large reversible breathing of its structure. This strategy based on (i) a consideration of a more general expression of the partition functions and (ii) a rigorous homogenization of the MD and MC parts allows us to capture the structural transition of the MIL-53(Cr) MOF-type solid in relation to the CO2 pressure. Further, we show that combining this revisited hybrid osmotic Monte Carlo (HOMC) approach to a newly developed “phase mixture” model, which is based on the existence of a pressure domain where several structural forms of MIL-53(Cr) are present, is an efficient way to accurately predict the adsorption behavior of this solid in the whole range of pressures. 1. Introduction The emergence of the metal organic frameworks (MOFs) systems constitutes one of the most exciting developments in recent nanoporous materials science.1 These solids built up from inorganic clusters and organic linkers have potential applications in many multidisciplinary areas including physics,2 chemistry,3 and biology.4 Among this class of solids, the MIL-53(Cr) (MIL: Materials of Institut Lavoisier) system is constituted of cornersharing CrO4(OH)2 octahedra chains, interconnected by transbenzenedicarboxylate ligands creating a 3D framework defining 1D diamond-shaped pores.5 This material presents the fascinating ability to “breathe” upon adsorption of various guest molecules.6,7 Such spectacular behavior, which raises both fundamental and application interests, has not been fully elucidated yet at the microscopic scale. It has been, however, evidenced that this breathing phenomenon is associated with a reversible structural switching between a narrow (Np) and a large pore (Lp) form, which implies large atomic displacements corresponding to a unit cell volume change from 1072 (Np) to 1486 Å3 (Lp), that is, 38% of the cell variation.8 This MOF solid initially present in the Lp form switches to the Np version upon CO2 adsorption at very low pressure before undergoing a second structural transition from the Np to the Lp forms starting at a pressure of ∼4 bar. These consecutive phase transitions between structures of different symmetry (orthorhombic and monoclinic for Lp and Np, respectively), which occur at a given CO2 pressure, are of “microscopic” type. Further, in situ X-ray powder diffraction measurements have revealed that these transitions are accompanied by a phase mixture (Lp, Np) domain whose composition varies with the pressure to obtain the purely Np or Lp form.9 The evolution of the phase mixture composition * To whom correspondence [email protected]. † Universite´ Montpellier 2. ‡ Universite´ Rennes 1.

should

be

addressed.

E-mail:

as a function of the pressure can be seen as a “macroscopic” transition, which is associated with a heterogeneous mixture of crystallites being either in the Np or Lp forms whether or not they are in contact with CO2. Molecular modeling thus appears as an invaluable tool to fully characterize both types of transitions. The challenges of our simulation approach are thus two-fold, (i) capturing the two-step microscopic structural transitions as a function of the CO2 pressure using a revisited hybrid Monte Carlo method based on a previously validated force field and a rigorous consideration of both expression and homogenization of the partition functions and (ii) predicting the CO2 adsorption isotherm in a wide range of pressures from an accurate description of the (Lp, Np) phase mixture composition along the macroscopic transition. 2. Theory and Methods Conventional modeling tools including molecular dynamics (MD) and grand canonical Monte Carlo (GCMC) techniques are not able to follow such a guest-assisted structural transition if one aims at rigorously mimicing the considered experimental conditions. On one hand, the GCMC simulations allow a fluctuation of the number of adsorbate molecules for a given chemical potential while usually considering the host framework rigid or with only tiny local structural rearrangements. Such an approach has been used, for instance, in a first approximation to calculate separately the CO2 adsorption isotherms for both the rigid Np and Lp forms of MIL-53(Cr) in order to build a simulated composite isotherm which is favorably compared to the experimental data.7 Snur et al.10 have followed another strategy consisting of switching between nearby rigid structures of silicalite-1 (∼3% of expansion/contraction of the unit cell volume) during the GCMC runs, which allowed them to successfully reproduce the shape of the experimental isotherms for the adsorption of benzene in this zeolite. However, such methodology cannot be applied in our case as the Np and Lp

10.1021/jp911484g  2010 American Chemical Society Published on Web 03/23/2010

Structural Transitions of a Porous MOF structures of MIL-53(Cr) strongly differ, as mentioned above, which renders complex the rescaling of the adsorbate molecules to match the swap of the host framework during the simulations. Indeed, we attempted to improve the efficiency of such MC simulations by coupling the anisotropic volume change in a triclinic symmetry with several intramolecular moves for the framework, including bond, valence, and torsion changes, rotation of (i) the phenyl rings around their centers of mass and (ii) between two phenyl rings around the axis connecting their centers of mass. Such effort failed to capture the structural switching of the MIL-53(Cr) system upon CO2 adsorption even using a very long simulation time. This was most probably due to the existence of cooperative intramolecular motions (rotation and torsion) within the framework which are not considered in a classical NpT move implemented in the MC simulations, the rescaling of the coordinates of the atoms implying a significant modification of the bond lengths and leading to a rejection of such a move. On the other hand, the MD approach is usually conducted for a fixed number of adsorbate molecules. As a typical illustration, our previous MD calculations based on our own force fields for the MIL-53(Cr) framework allowed us to predict which structural forms are present for a given number of CO2 molecules, while it was not possible to follow the structural switching when one continuously varied the concentration of the adsorbates.9 To overcome such a technical limitation, which enabled us to model the guest-assisted structural transition occurring under specific experimental conditions, the most appropriate approach would consist of combining MD and GCMC techniques in a hybrid osmotic Monte Carlo (HOMC) scheme to get an osmotic statistical ensemble (µNσT). To that purpose, two strategies can be envisaged. These include the following. (i) Implementation of a GCMD11,12 method in the osmotic ensemble (µN1σT),13 where N1 corresponds to the number of host frameworks (1 in our case for the MIL-53(Cr) solid). This approach has been mainly used in the past to probe the dynamics of various systems and thus extract diffusion, viscosity, and permeation data.14,15 (ii) The creation of a MD move which is part of the Markov chain in a HOMC route, allowing the volume fluctuation of the framework, which is crucial to capture such structural changes. This methodology allows an efficient exploration of the configuration space16 by means of the MD trial move which substitutes the volume change usually considered in the MC simulations. Further, one can also consider the anisotropic changes of both the shape and volume of the framework. Here, this latter approach has been selected to take up the great challenge to explore the breathing of the MIL-53(Cr) solid induced by the adsorption of guest molecules, corresponding to a magnitude of expansion/contraction of its unit cell volume that has never been addressed so far in any nanoporous system. One can also note that the complexity of our system also arises from a significant modification of the space group, that is, from an orthorhombic to a monoclinic symmetry, associated with this structural transition. This work significantly differs from the previous contributions reported by Chempath et al.17 and Banaszak et al.,18 where a homogenization of the MC and MD partition functions was not required as the investigated adsorbates, that is, n-butane in silicate-1 and ethylene in amorphous polyethylene, were considered flexible without any intramolecular constraint. Here, as the force field to describe the flexibility of the MIL-53(Cr) framework involves both bonded (constraints and intramolecular terms) and nonbonded terms,

J. Phys. Chem. C, Vol. 114, No. 14, 2010 6497 such simplification cannot be envisaged. Our first step thus consisted of proposing a more generalized expression of the MC partition function by including an extension of the model previously developed by Faller et al.16 It was then followed by a rigorous homogenization of the MC and MD partition functions in order to ensure a consistency between the configuration and momentum spaces and to get a reliable probability of acceptance for both the insertion/deletion and the change volume trial moves. The exploration of the configuration space (q) was conducted by considering the following: (i) A rotational contribution for the CO2 molecule via an inertia moment. The rotational partition function is considered only in some typical cases as, for instance, in liquid crystals, where this contribution plays a predominant role to accurately describe the different phase transitions and the surface anchoring. More generally, only a translational term is considered from the de Broglie wavelength.19 (ii) An intramolecular constraint bond for the MIL-53(Cr) solid by introducing a constraint Jacobian. Nobuhiro et al.20 have previously emphasized the importance of such an additional term to properly reproduce the structural behavior of polymer chains. This is even more relevant in our case as the MIL-53(Cr) exhibits a larger number of bond constraints than most of the polymers. The more general expressions considered for the MC and MD partition functions are reported in eqs 1 and 2. While for the MC partition function only the potential energy (U(rN)) is considered (eq 1), for the MD one (eq 2), the total Hamiltonian (Η(rN,pN)) needs to be included. ∞

2

MC QµN ) 1σT

i)1 N)0

2

MD QµN 1σT

)

[ ] qNi i

∑ ∑ exp(βµN) N ! ∫ |J | ∞

∑∑

i)1 N)0

i 1/2

exp(-β(U(rN, pN) +

i

exp(βµN)

[σ]V))drNdV

(1)

∫ exp(-β(H(rN, pN) + [σ]V))drNdpNdV (2)

H(rN, pN) ) U(rN) + K(pN) + K'(ωNi) K'(ωN1) ) 0 J2 ) 1

and

U(rN) ) UINTRA(rN) + UINTER(rN)

(3) (4)

In these equations, β corresponds to the reverse temperature, [σ] is the target determinant of the anisotropic pressure tensor, and ωN2 is the angular velocity of the adsorbate. U and H are the potential energy and Hamiltonian of the system respectively, V is the volume of the investigated framework,(rN,pN) are the vectors position and linear momentum of the system containing N atoms, and µ is the chemical potential. K and K′ are the translational and rotational kinetic energies, respectively, whose expressions have been previously reported,20 while in eq 4, UINTRA and UINTER are the intramolecular and intermolecular energies, respectively. Finally, N denotes the total number of particles present in the system corresponding to the sum N1 + N2, where i ) 1 and 2 are the labels for the framework and the adsorbate molecules, respectively, and qNi i are the corresponding total partition functions (i ) 1,2), which are defined by eq 5.

6498

J. Phys. Chem. C, Vol. 114, No. 14, 2010

Ghoufi and Maurin

qi ) qirotqitransqic

(5)

In this equation rot, trans, and c correspond to the rotational, translational, and Jacobian constraint contributions. As we introduced a bond constraint in the MC partition function, pNand rN are dependent variables, which means that a classical decomposition of the kinetic and the potential energy terms can not be realized. It thus requires the consideration of the determinant of the constraint matrix J.20,21 The expression of the MC partition function (see eq 1) is consistent with the one determined by Banazak et al.18 for the osmotic statistical ensemble. Further, to ensure a consistency between the MC and MD parts, we considered translational (CO2, MIL-53(Cr) framework) and rotational (CO2) contributions in the MC partition function. From the general expression of the MC partition function provided in eq 1 and using the classical expressions for each of these terms, the expressions of the eqs 6 and 7 can be obtained. One should note that for CO2, the constraint Jacobian was neglected as this adsorbate is treated as rigid in both GCMC and MD parts, implying that the internal conformation changes for each trial move (MD, translation, rotation, insertion, and deletion) vanish, thus leading to the absence of any energy variation of the constraint Jacobian for CO2. 2

MC QµN 1σT

)

∞

∑∑

ZitZir

exp(βµN)

i)1 N)0

∫

Zc1 )

VNi Ni!Λ3Ni |J1 | 1/2 |J01 | 1/2

Zr2 )

(

Zic

N

)

where ∆Gofn(rN1,pN1) corresponds to the energy difference associated with the Jacobian matrix (see Appendix A of ref 20). In our case, we used an algorithm which does not integrate the volume but considers the rotational part of the adsorbate molecules. Indeed, the probability of accepting the MD move is given by eq 10. It is worth noting that the same probability of acceptance is obtained if one considers the MD partition function expressed in eq 3. Indeed, the following equality Rnfo ) Rofn is respected.

{()

π(o f n) ) min 1,

4πI2kbT h

)

N/2

Zic ) 1

(7)

exp(-β(∆Uofn(rN) +

′ (ωN2) + [σ]∆Vofn + ∆Kofn(pN) + ∆Kofn

}

(7a) F(n)R(o f n) F(o)R(n f o)

)

(8)

In this latter equation, R(nfo) corresponds to the probability to generate a new trial state (n) starting with a previous old one (o). One can note that the translational and rotational contributions vanish when one expresses the probability of acceptance.

(10)

In this equation, ∆V ) Vn - Vo corresponds to the variation of the volume between the old (o) and the new configurations (n), which is created by the MD move. For the trial insertion/deletion move, the insertion and orientation bias21,22 was considered. This method consists of selecting in a first step an appropriate location for the center of mass of the inserted molecule by n test insertions of a Lennard-Jones (LJ) particle considering only the LJ interactions (ULJ); m orientations are then further tested for the selected location, and the resulting electrostatic and LJ interactions are computed (Uext). By using the MC partition function defined by eq 1, one obtains the probability of accepting the insertion as defined in eq 11, where ∆Ucorr represents the variation of the long-range corrections corresponding to the insertion of the molecule. Wi and Wo are the Rosenbluth factor20 of the statistical bias of insertion (i) and orientation (o). Further, the probability of accepting the deletion move can be thus expressed by eq 12.

π(N2 f N2 + 1) )

{

min 1,

4π2I2kbT |J21 | 1/2 V exp(βµ)WiWo × 2 1/2 Λ3(N2 + 1) h2 |J1,o |

MC QµN 1σT

(

N

}

(11)

}

(12)

exp(-β(Uext(rN) + ULJ(rN) + ∆Ucorr))

Zt1Zt2Zr1Zc1 exp(βµN) exp(-β(U(rN) + [σ]V))

π(o f n) ) min 1,

Vn Vo

N

In these equations, i corresponds to the (framework/adsorbate) pair, Λ is the de Broglie wavelength, I2 is the inertia moment of the CO2 molecule, and h is Planck’s constant. The superscript c is related to the bond constraint part. J1 and J10 are the real Jacobian for the framework (superscript 1) and a Jacobian for a conformation of reference, respectively. The so-obtained expression correspond to the particular case where the adsorbate molecule is linear and the components of the inertia moment are the same along the x, y, and z axes. The probability of accepting the MD trial move is thus given by considering the probability density (see eq 7a) and the detailed balance expressed by eq 8.

MC FµN 1σT

∆Gofn(rN1, pN1))) (9)

exp(-β(U(r , p ) +

i)2

with

′ (ωN2) + R(o f n) ) exp(-β(∆Kofn(pN) + ∆Kofn

∆Gofn(rN1, pN1)))

[σ]V))drNdV (6) Zit )

To consider the rotational and the Jacobian part when one switches from the MC to the MD part, we significantly modified the procedure given by Faller et al.,16 who introduced the kinetic part in the probability to generate a new configuration given by eq 9

{

π(N2 f N2 - 1) ) min 1,

exp(-βµ)

2 1/2 Λ3N2 h2 | |J1,o × V 4π2I k T |J2 | 1/2 II b 1

exp(β(Uext(rN) + ULJ(rN) + ∆Ucorr)) WiWo

3. Simulation Details In our simulations performed at 300 K, a Berendsen thermostat and an anisotropic barostat (with the following param-

Structural Transitions of a Porous MOF

Figure 1. Density of CO2 at 228 K and 0.1 bar as function of the MC steps calculated using the MC (black line) and HMC (red line) routes. The dashed line represents the experimental data taken in ref 28.

eters, τT ) 1.0 ps and τp ) 5.0 ps) were coupled with the QUATERNION and SHAKE-RATTLE algorithms and the velocity Verlet integrator.23–25 Recently, it has been shown that the weak coupling of the Berendsen integration algorithm (thermostat and barostat) is symplectic and reversible.26 This conclusion can be also drawn for the anisotropic Berendsen algorithm selected in this study, which allows a change of the size and the shape of the framework.23 The MIL53-(Cr) framework was described by our previously validated force field,9 while the CO2 molecule was represented by the rigid model developed by Harris and Yung.27 A simulation box consisting of 32 unit cells of MIL-53(Cr) was considered in both MC and MD simulations. The electrostatic interactions were calculated using the Ewald summation, while the van der Walls interactions were truncated at 12 Å.20 For the MC procedure, the ratios for each trial move were defined as follows: 0.1995 for the translation, 0.1995 for the rotation, 0.6 for the insertion/deletion, and 0.001 for the MD move (corresponding to 2/N, usually considered as the frequency of the volume change). The MD runs corresponded to 5000 steps (NMD) using a time step of 0.001 ps (tMD). We considered more than 70000 HOMC cycles corresponding to 100 ns of simulation time. To validate our findings, several ratios for each trial move and (NMD,tMD) couples have been tested in order to confirm that the simulations do not depend on the initial conditions. In addition, we have conducted a very long time simulation (3 × 105 cycles, that is, ∼450 ns on 16 processors) to confirm the stability of our algorithm. Validation of the Methodology. In a first step, we have tested the robustness of the MD move on a typical case which consisted of calculating the density of CO2. The values extracted from MC simulations performed in the NpT ensemble (involving the usual volume change) and HMC (hybrid Monte Carlo using a MD move in the classical (NpT) MC procedure) routes are F )1127.5 ( 26.1 and 1132.5 ( 6.2 kg · m-3, respectively, which both concur well with the experimental data F ) 1136.5 kg · m-3 (cf. Figure 1).28 It further emphasized that the MD move can efficiently substitute the volume change move usually considered in the MC simulations. Besides, as already mentioned by Faller et al.,16 one observed that the statistical deviation of F is smaller in the HMC scheme compared to that of the MC approach. Additionally, we have run in parallel GCMC and H-GCMC (using a MD move in a classical GCMC procedure) simulations, where CO2 was treated as rigid, to highlight that the creation/ annihilation trial moves are not perturbed by the introduction of the MD one. This later MD move consisted of sampling the phase space instead of only the configuration space. Figure 2 reports the resulting adsorbed amount of CO2 in MIL-53(Cr) as a function of the MC steps calculated at 300 K and 4.7 bar.

J. Phys. Chem. C, Vol. 114, No. 14, 2010 6499

Figure 2. CO2 adsorbed amount in MIL-53(Cr) obtained at 300 K and 4.7 bar as function of the number of MC steps. GCMC simulations (black dashed line) and H-GCMC scheme (red dotted line).

It is clearly observed that both approaches lead to the same mean values of 3.0 CO2/u · c This very good accordance allowed us to ensure that in the new scheme that we developed, (i) the ergodicity condition is respected and (ii) the probability of acceptance is accurately implemented. A further step consisted of computing the unit cell volume of MIL-53(Cr) and the CO2 adsorbed amount using our HOMC approach under the same conditions of temperature and pressure as those mentioned above. Figure 3a reports the evolution of the simulated data as a function of the HOMC simulation time. It can be observed that they converge toward 1055 Å3 and 3.0 CO2 molecules/ u · c, respectively, which reproduce very well the experimental data (1072 Å3 and 3.0 CO2 molecules/u · c). Moreover, our simulations confirm that the MIL-53(Cr) adopts its Np form in such conditions. To further test the validity of our HOMC approach, we calculated the correlation function for the density of CO2 confined within the pore of MIL-53(Cr). Figure 3b shows obviously that the resulting decorrelation is reached quickly (200 HOMC steps). We have also reported the temperature (calculated from the velocity) and the pressure (calculated from the virial expression) (cf. Figure 3c) as function of the MD simulation time. As can be observed, both mechanical and thermal equilibrium are ensured along the simulation even during the transition. We have also evidenced that the structural transformation can be related to a first-order transition which occurs after 200 ps, as shown in the energy profile reported in Figure 3d. Finally, in order to get a thermodynamic validation of our approach, we then simulated the adsorption enthalpy at 300 K for both the Lp and Np forms, which were then compared to those previously determined by microcalorimetry.8 To that purpose, we modified the usual formula reported by Parsonage and Nicholson for the grand canonical ensemble29 in order to take into account the fluctuation of the volume, leading to eq 13, which differs by the pV term in the total Hamiltonian. The resulting adsorption enthalpies of 18.2 and 33.7 kJ mol-1 for the Lp and Np forms, respectively, are in very good agreement with the experimental data (18 and 35 kJ mol-1) previously reported.8

∆h˙ ) kbT -

〈(U(r) + pV)N〉N1µσT - 〈(U(r) + pV)〉N1µσT〈N〉N1µσT 〈N2〉N1µσT - 〈N〉N1µσT〈N〉N1µσT

(13) Indeed, we have shown that the thermal, mechanical, and thermodynamic equilibria are respected, which fully validates our HOMC approach. Results and Discussion. The evolution of the unit cell volume of MIL-53(Cr) was first calculated up to 15 bar in order to span

6500

J. Phys. Chem. C, Vol. 114, No. 14, 2010

Ghoufi and Maurin

Figure 3. (a) Evolution of the simulated unit cell volume of MIL-53(Cr) (right axis) and of the number of adsorbed CO2 molecules per unit cell (left axis) as a function of the HOMC simulation time calculated at a pressure of 4.7 bar and 300 K. The dot and dot-dash lines represent the experimental amount of CO2 molecules and the unit cell volume, respectively, obtained under the same conditions of temperature and pressure. (b) Correlation function of the CO2 density calculated using the HOMC method. (c) Temperature (red line and left axis) and pressure (blue line and right axis) as function of the MD time. (d) Total potential energy (sum of the intramolecular and intermolecular terms) as a function of the MD time.

Figure 4. Evolution of the HOMC simulated unit cell volume of MIL53(Cr) along the whole adsorption process as a function of the pressure (open circle symbols) on a semilogarithmic scale. The solid circle symbols represent the experimental unit cell volume obtained for both Np and Lp forms.

the whole pressure range previously investigated by experiments (Figure 4). One observes that below 0.35 bar, the structure remains in its initial Lp form with a unit cell volume of 1496 Å3, very close to the one experimentally observed for the same range of pressures (1486 Å3).8 In contrast, a drop in the unit cell volume occurs for 0.35 bar, corresponding to a loading of 2.3 CO2/u · c, to attain a value of 1055 Å3. Our HOMC approach is thus able to capture this first structural switching Lp f Np for a pressure and CO2 loading very similar to that obtained from previous combined in situ XRPD and manometry experiments (p ) 0.30 bar and ∼2 CO2/u · c).8 When the pressure increases, one can observe only a small fluctuation of the unit cell volume up to 15 bar, which means that the MIL-53(Cr) structure remains in its Np version for the range of pressures considered. Our calculations thus fail to predict the reopening of the MIL-53(Cr) structure, that is, the second structural microscopic transition Np f Lp experimentally evidenced for a pressure of 6 bar, as illustrated in the experimental adsorption isotherm reported in Figure 5a. The Lp f Np transition failure might not be attributed to the force field as our previous MD

simulations have shown that it successfully reproduces the reopening of the MIL-53(Cr) structure for CO2 loading higher than 4 molecules per unit cell.9 Indeed, the absence of the Lp f Np transition in the present model should be due to the difficulty in exploring the full phase space, which can be explained by a highly orientational ordered arrangement of the confined CO2 molecules into the channel of the Np form. Indeed, we have previously shown that when absorbed in the Np form, the CO2 molecules are strictly aligned along the channel with a strong energetic double interaction with the opposing walls of the same pore.30,31 Additionally, the calculation of the order parameter32 confirms this highly confined geometry of CO2 in the Np form compared to the one in the Lp structure, highlighting for the first time a clear disorder/order transition of CO2, as reported in Figure 6. In that way, the insertion of additional molecules becomes more complex even by using the insertion statistical bias implemented in our HOMC approach or by increasing the number of trial orientations. Figure 5, which reports both experimental and HOMC simulated CO2 adsorption isotherms at 300 K, confirms such findings. While the experimental data are reasonably reproduced up to 4 bar, our HOMC simulations lead to a plateau in the adsorption isotherm, corresponding to a CO2 adsorbed amount close to the saturation capacity experimentally observed for the Np form. It thus appears that an accurate prediction of the adsorption isotherm does not seem possible without including additional information in our computational tool. We thus developed a complementary model based on previous experimental conclusions which clearly evidenced a phase mixture of the Np and Lp forms in the transition zone ranging from 4 to 8 bar.8 From the shape of the isotherms in the transition zone, it was thus possible to define the eq 14. sim sim N ) XNpNNp + YLpNLp

(14)

Structural Transitions of a Porous MOF

J. Phys. Chem. C, Vol. 114, No. 14, 2010 6501

Figure 5. Adsorption isotherms of CO2 in MIL-53(Cr) at 300 K: (a) experimental data (b), H-GCMC simulations (0), rigid Lp (4), and Np (3); (b) experimental data (b) and the combination of H-GCMC and “phase mixture” model (O).

Figure 6. Evolution of the order parameter (left axis and solid line) for the CO2 molecule and of the unit cell volume (right axis and dashed line) as function of the MD simulation time.

where for a given pressure range in the transition zone, N represents the total CO2 adsorbed amount, XNp and YLp are the sim sim and NLp correspond fractions of the Np and Lp forms, and NNp to the CO2 loading obtained from H-GCMC simulations performed on the rigid Np and Lp forms, respectively, for the same pressure. Considering that the sum XNp + YLp equals 1, eq 14 can be rewritten as follows (eq 15), where N only depends on the fraction of the Np form, XNp. sim sim N ) XNpNNp + (1 - XNp)NLp

(15)

In that way, we have two unknowns XNp and N for only one equation (eq 14). To solve this system, we first performed m GCMC simulations, considering each rigid form (Lp and Np) and various pressures (from p1 to pm) in the transition pressure zone of the isotherm with a step defined by (p1 - pm)/(m - 1). sim sim It was thus possible to compute NNp and NLp for all of these considered pressures. A set of m Np fractions XNp was further chosen in the range of [0-1], and we calculated for each of these XNp fractions m possible values of N using the expression of eq 15. In that way, one obtains (m × m) possible values of N. From the resulting square matrix (mm) (we considered as many GCMC simulations as XNp fractions), if one assumes a linear increasing evolution of N as a function of the pressure, as would be expected from the shape of the experimental isotherm when the transition occurs, the most simple solution is the one corresponding to the component of the matrix verifying i ) j, that is, its diagonal. One should bear in mind that this “phase mixture model” is built as previous experimental exploration has shown that the MIL-53(Cr) sample shows a

Figure 7. Fraction of the Np (solid line) and Lp (dashed line) forms of MIL-53(Cr) as a function of the pressure (semilogarithm plot).

heterogeneous structural behavior upon CO2 adsorption, which implies the existence of a phase mixture (Lp, Np) domain whose composition varies with the pressure. Indeed, this model does not have a particular physical meaning but aims to take into account experimental evidence that was neglected in most of the models developed recently33 which assumed only one type of crystallite undergoing a phase transition from both Lp to Np and Np to Lp at a given pressure. This phase mixture model was then employed to explore the region where the first microscopic structural transition Np f Lp occurs by considering m ) 10 in the range of [0.2-0.4] bar. Out of these two domains of pressure, only the purely Np or Lp are present, which means that the Lp/Np ratio is 0 or 1 depending on the pressure range. The so-extracted fractions of the Np and Lp forms as the function of pressure are shown in Figure 7. One can note that the two structural transitions Lp f Np and Np f Lp are clearly pointed out and that both very low and intermediate pressure ranges unambiguously show a phase mixture containing both Np and Lp forms. Finally, the combination of these data extracted from our phase mixture model and the H-GCMC adsorption isotherms calculated from the rigid Np and Lp forms leads to a resulting simulated adsorption isotherm which reproduces very well the experimental data, as shown in Figure 5b. This observation clearly emphasizes that investigating such a complex host/guest system cannot be envisaged without taking into account the existence of the phase mixture. 4. Conclusions The behavior of the MIL-53(Cr) solid upon CO2 adsorption has been decomposed into two “microscopic” transitions, corresponding to consecutive structural switchings between a

6502

J. Phys. Chem. C, Vol. 114, No. 14, 2010

Lp and a Np form. Both are accompanied by a phase mixture (Np, Lp) domain whose composition varies with the pressure. Our HOMC approach based on a more general expression of the partition functions and a rigorous homogenization of the MD and MC parts successfully reproduces the first microscopic transition Np f Lp occurring at low pressure, while it fails to capture the second one. Further, we have shown that the experimental adsorption isotherm can be only reproduced in the whole range of pressures by including a complementary macroscopic “phase mixture” model which consists of taking into account the existence of both Np and Lp forms in the transition zone. One should notice that the revisited HOMC method that we described here can be easily transferred to the investigation of the adsorption in other solids such as in zeolite (local relaxation), MCM-41 (roughness effect), or heterogeneous systems. References and Notes (1) Fe´rey, G. Chem. Soc. ReV. 2008, 37, 191. (2) Maspoch, D.; Ruiz-Molina, D.; Wurst, K.; Domingo, N.; Cavallini, M.; Biscarini, F.; Tejada, J.; Rovira, C.; Veciana, J. Nat. Mater. 2003, 2, 190. (3) Kitagawa, S.; Kitaura, R.; Noro, S. Angew. Chem., Int. Ed. 2004, 43, 2334. (4) Horcajada, P.; Serre, C.; Maurin, G.; Ramsahye, N.; Balas, F.; Vallet-Reghi, M.; Sebban, M.; Taulelle, F.; Fe´rey, G. J. Am. Chem. Soc. 2008, 130, 6774. (5) Serre, C.; Millange, F.; Thouvenot, C.; Nogue`s, M.; Marsolier, G.; Loue¨r, D.; Fe´rey, G. J. Am. Chem. Soc. 2002, 124, 13519. (6) Llewellyn, P. L.; Maurin, G.; Devic, T.; Loera-Serna, S.; Rosenbach, N.; Serre, C.; Bourrelly, S.; Horcajada, P.; Filinchuk, Y.; Fe´rey, G. J. Am. Chem. Soc. 2008, 130, 12808. (7) Ramsahye, N. A.; Maurin, G.; Bourrelly, S.; Llewellyn, P. L.; Loiseau, T.; Serre, C.; Fe´rey, G. Chem. Commun. 2007, 10, 3261. (8) Serre, C.; Bourrelly, S.; Vimont, A.; Ramsahye, N. A.; Maurin, G.; Llewellyn, P. L.; Daturi, M.; Filinchuk, Y.; Leynaud, O.; Barnes, P.; Fe´rey, G. AdV. Mater. 2007, 19, 2246. (9) Salles, F.; Ghoufi, A.; Maurin, G.; Bell, R. G.; Mellot-Draznieks, C.; Llewellyn, P. L.; Serre, C.; Fe´rey, G. Angew. Chem., Int. Ed. 2008, 47, 8487.

Ghoufi and Maurin (10) Snurr, R. Q.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1994, 98, 5111. (11) Duane, S.; Kennedy, A. D.; Pendleton, B. J.; Roweth, D. Phys. ReV. B 1987, 195, 216. (12) Mehlig, B.; Heermann, D. W.; Forrest, B. Phys. ReV. B 1992, 45, 679. (13) Escobedo, F. A. J. Chem. Phys. 1998, 108, 8760. (14) Duren, T.; Keil, F. J.; Seaton, N. A. Mol. Phys. 2002, 100, 3741. (15) Nagumo, R.; Takaba, H.; Nakao, S. I. J. Chem. Eng. Jpn. 2007, 40, 1045. (16) Faller, R.; De Pablo, J. J. J. Chem. Phys. 2002, 116, 55. (17) Chempath, S.; Clark, L.; Snurr, R. Q. J. Chem. Phys. 2003, 118, 7635. (18) Banaszak, B. J.; Faller, R.; De Pablo, J. J. J. Chem. Phys. 2004, 120, 11304. (19) Frenkel D.; Smit, B. Understanding molecular simulation; Academic Press:New York, 1996. (20) Nobuhiro, G.; Scheraga, H. A. Macromolecules 1970, 9, 535. (21) Ghoufi, A.; Goujon, F.; Lachet, V.; Malfreyt, P. J. Chem. Phys. 2008, 128, 154718. (22) Berendsen, H. J. C.; Postma, J. M.; Van Gunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (23) Miller, T. F.; Eleftheriou, M.; Pattnaik, P.; Nirango, A.; Newns, D.; Martyna, G. J. J. Chem. Phys. 2002, 116, 8649. (24) Melchionna, S.; Ciccoti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (25) Swope, W. C.; Andersen, H. C.; Berens, P. H.; Wilson, K. R. J. Chem. Phys. 1982, 76, 637. (26) Morishita, T. J. Chem. Phys. 2000, 113, 2976. (27) Harris, J. G.; Hung, K. J. Phys. Chem. 1995, 99, 12021. (28) Data taken from the saturation properties of carbon dioxide at the http://webbook.nist.gov. (29) Nicholson, D.; Parsonage, N. Computer Simulation and the Statistical Mechanics of Adsorption; Academic Press: London, 1982. (30) Ramsahye, N.; Maurin, G.; Bourrelly, S.; Llewellyn, P. L.; Serre, C.; Loiseau, T.; Devic, T.; Fe´rey, G. J. Phys. Chem. C 2008, 112, 514. (31) Salles, F.; Jobic, H.; Ghoufi, A.; Llewellyn, P. L.; Serre, C.; Bourrelly, S.; Fe´rey, G.; Maurin, G. Angew. Chem., Int. Ed. 2009, 42, 5044. (32) De Miguel, E.; Rull, L. F.; Chalam, M. K.; Gubbins, K. Mol. Phys. 1991, 74, 405. (33) Coudert, F. X.; Mellot-Draznieks, C.; Fuchs, A. H.; Boutin, A. J. Am. Chem. Soc. 2009, 131, 344.

JP911484G