First-Principles-Derived Force Field for Copper Paddle-Wheel-Based

Aug 12, 2010 - In the case of the dicopper formate model, we performed 20 independent GA runs. (population size 200; 100 generations) and a final “s...
0 downloads 0 Views 3MB Size
Download PDF
14402

J. Phys. Chem. C 2010, 114, 14402–14409

First-Principles-Derived Force Field for Copper Paddle-Wheel-Based Metal-Organic Frameworks Maxim Tafipolsky, Saeed Amirjalayer, and Rochus Schmid* Lehrstuhl fu¨r Anorganische Chemie 2, Organometallics and Materials Chemistry, Ruhr-UniVersita¨t Bochum, UniVersita¨tsstrasse 150, D-44780 Bochum, Germany ReceiVed: May 15, 2010; ReVised Manuscript ReceiVed: July 19, 2010

We present a fully flexible and ab initio-derived molecular mechanics force field for the ubiquitous copper paddle-wheel building block Cu2(O2C)4 in metal-organic frameworks. The force field expression is based on the established MM3 force field, extended by additional cross terms and specific bond-stretching and anglebending terms for the square-planar CuO4 coordination environment. Using reference data computed at the DFT level for nonperiodic reference systems, the parametrization is performed using an automated genetic algorithm optimization strategy in order to reproduce structure and low normal modes of the model systems. It is validated on the much investigated Cu-btc (HKUST-1) metal-organic framework. Beyond the structure, lattice-dynamic-dependent properties such as the bulk modulus and the observed negative thermal expansion effect of Cu-btc are quantitatively predicted by the force field without recourse with respect to experimental data. In connection with available parametrizations of various organic linkers, it can be useful for aiding the structure determination of known MOFs, but it also paves the way for the computational prescreening of yet unknown copper paddle-wheel-based frameworks. Introduction Metal-organic frameworks (MOFs), being a relatively new class of nanoporous materials, have already attracted considerable interest owing to their potential industrial applications in sensing, gas storage, chemical separations, and catalysis.1-3 In the meantime, a wide variety of such porous coordination polymers are known, with new structures appearing in the literature at a high rate. This is partly due to the wide variety of different inorganic connectors, the so-called secondary building units (SBUs).4 However, a large number of new network topologies and structural variations result from new or modified organic linkers combined with well-known SBUs. In this context, one type of SBU, namely, the Cu2(O2C)4 paddlewheel metal dimer, forming a four-coordinate square-planar vertex is extensively used. These systems can be synthesized relatively easily, and the remaining axial coordination site, usually referred to as an open metal site,5 can be advantageous for gas storage or in catalysis. One of the earliest and most investigated MOFs is the copper paddle-wheel-based Cu3(btc)2 (btc is 1,3,5-benzenetricarboxylate), abbreviated in the following text as Cu-btc (Figure 1).6 Many other copper paddle-wheelbased systems have been synthesized, as, for example, systems with a high hydrogen storage capacity.7 Another structurally interesting example is the chiral networks that could be realized by using a homochiral linker together with the copper paddlewheel SBU.8 The possibility to tune the structure and dimensionality of such frameworks by systematic modifications of the linker was exemplified on this inorganic node.9 Note that an additional class of porous coordination polymers can be formed by using ditopic nitrogen donor ligands as pillars connecting 2D sheets of paddle-wheel SBUs linked by dicarboxylic acids.10,11 Here, however, the open metal site is blocked. For a deeper understanding and a rational design of these porous materials, theoretical modeling has been of pivotal * Corresponding author. E-mail: [email protected].

Figure 1. Copper paddle-wheel-based Cu2(btc)3.

importance from the beginning.12-14 In particular, atomistic simulations allow a better understanding of the host-guest interactions, which are most relevant to technical applications. The methodology of Monte Carlo simulations in the grand canonical ensemble (GCMC) has been developed in the field of zeolite research,15,16 and the simulation protocols could be transferred to MOF systems.13 By analogy to the rather rigid zeolite frameworks, the network geometry of MOFs was assumed to be rigid and was taken from experimental data in the majority of investigations. However, one of the key points of MOFs is the conformational flexibility of both organic linkers and inorganic SBUs. In the case of guest molecule adsorption, this leads to very interesting dynamic behavior such as the socalled “breathing” effect, where the pore volume changes significantly,17 or the “gate opening” effect with strong hysteresis

10.1021/jp104441d  2010 American Chemical Society Published on Web 08/12/2010

Force Field for Cu Metal-Organic Frameworks in the adsorption isotherms.18 It is clear that rigid framework models are of no use in these cases. Unfortunately, quantum mechanical methods are numerically too demanding,19 and classical molecular mechanics force field (FF) methods lack the parameters needed for the inorganic SBUs.14 Sometimes, generic force fields with an atom-based parametrization, such as the UFF method,20 are employed. However, the accuracy of these methods for different systems is not clear or validated and, more importantly, cannot be systematically improved. Thus, more recently, fully flexible FFs for selected MOFs have been developed. Most of them are parametrized for the Zn4O(O2C)6 SBU in MOF-5 and the IRMOF family and were employed for the investigation of framework-dynamic-related properties such as (negative) thermal expansion, elastic constants, vibrational spectra, and thermal conductivity.21-26 Another system for which flexible FFs have been developed are the breathing MIL-53 systems.27 Here, the shrinking and expansion of the unit cell with increasing guest molecule loading could be reproduced. We have recently developed an automated strategy utilizing a genetic algorithm (GA) optimization method to derive accurate FF parameters in a consistent way from first-principlescalculated reference data.26 With this, we could parametrize an accurate FF for the Zn4O(O2C)6 SBU, which is able to reproduce the elastic constants of MOF-5 determined by periodic DFT methods.28 By employing the same parametrization strategy, we developed a first-principles-derived FF for the boron-containing 3D-periodic covalent organic framework COF-102.29 To our knowledge, such a fully flexible FF was used here for the first time to compute relative strain energies of both known and unknown network topologies, explaining the observed structure.30 This FF has meanwhile been used to investigate yet unknown COFs for their potential in hydrogen storage.31 For the above-mentioned copper paddle-wheel SBU, there is to our knowledge no accurately parametrized flexible force field available, presumably because of the difficulty in describing the coordination environment in this deformed square-planar Cu(II) system. Until now, theoretical investigations have concentrated on the Cu-btc system, with a particular focus on guest molecule adsorption and transport within the pores.32-41 To the best of our knowledge, all of these molecular simulations were performed by keeping the framework rigid. As a continuation of our previous work on the Zn4O(O2C)6 SBU, we present here a first-principles-derived flexible force field for the dicopper(II) tetracarboxylate core. Using the building block approach introduced previously for the case of MOF-5,26 this FF can be used in connection with parameter sets for di-, tri-, or tetracarboxylic linkers to simulate any copper paddle-wheel-based MOF. Note that the new force field currently excludes axial coordination to the open metal site. In a number of adsorption studies, the coordinative bonding to the copper atom was described by a combination of (partly readjusted) electrostatic charge-charge and van der Waals interactions.38 However, this coordination is clearly due to an orbital interaction with a strong directionality, which cannot be captured by two-body potentials. Preliminary tests of such a nonbonded FF14 for water ligation resulted in pronounced structural artifacts as a result of the lack of directionality of the Cu-Owater bond. Thus, the FF is intentionally restricted to the case of ligand-free “activated” MOFs. We are currently working on an extension of our FF to include ligand coordination in the future. The final FF is used to investigate properties such as elastic constants and the thermal expansion of the experimentally well characterized Cu-btc, which depend on the flexibility and dynamics of the framework. The results demonstrate the

J. Phys. Chem. C, Vol. 114, No. 34, 2010 14403

Figure 2. Model system used for the force field parametrization. Important bond lengths and angles are given in angstroms and degrees. Computed atomic charges for each atom type are shown in parentheses.

potential of such strictly bottom-up-derived FFs to predict even deformed network structures, lifting the constriction on experimentally observed frameworks. Computational Details The model systems used for the FF parametrization were fully optimized by density functional theory (DFT) using the hybrid functional B3LYP42-44 with the Gaussian 03 program package.45 The augmented correlation-consistent basis sets aug-cc-pVDZ46 were employed for all atoms, if not stated otherwise. Note that we found the augmented basis sets to be necessary to describe the conjugation of anionic carboxylates.47 For consistency, we used them throughout. For the Cu atoms, energy-consistent pseudopotentials denoted as (aug-)cc-pVDZ-PP,48-50 obtained from the EMSL Basis Set Exchange site,51,52 were employed. Because of the presence of low vibrational modes (on the order of 10 cm-1), we used a tight convergence criterion (opt ) tight) and a finer grid for the integration (int ) ultrafine). The optimized structures were confirmed to be true minima by vibrational frequency calculations. For all force field calculations, an adapted version of the Tinker suite of molecular mechanics programs was employed.53 Details of our GA implementation for force field generation can be found in ref 26. For the periodic Cu-btc, a full crystal energy minimization was performed without any symmetry constraint (space group P1). A cutoff value of 12 Å for the van der Waals interactions (switched to zero at the cutoff radius) was used, and the electrostatic interactions were computed by a smooth particle mesh Ewald summation with a real space cutoff value of 12 Å. Molecular dynamics simulations in the NPT ensemble were performed with an integration time step of 1.0 fs, a 0.1 ps thermostat relaxation time, and a 2.0 ps barostat relaxation time using the Berendsen bath coupling method54 as implemented in Tinker. The infrared spectrum of the periodic Cu-btc was calculated by a Fourier transform of the dipole moment time derivative autocorrelation function, derived from an MD simulation at room temperature in the NVE ensemble, following previous work.25 Results and Discussion Reference Systems. On the basis of our previous experience, where the necessary FF parameters for MOF-5 could be successfully derived using the basic zinc formate as a model system,26 we used the most simple model system for a dicopper paddle-wheel unit, namely, the dicopper tetraformate shown in Figure 2. Dinuclear paddle-wheel carboxylates are ubiquitous in metal-organic chemistry.55 A number of second-row transition-metal-containing compounds have been investigated both experimentally and theoretically, with a special emphasis on

14404

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

Tafipolsky et al.

TABLE 1: Selected Structural Data for Dinuclear Copper Tetracarboxylates formate: Cu2(O2CH)4 method

B3LYP

r(Cu-Cu), Å r(Cu-Ocarb), Å r(Cu-Ow), Å r(C-O), Å ∠O-C-O, deg

2.518 1.967 1.265 127.8

a

a

MP2

b

formate: Cu2(O2CH)4(H2O)2 a

FF

B3LYP

2.504 1.935

2.519 1.967

1.274 127.7

1.265 127.8

2.604 1.980, 2.001 2.270 1.262, 1.265 127.7

acetate: Cu2(O2CCH3)4(H2O)2 B3LYPa

exp. (XRD)c

2.564 1.974, 1.998 2.302 1.268, 1.272 125.4

2.613 1.947,d 1.990 2.161 1.260d 125d

Triplet state; UB3LYP/aug-cc-pVDZ-PP. b UMP2/aug-cc-pVDZ-PP. c Experimental values from ref 69. d Averaged values.

the metal-metal interactions. For example, a DFT investigation of molecular structures and vibrational spectra of dinuclear tetracarboxylates (formate and acetate) of molybdenum and rhodium confirms the reliability of the DFT method.56 Much less attention has been devoted to the theoretical investigation of the first-row transition-metal carboxylates (excluding chromium) such as copper and zinc compounds, with the general belief that no metal-metal bond is formed in their dinuclear carboxylates. Dimeric copper(II) carboxylates have been known for a long time, since the discovery of unusual magnetic properties of copper acetate monohydrate, Cu2(O2CCH3)4(H2O)2, whose antiferromagnetic character was explained by its dimeric structure,57 which was predicted58 before it could be confirmed.59 Magnetic coupling between metal centers in this particular compound was quite early a subject of theoretical study.60 Its triplet excited state is very close to the open-shell singlet ground state, and the energy splitting is a very sensitive probe for judging the performance of different levels of theory.61-64 The geometry of Cu2(O2CH)4 was optimized using the spinunrestricted hybrid DFT method (UB3LYP) and applying the so-called broken-spatial-symmetry (BS) approach.65,66 Because of the antiferromagnetic exchange coupling, the high-spin (triplet) excited state of the dimeric copper(II) acetate is experimentally found to be only 292.2 cm-1 above the openshell singlet ground state (Ci symmetry),67 whereas our calculations (UB3LYP/aug-cc-pVDZ-PP) gave a value of 372 cm-1 for the Cu2(O2CCH3)4(H2O)2 acetate, taken as 2 times the (geometry-optimized) energy difference between the high-spin and the BS states, 2(EHS - EBS). (See also ref 61.) Despite the accuracy of this energy separation, the open-shell singlet ground state shows sizable spin contamination (S2 ) 0.98 for the formate model), which might affect its use as a reference.68 Because both the geometries and the vibrational frequencies of these two states are nearly indistinguishable for the purpose of our force field parametrization, we used the triplet state as the reference because it is also more straightforward to calculate. In addition, the triplet state is appreciably populated at room temperature. Note, however, that the closed-shell singlet state was calculated to lie energetically far above the triplet state (12860 cm-1) and should not be used as a reference. In Figure 2, the optimized structure with relevant geometrical data is shown. To assess the quality of our DFT reference structure, in Table 1 a comparison with experimental data was performed. To facilitate this, the computed results for the water-coordinated formate and acetate are also given. An important point when comparing experimental data and theoretical results determined in vacuum is the hydrogen bond interactions between the axial water ligand and the carboxylate oxygens. For the theoretically optimized structures with axial water, the Cu-Cu-Owater angle is significantly bent in order to allow for an intramolecular interaction with the carboxylate oxygens, which is not observed

in the neutron diffraction study of the dicopper acetate.69 It is purely an artifact of the isolated molecule calculation, where the presence of other surrounding molecules is ignored and leads to a reduction in symmetry and two different Cu-Ocarb bond lengths. However, the different intermolecular hydrogen bonding interactions in the experimental structures also result in a scattering of the observed Cu-Ocarb bond distances, making it hard to compare properly. We note here that in the assynthesized Cu-btc material6 no bending of the Cu-Cu-Owater is seen and all Cu-Ocarb distances are of the same length (within experimental error). Obviously, an extensive hydrogen-bonded network with some additional water molecules is formed. As observed previously,26 the metal-oxygen bond lengths are somewhat overestimated on the DFT level as compared to the corresponding MP2 results (UMP2/aug-cc-pVDZ-PP). The most noticeable change on axial coordination with water is the fact that both the Cu-Cu distance and the Cu-Ocarb distances are elongated by ca. 0.08 Å and 0.01-0.03 Å, respectively. This effect has recently been observed for the Cu-btc network.70 Our calculations on the paddle-wheel unit without axial water give the Cu-Cu distance in good agreement with that found in the fully desolvated Cu-btc framework (2.50 Å70,71). Note that the Cu-Ocarb bonds in solid copper acetate monohydrate, which are not perturbed by hydrogen bonds, are similar in length (1.942(1) and 1.952(1) Å69) to those found in the Cu-btc framework (1.952(3) Å6). From this comparison, we can conclude that effects such as the systematic error of the given level of theory, the simplified formate model system, or additional hydrogen bonding on the structural parameters are small. Thus, the formate model computed on the DFT level can serve well as a reference for our force field parametrization, as we seek for the minimum structure and curvature of an “unbiased” potential energy surface. As in our previous work,26 we intentionally do not correct for deficiencies of the DFT method but fit to reproduce this higher-level method in the sense of a multiscale simulation approach. Accordingly, the calculated Hessian matrix elements were taken unscaled for the parameter fit. Following the protocol developed previously,26 the effective atomic charges were obtained from the fitting of the electrostatic potential (ESP) based on the Merz-Kollman sampling scheme72 using ca. 2000 grid points per atom. (For a critical comparison of methods for charge derivation, see ref 14.) The van der Waals exclusion radius of 2.26 Å for the Cu atom was taken from ref 73. The resulting atomic charges are given in parentheses in Figure 2. Our calculated charges are in good agreement with those obtained in previous studies with slightly different basis sets and fitting protocols.35,41 To describe the Cu-btc system, additional parameters, which are missing in the standard MM3 parameter set, were derived from a sodium benzoate model, which was optimized on the same B3LYP/aug-cc-pVDZ level of theory, following our previous study on carboxylate force fields.47

Force Field for Cu Metal-Organic Frameworks

J. Phys. Chem. C, Vol. 114, No. 34, 2010 14405

TABLE 2: Force Field Parameters for the Copper Paddle-Wheel Core and Carboxylate Linkera bond stretch

reference distance, Å

force constant, mdyn/Å

Cu-Cu Cu-Ocarb Ccarb-Ocarb Ccarb-Cph Cph-Cph

2.656 1.976 1.262 1.503 1.38b

0.25 1.27 9.71 4.445 6.56b

Morse potential R, Å-1 2.5

in-plane angle bending

reference angle, deg

force constant, mdyn Å/rad2

Ocarb-Cu-Ocarb Cu-Ocarb-Ccarb Ocarb-Ccarb-Ocarb Ocarb-Ccarb-Cph Cph-Cph-Ccarb

(90, 180)c 118.71 126.84 117.52 115.31

0.44c 0.26 1.61 0.999 0.723

out-of-plane bending (central-apex)

reference angle, deg

force constant, mdyn Å/rad2

Ccarb-Ocarb Ccarb-Cph

0.0 0.0

1.313 1.313

cross term (stretch-stretch and stretchbend)d Ocarb-Ccarb-Ocarb Ocarb-Cu-Ocarb Cu-Ocarb-Ccarb Ocarb-Ccarb-Cph Cph-Cph-Cph Ccarb-Cph-Cph

stretch1-bend, stretch2-bend, stretch1-stretch2, mdyn/rad mdyn/rad mdyn/Å 0.479

0.479

0.147 0.5 0.124 0.353

0.026 0.4 0.124 0.092

1.461 0.22 0.22 0.6 1.061 0.293

torsion

n

barrier height, kcal/mol

Cu-Ocarb-Ccarb-Ocarb Cph-Cph-Ccarb-Ocarb Cph-Cph-Cph-Ccarb Ccarb-Cph-Cph-H Cph-Cph-Cph-Cph Cph-Cph-Cph-H H-Cph-Cph-H

2 2 2 2 2 2 2

6.62 1.75 5.983b 5.983b 5.983b 5.385b 6.881b

van der Waals terme

radius, Å

well depth, kcal/mol

bond length reduction factor

Cu Ocarb Cph Ccarb H

2.26 1.82 1.96 1.94 1.62

0.296 0.059 0.056 0.056 0.020

0.923

a

Only parameters deviating from the original MM3 parameter set are given. b These values were adjusted by performing an SCF-MO calculation for the π system. c Specific bending potential for a square-planar conformation; see the text. d The reference distances and angles are identical to the regular bond stretch and angle bending parameters. e The parameters are taken from ref 73.

Parametrization. The parametrization strategy focuses on accuracy in reproducing the reference data at the expense of transferability. The primary step is the choice of the nonbonded parameters. van der Waals (vdW) interactions are described by the MM3 specific Buckingham potential74 and the usual Lorentz-Berthelot mixing rules, with all parameters taken unadjusted from ref 73. (For a listing, see Table 2.) These generic vdW parameters are available for the full periodic table

and are the only values derived or extrapolated from experimental data. Note that the primary intention of our force field is to reproduce the conformational flexibility of the network itself, and it is not specifically adjusted for the weaker intermolecular host-guest interactions. However, with this approach we could successfully reproduce the diffusivity of benzene MOF-5 in previous work.75 In the future, the parameter set might have to be fine-tuned selectively for intermolecular interactions, but doing this in a purely first-principles-based bottom-up strategy is computationally demanding and has until now been performed only for small gas molecules such as hydrogen or methane.14 All other parameters have been derived from the corresponding first-principles calculations. In contrast to the MM3 energy expression, where bond dipoles are used to describe electrostatic interactions, atomic point charges are used in our force field (with interactions between 1,4-connected atoms scaled by a factor of 0.5). To arrive at a charge-neutral Cu2(O2C)4 fragment, the small charge on the formate hydrogen was added to Cu and Ocarb, leading to the following charges used in the force field calculations: Cu, +1.2; Ocarb, -0.64; and Ccarb, +0.68. For the remaining atoms of the organic linkers, we use the following atomic charges: Cph, -0.12 and H, +0.12 for all Cph-H groups and a zero charge for all other Cph species. A general issue in modeling coordination compounds is the large partial charges needed to represent the electrostatic potential in the surrounding region. In particular, in the copper dimer the charge of +1.2 on the two metal atoms in close proximity of only about 2.5 Å leads to a strong electrostatic repulsion, which is an artifact of the atomic point charge model.76 This leads to carboxylate bending potentials that are too stiff to compensate for this repulsion and maintain the structure. To avoid this problem, we introduced an explicit Cu-Cu bond to exclude this 1,2-interaction. Note that this is only for modeling the system properly and has nothing to do with the true presence of a chemical bond between the metal atoms. The functional form of the energy expression for the bonded interactions is a modified form of the MM3(2000) force field of Allinger.74 We have modified the original MM3 stretch-bend cross term by introducing two different force constants for the two bonds that comprise a bond angle and implemented a new stretch-stretch cross term.26 The weaker Cu-Ocarb coordination bonds are described by a Morse potential Morse Estretch )

kstretch 2R2

[1 - e-R(r-r0)]2

(where kstretch is the force constant, r0 is the reference distance, and R ) 2.5 Å-1), which was found to be superior to the standard MM3 quartic term used for all other bonds.14 For the square-coordination environment of Cu(II) with angular minima at 90 and 180°, an improved Fourier-type angle-bending potential is used77 (discussed in more detail in ref 78)

1 Ebend(θ) ) kbend[1 + cos(θ)][1 + cos(2θ)] 2 where θ is the Ocarb-Cu-Ocarb angle and only kbend needs to be adjusted. To fit the necessary additional parameters, the DFT-optimized geometries of the reference systems and the Cartesian Hessian matrices were projected into a redundant set of internal coordinates defined by the force field connectivity. A genetic algorithm was used to optimize the FF parameters to reproduce

14406

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

Tafipolsky et al.

TABLE 3: Vibrational Modes (in cm-1) below 300 cm-1 for the Dicopper Tetraformate Model System UB3LYP/ aug-cc-pVDZ-PP

force field

deviation

mode degeneracy

58 67 106 154 164 222 225 232 251 274

91 124 118 161 164 219 179 229 229 290

-33 -58 -13 -8 0 +4 +45 +3 +22 -16

1 1 2 2 1 1 1 2 1 1

the reference data. For the objective function, both geometry and Hessian data (diagonal and selected off-diagonal elements) for the parametrized internal coordinates were included in a weighted root-mean-square deviation scheme. In the case of the dicopper formate model, we performed 20 independent GA runs (population size 200; 100 generations) and a final “shoot-out” run using the best individuals from the previous runs. This protocol reliably results in a very high final fitness. For further details on the fitting strategy, see ref 26. The linker model system was used to parametrize the Ccarb-Cph interaction including stretching, bending, out-of-plane bending, torsion, and cross terms. For torsional barriers, in addition to the curvature at the energy minimum, rotational energy barriers were also included in the fitting procedure. Following the MM3 scheme for aromatic systems, the parameters for the Cph-Cph bond stretch force constant and the 2-fold torsion potential were adjusted by performing a self-consistent-field molecular orbital calculation for the π system using a modified Pariser-ParrPople method.79 However, this was not done “on the fly”; the parameters were kept fixed. In Table 2, all parameters are given. Validation. At first, we directly compare the results of the dicopper formate model (Figure 2) with the DFT reference. As expected, the structure is reproduced very well as shown in Table 1. The low vibrational frequencies below 300 cm-1 are compared in Table 3. The largest relative deviations are observed for the two lowest modes, which can be described as a twist around the Cu · · · Cu axis and a D4h to D2d symmetry distortion from a square-planar to a tetrahedral copper coordination environment, respectively. These discrepancies in the normal modes persist, despite the low deviations between elements of the redundant internal coordinate Hessian, included in the fitting, from the force field and the DFT reference, respectively. However, a close inspection of the off-diagonal terms of the reference Hessian shows, for example, strong coupling between the Cu-Cu and the Ocarb-Ccarb stretch with a force constant above 1 mdyn/Å. This reveals the peculiar electronic interactions in a coordination complex, which cannot fully be captured by a force field without a stretch-stretch cross term for nonadjacent bonds. In addition, we presume that the large electrostatic interactions between highly charged atomic sites in close proximity might also lead to a potential that is too stiff.76 However, the overall agreement is gratifying in view of the fact that only a limited number of cross terms are included. In particular, the frequency of the normal mode best described as the Cu · · · Cu stretch (222/219 cm-1) and the soft in-plane deformations of the square SBU (154/161 cm-1 and 164/164 cm-1, bending of Ocarb-Cu-Ocarb), which are most relevant to the flexibility of the MOFs, are very well reproduced by the force field. In the next step, the FF was used to investigate the experimentally well characterized periodic Cu-btc, which rep-

resents a 3,4-connected network in the tbo topology.80 The optimized lattice size of 26.3833 Å (corresponding to the zero Kelvin limit) is in good agreement with the value of 26.3046(2)81 for a fully desolvated Cu-btc measured at low temperature (T ) 77 K). Our results suggest that the lattice constant is particularly sensitive to the Cu-O length, and the discrepancy can be explained by the overestimation of this distance by the DFT method (1.97 Å) as compared to the experimental value (1.92 Å).81 A similar observation has been made before in our study on MOF-5.26 We intentionally accept this slight deviation, which is due to the first-principles reference, and refrain from a manual “tuning” of the parameters after the automated GA fitting procedure. In the tbo network, the trigonal vertex has a 3-fold symmetry axis, whereas the four-coordinate vertex has reduced local D2d symmetry with two different angles of 109.5 and 70.5°. As a consequence, a deformation is needed to accommodate the square-planar copper paddle-wheel SBU, leading to the typical bowing of the linker as shown on the right in Figure 3. In contrast to the 19.5° deviation in the parent network, Ccarb-XCu2-Ccarb (where XCu2 is the centroid of the two Cu atoms) deviates only by 2.4° from orthogonality. The remaining deformation is transferred to an out-of-plane bending of the carboxylate “hinge” and a bending of the carboxylate out of the phenyl ring plane. However, relaxing a cut-out paddle-wheel unit with the adjacent phenyl rings (shown on the left in Figure 3) reveals that a strain energy of only 1.1 kcal/mol is stored in this fragment, which is well accessible by thermal fluctuations. A comparison of the Ccarb-XCu2-Ccarb and the XCu2-Ccarb-Cph angles with the corresponding experimental data (in parentheses in Figure 3) in ref 81 demonstrates the accuracy of the FF. For completeness, we also computed the IR spectrum of Cubtc using the dipole moment time derivative autocorrelation function in order to estimate the IR intensities of the power spectrum. In Figure 4, the result is compared with an experimental measurement for an activated Cu-btc sample, showing similar peaks as reported in the literature.82 Overall, even the smaller features of the spectrum above 400 cm-1 are well reproduced. In particular, the strongest peaks related to the Ocarb-Ccarb-Ocarb bending (760 cm-1) and the Ccarb-Ocarb stretching (asym, 1645 cm-1; sym, 1440 cm-1) are well reproduced. Note that a comparison of these higher normal modes, computed on the DFT and FF levels for the formate reference system, reveal the necessity of further cross terms also for nonadjacent bonds, in order to capture the coupling of the symmetric and asymmetric Ccarb-Ocarb stretching modes accurately.47 We speculate that this, and the Cph-Cph parameter taken unfitted from MM3, are the cause for the smaller and slightly shifted peaks in the region above 1200 cm-1 in the computed spectrum. However, these medium- and highfrequency normal modes are rather stiff and contribute only marginally to lattice deformation and dynamics compared to the low modes below 400 cm-1. To test the validity of our flexible force field for dynamic properties, we computed the thermal expansion coefficient of desolvated Cu-btc. In Figure 5, the temperature dependence of the lattice constant is shown and was obtained from a series of NPT ensemble simulations between 200 and 500 K at a pressure of 1 atm. At each temperature, the system was subjected to 100 ps of temperature equilibration at constant volume (NVT), followed by a 100 ps equilibration in the NPT ensemble and a 400 ps sampling run (NPT). Despite the above-discussed offset in the lattice parameter, our simulations give a volumetric thermal expansion coefficient of -1.07 × 10-5 K-1 and

Force Field for Cu Metal-Organic Frameworks

J. Phys. Chem. C, Vol. 114, No. 34, 2010 14407

Figure 3. Copper paddle-wheel fragment with adjacent phenyl rings cut out of Cu-btc shown on the right, demonstrating the deformation necessary to adopt to the tbo network topology. Ccarb-XCu2-Ccarb and XCu2-Ccarb-Cph angles are given in degrees; experimental data in parentheses are from ref 81.

In addition, we computed the elastic constants of the material, which can be important for its technical application. The bulk modulus, B, is calculated to be 25 GPa (computational details in ref 26), in excellent agreement with the experimental value of about 30 GPa for the nondesolvated Cu-btc.85 In accordance with the lower NTE, the predicted bulk modulus of Cu-btc is appreciably larger than that found in MOF-5 (14 GPa),26 indicating its higher stiffness. Note that in contrast to the discrepancy between experimental results determined by nanoindentation and theoretical results in the case of MOF-5,28 here good agreement between the measured and computed bulk moduli is found. The shear modulus (C11 - C12) of Cu-btc is computed to be 3.5 GPa, and the strain modulus C44 is 13 GPa. Conclusions Figure 4. Comparison of the computed IR spectrum (bottom) with the measured (top) spectrum for a desolvated, activated Cu-btc sample.83

Figure 5. Temperature-dependent variation of the simulated lattice constant of Cu-btc with standard deviations (0). The value at T ) 0 K results from the full optimization of Cu-btc. Experimental data are from synchrotron powder X-ray diffraction measurements (2).71

quantitatively reproduce84 the negative thermal expansion (NTE) behavior of the desolvated Cu-btc material as measured recently in the experimental studies71,81 (Figure 5), which is, however, significantly smaller than that found for MOF-5.22 This result is a clear indication of the accuracy of the first-principles-derived FF, which is able to reproduce such a complex global framework motion without the use of experimental data in the parametrization.

The MM3-based force field presented here is, to our knowledge, the first ab-initio-derived fully flexible force field for copper paddle-wheel-based MOFs. The energy expression is extended with specific cross terms and bond-stretching and angle-bending terms, which are suited to the modeling of the pseudosquare copper coordination environment. Following our previous work on Zn4O-based MOFs,26 we used a genetic algorithm parameter fitting method. Structure and Hessian matrix reference data of nonperiodic models were calculated on the DFT level (B3LYP/aug-cc-pVDZ level), which was verified to be sufficiently accurate in comparison to experimental data and higher-level theory. For validation, the resulting force field was employed to compute the structure, IR spectrum, elastic constants, and thermal expansion of the well-known Cu-btc (HKUST-1).6 An indication of the accuracy of the force field is the fact that the measured NTE71,81 is quantitatively reproduced. The origin of this effect is the delicate global motion of the framework and shows that in particular the low normal modes responsible for the network deformation are well parametrized. In addition, the measured bulk modulus of Cubtc85 is well reproduced by the force field. In conclusion, we could show that an accurate force field can be used to compute network structures and complex lattice dynamics of copper paddle-wheel-based metal-organic frameworks. In combination with parametrizations for other organic linkers, it can be employed for arbitrary systems beyond the ones investigated here, such as those where recently a symmetrypreserving isomerism has been observed.86 We are currently

14408

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

using the new force field to explore possible topologies of 3,4connected networks such as Cu-btc. With its predictive potential, it is useful both in the structural characterization of alreadyknown MOFs and in computing yet unknown structures (e.g., those for a theoretical prescreening of host-guest interactions). Acknowledgment. The Deutsche Forschungsgemeinschaft (SFB 585 and SPP 1362) and the Alfried Krupp von Bohlen und Halbach Stiftung are gratefully acknowledged for their financial support of this project. References and Notes (1) Rowsell, J. L. C.; Yaghi, O. M. Microporous Mesoporous Mater. 2004, 73, 3. (2) Ferey, G. Chem. Soc. ReV. 2008, 37, 191. (3) Kitagawa, S.; Kitaura, R.; Noro, S. Angew. Chem., Int. Ed. 2004, 43, 2334. (4) Tranchemontagne, D. J.; Mendoza-Cortes, J. L.; O’Keeffe, M.; Yaghi, O. M. Chem. Soc. ReV. 2009, 38, 1257. (5) Murray, L. J.; Dinca, M.; Long, J. R. Chem. Soc. ReV. 2009, 38, 1294. (6) Chui, S. S. Y.; Lo, S. M. F.; Charmant, J. P. H.; Orpen, A. G.; Williams, I. D. Science 1999, 283, 1148. (7) Lin, X.; Telepeni, I.; Blake, A. J.; Dailly, A.; Brown, C. M.; Simmons, J. M.; Zoppi, M.; Walker, G. S.; Thomas, K. M.; Mays, T. J.; Hubberstey, P.; Champness, N. R.; Schroeder, M. J. Am. Chem. Soc. 2009, 131, 2159. (8) Ma, L.; Lin, W. J. Am. Chem. Soc. 2008, 130, 13834. (9) Furukawa, H.; Kim, J.; Ockwig, N. W.; O’Keeffe, M.; Yaghi, O. M. J. Am. Chem. Soc. 2008, 130, 11650. (10) Seki, K.; Mori, W. J. Phys. Chem. B 2002, 106, 1380. (11) Kitaura, R.; Iwahori, F.; Matsuda, R.; Kitagawa, S.; Kubota, Y.; Takata, M.; Kobayashi, T. C. Inorg. Chem. 2004, 43, 6522. (12) Keskin, S.; Liu, J.; Rankin, R. B.; Johnson, J. K.; Sholl, D. S. Ind. Eng. Chem. Res. 2009, 48, 2355. (13) Duren, T.; Bae, Y. S.; Snurr, R. Q. Chem. Soc. ReV. 2009, 38, 1237. (14) Tafipolsky, M.; Amirjalayer, S.; Schmid, R. Microporous Mesoporous Mater. 2010, 129, 304. (15) Demontis, P.; Suffritti, G. B. Chem. ReV. 1997, 97, 2845. (16) Smit, B.; Maesen, T. L. M. Chem. ReV. 2008, 108, 4125. (17) Ferey, G.; Serre, C. Chem. Soc. ReV. 2009, 38, 1380. (18) Tanaka, D.; Nakagawa, K.; Higuchi, M.; Horike, S.; Kubota, Y.; Kobayashi, L. C.; Takata, M.; Kitagawa, S. Angew. Chem., Int. Ed. 2008, 47, 3914. (19) In principle, for example, periodic plane wave DFT calculations can and are used to model the MOFs; see ref 14. However, with their large unit cells these computations are far from routine; in particular, the computation of a large number of configurations is still out of reach. (20) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. (21) Greathouse, J. A.; Allendorf, M. D. J. Phys. Chem. C 2008, 112, 5795. (22) Dubbeldam, D.; Walton, K. S.; Ellis, D. E.; Snurr, R. Q. Angew. Chem., Int. Ed. 2007, 46, 4496. (23) Han, S. S.; Goddard, W. A. J. Phys. Chem. C 2007, 111, 15185. (24) Huang, B. L.; McGaughey, A. J. H.; Kaviany, M. Int. J. Heat Mass Transfer 2007, 50, 393. (25) Tafipolsky, M.; Amirjalayer, S.; Schmid, R. J. Comput. Chem. 2007, 28, 1169. (26) Tafipolsky, M.; Schmid, R. J. Phys. Chem. B 2009, 113, 1341. (27) Salles, F.; Ghoufi, A.; Maurin, G.; Bell, R. G.; Mellot-Draznieks, C.; Ferey, G. Angew. Chem., Int. Ed. 2008, 47, 8487. (28) Bahr, D. F.; Reid, J. A.; Mook, W. M.; Bauer, C. A.; Stumpf, R.; Skulan, A. J.; Moody, N. R.; Simmons, B. A.; Shindel, M. M.; Allendorf, M. D. Phys. ReV. B 2007, 76, 184106. (29) El-Kaderi, H. M.; Hunt, J. R.; Mendoza-Cortes, J. L.; Cote, A. P.; Taylor, R. E.; O’Keeffe, M.; Yaghi, O. M. Science 2007, 316, 268. (30) Schmid, R.; Tafipolsky, M. J. Am. Chem. Soc. 2008, 130, 12600. (31) Klontzas, E.; Tylianakis, E.; Froudakis, G. E. Nano Lett. 2010, 10, 452. (32) Vishnyakov, A.; Ravikovitch, P. I.; Neimark, A. V.; Bulow, M.; Wang, Q. M. Nano Lett. 2003, 3, 713. (33) Skoulidas, A. I. J. Am. Chem. Soc. 2004, 126, 1356. (34) Liu, B.; Smit, B. Langmuir 2009, 25, 5918. (35) Liu, J. C.; Rankin, R. B.; Johnson, J. K. Mol. Simul. 2009, 35, 60. (36) Yazaydin, A. O.; Benin, A. I.; Faheem, S. A.; Jakubczak, P.; Low, J. J.; Willis, R. R.; Snurr, R. Q. Chem. Mater. 2009, 21, 1425. (37) Babarao, R.; Jiang, J. W.; Sandler, S. I. Langmuir 2009, 25, 5239.

Tafipolsky et al. (38) Castillo, J. M.; Vlugt, T. J. H.; Calero, S. J. Phys. Chem. C 2008, 112, 15934. (39) Chmelik, C.; Karger, J.; Wiebcke, M.; Caro, J.; van Baten, J. M.; Krishna, R. Microporous Mesoporous Mater. 2009, 117, 22. (40) Garcia-Perez, E.; Gascon, J.; Morales-Florez, V.; Castillo, J. M.; Kapteijn, F.; Calero, S. Langmuir 2009, 25, 1725. (41) Yang, Q.-Y.; Zhong, C.-L. J. Phys. Chem. B 2006, 110, 17776. (42) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (43) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (44) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (46) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (47) Tafipolsky, M.; Schmid, R. J. Chem. Theory Comput. 2009, 5, 2822. (48) Wilson, A. K.; Woon, D. E.; Peterson, K. A.; Dunning, T. H. J. Chem. Phys. 1999, 110, 7667. (49) Figgen, D.; Rauhut, G.; Dolg, M.; Stoll, H. Chem. Phys. 2005, 311, 227. (50) Peterson, K. A.; Puzzarini, C. Theor. Chem. Acc. 2005, 114, 283. (51) Feller, D. J. Comput. Chem. 1996, 17, 1571. (52) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. J. Chem. Inf. Model. 2007, 47, 1045. (53) Tinker - Software Tools for Molecular Design, version 4.2; Ponder, W.; Ren, P.; Pappu, R. V.; Hart, R. K.; Hodgson, M. E.; Cistola, D. P.; Kundrot, C. E.; Richards, F. M. Washington University School of Medicine, June 2004; http://dasher.wustl.edu/tinker/. (54) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (55) Cotton, F. A.; Murillo, C. A.; Walton, R. A. Multiple Bonds between Metal Atoms, 3rd ed.; Springer Science: New York, 2005. (56) Cotton, F. A.; Feng, X. J. J. Am. Chem. Soc. 1998, 120, 3387. (57) Melnik, M.; Kabesova, M.; Koman, M.; Macaskova, L.; Garaj, J.; Holloway, C. E.; Valent, A. J. Coord. Chem. 1998, 45, 147. (58) Bleaney, B.; Bowers, K. D. Proc. R. Soc. London, Ser. A 1952, 214, 451. (59) Van Niekerk, J. B.; Schoening, F. R. L. Acta Crystallogr. 1953, 6, 227. (60) de Loth, P.; Cassoux, P.; Daudey, J. P.; Malrieulb, J. P. J. Am. Chem. Soc. 1981, 103, 4007. (61) Rivero, P.; Moreira, I. d. P. R.; Illas, F.; Scuseria, G. E. J. Chem. Phys. 2008, 129, 184110. (62) Rodriguez-Fortea, A.; Alemany, P.; Alvarez, S.; Ruiz, E. Chem.sEur. J. 2001, 7, 627. (63) Ukai, T.; Nakata, K.; Yamanaka, S.; Takada, T.; Yamaguchi, K. Mol. Phys. 2007, 105, 2667. (64) Ali, M. E.; Datta, S. N. THEOCHEM 2006, 775, 19. (65) Noodleman, L. J. Chem. Phys. 1981, 74, 5737. (66) Noodleman, L.; Case, D. A. AdV. Inorg. Chem. 1992, 38, 423. (67) Elmali, A. Turk. J. Phys. 2000, 24, 667. (68) Perdew, J. P.; Ruzsinszky, A.; Constantin, L. A.; Sun, J.; Csonka, G. I. J. Chem. Theory Comput. 2009, 5, 902. (69) Brown, G. M.; Chidambaram, R. Acta Crystallogr., Sect. B 1973, 29, 2393. (70) Prestipino, C.; Regli, L.; Vitillo, J. G.; Bonino, F.; Damin, A.; Lamberti, C.; Zecchina, A.; Solari, P. L.; Kongshaug, K. O.; Bordiga, S. Chem. Mater. 2006, 18, 1337. (71) Wu, Y.; Kobayashi, A.; Halder, G. J.; Peterson, V. K.; Chapman, K. W.; Lock, N.; Southon, P. D.; Kepert, C. J. Angew. Chem., Int. Ed. 2008, 47, 8929. (72) Besler, B. H.; Merz, K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431. (73) Allinger, N. L.; Zhou, X. F.; Bergsma, J. THEOCHEM 1994, 118, 69. (74) Allinger, N. L.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. Soc. 1989, 111, 8551. (75) Amirjalayer, S.; Tafipolsky, M.; Schmid, R. Angew. Chem., Int. Ed. 2007, 46, 463.

Force Field for Cu Metal-Organic Frameworks (76) Piquemal, J.-P.; Gresh, N.; Giessner-Prettre, C. J. Phys. Chem. A 2003, 107, 10353. (77) Allured, V. S.; Kelly, C. M.; Landis, C. R. J. Am. Chem. Soc. 1991, 113, 1. (78) Rappe, A. K.; Bormann-Rochotte, L. M.; Wiser, D. C.; Hart, J. R.; Pietsch, M. A.; Casewit, C. J.; Skiff, W. M. Mol. Phys. 2007, 105, 301. (79) Allinger, N. L.; Li, F.; Yan, L. Q.; Tai, J. C. J. Comput. Chem. 1990, 11, 868. (80) Delgado-Friedrichs, O.; O’Keeffe, M.; Yaghi, O. M. Acta Crystallogr., Sect. A 2006, 62, 350. (81) Peterson, V. K.; Liu, Y.; Brown, C. M.; Kepert, C. J. J. Am. Chem. Soc. 2006, 128, 15578.

J. Phys. Chem. C, Vol. 114, No. 34, 2010 14409 (82) Rowsell, J. L. C.; Yaghi, O. M. J. Am. Chem. Soc. 2006, 128, 1304. (83) Zhang, X.; Fischer, R. A. Unpublished results. (84) The S-PXRD values plotted in Figure 5 are taken from Table S1 of the Supporting Information in ref 71. From this data, an averaged linear volumetric expansion coefficient of-1.06 × 10-5 K-1 can be derived by linear regression. (85) Chapman, K. W.; Halder, G. J.; Chupas, P. J. J. Am. Chem. Soc. 2008, 130, 10524. (86) Sun, D.; Ma, S.; Simmons, J. M.; Li, J.-R.; Yuan, D.; Zhou, H.-C. Chem. Commun. 2010, 46, 1329.

JP104441D