Mechanisms of Plastic Deformation of Metal–Organic Framework-5

Oct 28, 2015 - The Impact of Mesopores on the Mechanical Stability of HKUST-1: A Multiscale Investigation. Johannes P. Dürholt , Julian Keupp , and R...
0 downloads 3 Views 9MB Size
Download PDF
Article pubs.acs.org/JPCC

Mechanisms of Plastic Deformation of Metal−Organic Framework‑5 Kiettipong Banlusan, Edwin Antillon, and Alejandro Strachan* School of Materials Engineering and Birck Nanotechnology Center Purdue University, West Lafayette, Indiana 47907, United States S Supporting Information *

ABSTRACT: We use large-scale molecular dynamics simulations to investigate the mechanisms responsible for plastic deformation in metal−organic framework-5 (MOF-5). Simulations of uniaxial compression along [001], [101], and [111] directions reveal that structural collapse of {001} planes is responsible for irreversible deformation. The process involves slip along either one of the two ⟨100⟩ directions on the collapsing plane; this local shear process is due to the flexibility of the connection between of Zn−O clusters and 1,4-benzenedicarboxylate ligands. Thus, the collapse is driven both by compressive and shear stresses, and this fact explains the anisotropy in the mechanical response of this cubic crystal. The development of shear-collapse bands follows a nucleation and growth process with nuclei elongated along the slip direction and their subsequent growth in the directions normal to the slip and at much slower rates. This process is reminiscent of the glide of screw dislocations. Compression along the [101] and [111] directions led to intersection of active shear-collapse bands and the activation of multiple ⟨001⟩{100} systems. We also find that partially collapsed planes reduce the stiffness of the structures, an observation that can explain discrepancies between experimental and theoretical stiffness predictions.

1. INTRODUCTION Metal−organic frameworks (MOFs) are a relatively new class of materials that consists of short organic chains bridging between metal complexes; they exhibit an impressive diversity of chemistries and structures with very high porosity at the molecular scale. These characteristics make them attractive for a wide range of applications, including gas storage, sequestration, separation, catalysis, and biomedicine.1−8 Mechanical damping and shock absorption have also been proposed and investigated as possible applications for MOFs making use of the collapse of their nanoscale pores to dissipate mechanical energy.9 In addition to such applications, understanding the mechanical response of these materials to external loads is critical for their synthesis and processing where they are subjected to significant stresses.8 Despite the importance of these emerging materials, significant questions regarding their mechanical response remain unanswered. Specifically, the molecular-level mechanisms that govern their inelastic deformation are unknown10,11 beyond the processes of pore collapse at the unit cell level.12 From a fundamental science stand point uncovering the mechanisms of plastic deformation in MOFs is intriguing as their unique structures is likely to result in processes different from those known to operate in other classes of material, such as dislocation slip and twinning in metals or shear banding in amorphous systems. Our large-scale MD simulations confirm this to be the case, with inelastic deformation driven both by compressive and shear loads and showing a nucleation and propagation akin to that in dislocations but with important differences. © XXXX American Chemical Society

The elastic response of MOFs has been characterized experimentally via indentation11 and theoretically using density functional theory (DFT).11,13−17 Nanoindentation experiments resulted in stiffness of MOF-5 about an order of magnitude smaller than the theoretical predictions; the authors speculated that plastic deformation could account for the smaller experimental values.11 Such poorly accounted for discrepancies between theory and experiments even in the relatively simple case of small deformations attest to the importance of development of a predictive understanding of this class of materials. Both experiments and theory are beginning to shed light on the response of MOFs beyond the elastic limit. This process can be driven by mechanical loads and by gas adsorption/ desorption. The deformation of individual crystallites of zeolitic-imidazolate framework (ZIF-8) under compression was studied in situ under a TEM. The authors demonstrated very large strains and amorphization for deformations over 25%. Interestingly, the instantaneous Young’s modulus decreases after the initial deformation and then gradually increases when the process of amorphization starts.10 The pressure-induced amorphization of ZIF-8 was also studied using molecular dynamics (MD) simulations. The simulations show that ZIF-8 softens under hydrostatic compression due to the reduction of shear modulus.18 In addition, DFT and tightReceived: June 8, 2015 Revised: September 30, 2015

A

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

and c = [101]. This 848-atom unit cell was replicated 4 × 6 × 9 times resulting in a supercell with 183168 atoms and 15.26 × 16.19 × 34.34 nm3 dimensions. Similarly, simulations of compression along [111] start with a 636 atoms unit cell with cell unit vectors: a = [101̅], b = [1̅21̅], c = [111], which was replicated 8 × 4 × 7 times to result in a supercell containing 142,464 atoms with dimension of 15.26 × 13.21 × 32.71 nm3. All simulation cells are periodic in x, y, and z directions. 2.2. Simulations of Uniaxial Deformation. All MD simulations were performed using the LAMMPS code from Sandia National Laboratories.29 Prior to deformation, the simulation cells were equilibrated via isothermal−isobaric (NPT) simulations for 50 ps at 300 K and zero pressure. We used a Nosé−Hoover thermostat (with a coupling constant of 0.1 ps) and barostat (with a coupling constant of 5 ps). Uniaxial compression of the equilibrated systems was performed by continuously changing the dimensions of simulation box at every step of the MD simulation with an engineering strain rate of −3 × 109 s−1. The lateral dimensions of the cell were allowed to relax to maintain zero stress using a Nosé−Hoover barostat. The simulations were performed at 300 K and up to 150 ps, resulting in the maximum strain of 0.45. 2.3. Force Field. Atomic interactions are described using the reactive force field ReaxFF developed to describe the interaction between MOFs and water.30 This potential was parametrized using ab initio electronic structure calculations of cohesive energies and elastic constants of four ZnO crystal structures and their interactions with water31 combined with hydrocarbon parametrizarion.30 Predicted lattice constants and water stability of MOFs are in good agreement with experiments.30 An advantage of ReaxFF over nonreactive valence bond potential is its ability to describe chemical reactions involving bond breaking and formation. While ReaxFF has been shown to describe complex chemical processes,32−34 it is not without approximations. As with any MD study, the functional form and parameters used in the simulations affect atomistic trajectories and could influence results and conclusions. Thus, it is critical to establish whether the main results are robust with respect to the uncertainties in the interatomic potential used for the simulations. To accomplish this, we performed additional simulations with the consistent valence force field (CVFF) tuned to describe MOFs in ref 35. As shown in detail in section S2 of the Supporting Information (SI), both ReaxFF and CVFF lead to the same overall picture of plastic deformation, including nucleation and propagation of strain localization and the general shape of the stress−strain curves. In addition, the relaxation mechanism predicted by the force field calculations at the unit cell level was tested using density functional theory (DFT) calculations. These calculations, described in detail in section S1 of the SI, were performed using SeqQuest, a localized basis set DFT code developed at Sandia National Laboratories.36 The examples of LAMMPS input files for MD simulations with ReaxFF and CVFF force fields and SeqQuest input file for DFT calculations, described in detail in section S4 of the SI, are provided as separate supplementary files. These MD and DFT calculations can be performed online without downloading or installing any software using the LAMMPS37 and SeqQuest tools available in nanoHUB.38 Section S1 in the SI also compares the elastic constant predicted by ReaxFF and CVFF with available experiments and DFT predictions. 2.4. MD Analysis: Characterization of Strain Localization. In order to capture the mechanisms of deformation of

binding methods were used to study stacking faults in ZIF-8 showing a relatively high energy barrier to a low-energy intrinsic stacking fault; the authors suggest that this may explain the tendency to amorphization.19 MD simulations also revealed the process responsible for void collapse under hydrostatic compression in several MOFs at the unit cell level.12 Collapse is accompanied by rearrangement of the linker molecules and by shear. MOFs can undergo structural deformation during adsorption/desorption due to their interaction with the guest molecules. For example, reversible transitions between largepore and narrow-pore states have been experimentally observed in the flexible structures of chromium20 and vanadium21 benzenedicarboxylates. Monte Carlo simulations revealed that this transition is driven by the shear on the layers of BDC segments.22,23 Irreversible structural collapse was reported in MOF-5 for a critical amount of water adsorption.24 While understanding processes at the unit cell is a key first step, mechanical deformation beyond the elastic limit rarely occurs homogeneously. For example, metals localize plastic deformation on slip planes via the glide of dislocations,25,26 and many amorphous solids27 and molecular crystals28 localize strain in shear bands. The mechanisms by which inelastic deformation in porous MOFs nucleates and propagates have not been characterized, and they are the focus of this paper. Our large-scale MD simulations of uniaxial deformation of MOF-5 show that plastic deformation localizes among easily collapsible planes, but the nucleation event involves the collapse of one or few unit cells and not of the entire plane.

2. SIMULATION DETAILS 2.1. Systems of Interest. The cubic unit cell MOF-5, illustrated in Figure 1b, consists of tetrahedral Zn4O clusters

Figure 1. (a) Stress−strain relations of MOF-5 loaded along [001], [101], and [111] directions. Abrupt reduction of stress (e.g., at strain of 2.5%−5% for [001] loading) corresponds to the structural collapse of MOF-5 during the plastic deformation. (b) Structure of a cubic unit cell of MOF-5 consisting of Zn−O clusters connected by organic BDC ligands with flexible bonds. (c) Deformation of MOF-5 unit cell by bending of the flexible connectors between Zn−O clusters and organic linkers upon compression along [001] direction.

connected by 1,4-benzenedicarboxylate (BDC) ligands. The chemical formula of this material is [Zn4O(BDC)3]8 with 424 atoms per unit cell. To gain insight into the anisotropy in mechanical response and deformation mechanisms of this material, we performed uniaxial compression simulations along the [001], [101], and [111] directions. For compression along [001], the simulation cell was prepared by replicating the cubic unit cell 6 × 6 × 12 times along x, y, and z directions resulting in 183168 atoms in a supercell with dimension of 16.19 × 16.19 × 32.38 nm3. For compression along [101], a unit cell with the following cell unit vectors was created: a = [101]̅ , b = [010], B

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Table 1. Schmid Factors and Compression Factors of Possible Active Slip-Collapse Systems of MOF-5 Loaded along Different Directions slip systems

loading direction [001]

loading direction [101]

loading direction [111]

slip plane

slip direction

Schmid

compress.

Schmid

compress.

Schmid

compress.

(001)

[100] [010] [100] [001] [010] [001]

0 0 0 0 0 0

1 1 0 0 0 0

1/2 0 0 0 0 1/2

1/√2 1/√2 0 0 1/√2 1/√2

1/3 1/3 1/3 1/3 1/3 1/3

1/√3 1/√3 1/√3 1/√3 1/√3 1/√3

(010) (100)

Figure 2. Atomistic snapshots show plastic deformation of MOF-5 during the compression along [001] direction at various strains corresponding to the values identified by the letters in Figure 1a. Colors indicate the percentage of local volume decrease compared to the initial value in asequilibrated systems. For visualizing clarification, the color scale bar on the left varies from 0% to 30%; however, the volume of any groups of atoms can reduce by more than 30%. The bottom panels show only groups of atoms with volume decrease more than 10% to clarify the nucleation and propagation of plasticity. The arrows in the bottom panel of Figure 2d illustrate the propagation in ⟨100⟩ and ⟨010⟩ directions on (001) plane of volume collapse.

level is shown in Figure 1c and involves collapse of {001} planes in the structure with a relative displacement between consecutive plans along one of the two ⟨100⟩ directions contained on the collapsing plane. In the case of loading along [001], the driving force for the structural collapse is the compressive stress on the (001) plane, σzz. Once a threshold value is reached, the collapse is abrupt. In the case of compression along [101] and [111] the active plane normals make angles of 45° and 54.74°, respectively, with the loading axis. Thus, both compressive and shear stresses act on the active plane. This shear stress facilitates the deformation of the unit cell shown in Figure 1c and results in a lower yield stress and a more gradual softening. In contrast to the behavior of MOFs, plastic deformation in fcc metals is dominated by the resolved shear stress on the (111) slip planes, and other components of the stress play a secondary role.25 Table 1 summarizes the fraction of the applied stress that results in compression normal to the active planes and the Schmid factors (the shear stress resolved along slip direction) for all slipcollapse systems and the three loading conditions studied here (see section S3 in the SI for more detail). The key questions that remain to be answered are whether the volume collapse responsible for deformation happens homogeneously or there is strain localization and how these processes nucleate and propagate throughout the sample.

MOF-5 we need to discriminate between local elastic and plastic processes and quantify the local plastic deformation. The perfect structure, see Figure 1b, is a simple cubic arrangement of Zn4O, each connected to its six neighbors via BDC segments. Thus, we compute local volume changes by tracking the evolution of this simple cubic lattice; each lattice site is obtained as the center of mass position of each Zn4O group and the nearest half of the six BDC segments of BDC. The volume associated with each lattice site is calculated by grouping its six nearest neighbors into three pairs of opposing sites. The resulting vectors u⇀1,u⇀2,u⇀3 can be viewed as the edges of parallelepiped surrounding the central site, and its volume can be calculated as V = |u⇀1·(u⇀2 × u⇀3)|.

3. STRESS−STRAIN RELATIONSHIPS AND YIELD CONDITIONS Figure 1a shows stress−strain curves for the MOF-5 samples loaded along the [001], [101], and [111] directions. The Young’s modulus along [001] is significantly higher than those along [101] and [111]. This can be understood in terms of the molecular structure of MOF-5; see Figure 1b. Loading along ⟨001⟩ directions results in the compression of the organic BDC segments, while deformation along the other two directions can be accommodated by, more compliant, bending of the structure. This also explains the more abrupt reduction in the stress when the yield point is reached. The molecular-level mechanism responsible for plastic deformation at the unit cell C

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 3. Snapshots of 2D layers representing the C.M. position of Zn4O[BDC]3 clusters in four consecutive (001) planes at the region indicated by a dashed-line square in Figure 2d for MOF-5 loaded along the [001] direction. The left panels show the side view (YZ plane) of analyzed layers. The middle and right panels show the top views (XY plane) of C/D and A/B planes, respectively. Solid arrows indicate the slip direction, and broken arrows represent the propagation direction of the plastic deformation.

4. MOLECULAR MECHANISMS OF PLASTIC DEFORMATION 4.1. Compression along [001] Direction. Figure 2 shows atomistic snapshots during the compression of MOF-5 along the [001] direction; the letters correspond to the strains marked in Figure 1a. For clarity, the bottom panels display only regions that have undergone volume collapse over 10%. Atoms are colored based on the local volume change defined in the simulation details section. As expected, the collapse occurs on (001) planes. We observe significant strain localization and, interestingly, interactions between regions of plastic deformation. Volume collapse at a given location releases the compressive stress in its local neighborhood and can result in a shielding effect at distances of several unit cells away from collapsed region. This can be seen in several locations in snapshot b corresponding to a strain of 3.75% and at later times. The propagation of the volume collapse in a given (001) plane releases the compressive stresses in neighboring planes and can even lead to the recovery of plasticity of those regions as illustrated in parts c−e of Figure 2. This phenomenon is even more apparent at lower strain rates and for larger simulation cells; see section S2 in the SI. The initial nucleation of local slip-collapse events within a defect-free material requires significant levels of stress (130 MPa from the ReaxFF simulation and 170 MPa according to CVFF; see section S2 in the SI) and leads to a significant relaxation in stress. The large relaxation and jagged stress− strain curves are common in small systems since a single plastic even represents a significant fraction of the total strain; such behavior has been observed in metallic nanowires and pillars.39 The subsequent nucleation events require much smaller stress, see points f to g in Figures 1a and 2 (see also Figures S9 and S11 in the SI), being facilitated by the heterogeneous nature of the material containing collapsed regions. Note that the

nucleation of a new void collapse in Figure 2g occurs in a plane right above a previously collapsed one. Nucleation and propagation processes. To uncover the molecular processes responsible for nucleation and propagation of localized plastic deformation, we follow the temporal evolution of four consecutive planes during their collapse in Figure 3. The left panels provide a side view where void collapse in four subsequent planes denoted A, B, C, and D. The top views on the right panels provide additional information about the nucleation and propagation process. At 3.75% strain, i.e., shortly after the yield point, the figure shows the nucleation of two slip-collapse bands; in both cases, the slip is along [100] and the bands are elongated along the slip direction and approximately two and a half unit cells in thickness. It is interesting to note that the nucleation of the slip-collapse bands occurs at opposite locations in the sample; this is due to shielding as discussed above. In both cases, a narrow strip of material undergoes slip with minimal volume collapse; the collapse trails the slip as the material continues to be deformed. At 5.63% strain local volume collapse is clear in the regions that have undergone slip and deformation bands propagate in the two directions normal to the slip; i.e., propagation is along ±[010]. Dashed lines separate the regions where plastic deformation has occurred from the intact material and help visualize the propagation of the deformation band. In the case of loading along [001] there are four active slip-collapse systems [010](001), [01̅0](001), [100](001), and [1̅00](001). All of these are equivalent and are seen in the simulations. The picture emerging from these results is very reminiscent of a screw dislocation, where the slip direction is parallel to the dislocation line that also separates slipped and nonslipped material. Two important differences should be noted. First, contrary to the slip-collapse bands, dislocation slip does not lead to volume changes and is, thus, not driven by normal D

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 4. Atomistic snapshots show plastic deformation of MOF-5 during compression along the [101] direction at various strains corresponding to the values identified by the letters in Figure 1a. The colors represent the local volume decrease with the color scale bar on the left varying from 0% to 35% for visualizing clarification; however, volume decreases of more than 35% can occur. The clusters of atoms with volume decrease more than 15% are illustrated in the bottom panels showing that the plastic deformation occurs on (001) and (100) planes.

stresses. Second, multiple dislocation slip events can occur on the same slip plane (this results in macroscopic slip steps on the surfaces of single crystals), while in the case of MOF-5 collapse can occur only once and additional slips would require breaking the MOF-5 network; this is not observed in our simulations. Void Collapse and Stiffness. These atomistic-level details can shed light into the discrepancies between nanoindentation measurements of stiffness for MOF-5 (2.7 GPa) and the much higher DFT predictions (∼21 GPa). Plastic deformation was proposed as a likely reason for the discrepancies,11 and our simulations support this hypothesis. After the initial collapse compressing the plastically deformed material is more compliant. For example, the effective Young’s moduli from points c to d and from f to g are 2.3 and 1.9 GPa, respectively, lower than the perfect crystal value of 6.7 GPa. This is because partially collapsed planes exhibit lower stiffness. The SI (Figure S10) contains additional data for a larger sample compressed at lower rates. 4.2. Compression along [101] Direction. As opposed to compression along [001], loading along [101] leads both to resolved compressive stress on the active {100} planes and also a resolved shear stress on the {100}⟨ 001⟩ systems. In this case, four slip-collapse systems exhibit compressive or shear loads: [100](001), [010](001), [010](100), and [001](100). However, only two of these, [100](001) and [001](100), are driven both by compressive stress and resolved shear stress, whereas [010](001) and [010](100) only experience compression; see Table 1. The (010) planes, whose normal is 90° from the load axis, have zero resolved shear and compressive stress and never becomes active. Figures 4 and 5 provide insight into the molecular processes of plastic deformation during compression along [101] direction. As expected, slip-collapse bands nucleate and propagate on both systems, [100](001) and [001](100), exhibiting compressive and shear loading. The intersection between slip-collapse bands results in 1D regions of large volume decrease, as indicated by red in Figure 4d−f. We also observe that for large deformations and after the initial microstructure with collapsed regions develops plastic deformation is likely to nucleate near the predeformed regions similar to the previous case, resulting in growth of plastically deformed regions as shown in Figure 4g,h.

Figure 5. Snapshots of 2D layers representing the C.M. position of Zn4O[BDC]3 clusters in two consecutive (100) planes (A/B planes) and two consecutive (001) planes (C/D planes) at the region indicated by a dashed-line square in Figure 4e for compression of MOF-5 along [101] direction. (a) The top panel shows four analyzed layers subjected to the compressive stress along z direction, resulting in resolved shear stresses on (100) and (001) planes. The middle and bottom panels show C/D and A/B layers viewed along the normals of (001) and (100) planes, respectively, with the intersection between the planes indicated by red/green lines for C/D layers, and pink/blue lines for A/B layers. (b) At strain of 22.5%, the plastic deformation is nucleated by slips along [001] direction on (100) plane (A/B layers) and [100] direction on (001) plane (C/D layers). Slip directions are indicated by solid-line arrows, while dashed-line arrows show the propagation directions, which are perpendicular to the slip directions. (c) Slip along [010] direction that is driven by the compressive stress normal to the slip planes occurs at the higher strain (33.75%).

Figure 5 shows the nucleation and propagation of slipcollapse on two sets of intersecting planes located in the region surrounded by dashed lines in Figure 4e. This analysis shows a more complex process than in the [001] case. Initially, the E

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 6. Atomistic snapshots show plastic deformation of MOF-5 during compression along the [111] direction at various strains corresponding to the values identified by the letters in Figure 1a. The colors represent the local volume decrease with the color scale bar on the left varying from 0% to 70% for visualizing clarification; however, volume decreases of more than 70% can occur. The groups of atoms with volume decrease of more than 35% are illustrated in the bottom panels showing that the plastic deformation occurs on {001} family planes.

Figure 7. Snapshots of 2D layers representing the C.M. position of Zn4O[BDC]3 groups in two consecutive {001} planes, viewed along their normal, at the region indicated in the bottom left panel for compression of MOF-5 along the [111] direction. The colors represent the local volume decrease consistent with Figure 6. Two slip directions can be activated as illustrated by the solid-line arrows.

structure while the compressive stresses on (001) and (100) facets increase. 4.3. Compression along [111] Direction. Contrary to the previous two cases, the six possible ⟨100⟩{001} slip-collapse systems experience the same level of compression and shear during compression along the [111] direction, as indicated by Table 1. Figure 6 shows strain localization of MOF-5 during compression along [111]. All three {001} planes become active. The higher ratio between compression and Schmid factors of the slip-collapse systems compared to the [101] case (see Table S3 and S4 in the SI) suggests that processes are more dominated by compressive stresses. Thus, compression along the [111] direction exhibits higher yield stress (54 MPa) than that of the [101] direction (38 MPa) as shown by stress−strain curves in Figure 1a. For each slip-collapse plane, the two possible slip directions exhibit identical shear stresses and, consequently, can occur with equal probability. We observe the coexistence of both

angle between the slip plane normals ((001) and (100) planes) and applied stress is 45°, resulting in a high resolved shear stress. Thus, [100](001) and [001](100) slip systems are preferably activated (see Figure 5b at strain = 22.5%). Unlike the [001] case, the active planes intersect as shown in Figure 5a and the two slip directions are not consistent with each other; both slip systems cannot coexist at the intersection without breaking the structure. Thus, we observe slip along the [010] direction; see Figure 5c at a strain of 33.75%. The [010](001) and [010](100) slip systems that do not experience resolved shear stresses and are not activated initially become active. We note that for large deformations the latter slip systems also become more attractive due to the significant deformation of the structure. The initial slip systems, [100](001) and [001] (100), make 45° angles with the compressive axis, leading the highest possible resolved shear stress. However, at high strains the angle between the slip plane normal and applied stress decreases significantly because of the bending of MOF-5 F

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C possible slip options, as can be seen in Figure 7. The figure shows the nucleation and propagation of slip-collapse on {001} plane of MOF-5 structure loaded along [111] direction. Two perpendicular slip directions with the same driving forces are activated on this slip plane. Groups of atoms at the intersection between the active slip planes can experience another perpendicular driving forces along the intersected planes, resulting in more volume decrease as indicated by red color at strain of 21.1% in Figure 7. The slip occurs along the first slip direction at strain of 21.1−22.9%. The second slip step is activated at higher strain, and it ends up with a collapsed structure with a relative displacement between two consecutive layers along the combination of two slip direction as shown in Figure 7 at strain of 34.5%.

cal predictions of the stiffness of this material, while the simulations focused on MOF-5, we believe, are applicable to a wide range of MOFs and ZIFs that share similar topologies. Finally, the mechanistic understanding of inelastic deformation presented here may enable the design of new MOFs with tailored mechanical properties.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b05446. Force field validation, elastic constant calculations, DFT calculations of structural collapse, additional simulations using CVFF force field, additional simulations of the larger systems and lower strain rate, calculations of the compressive stress and resolved shear stress factors (PDF) Input files (ZIP)

5. DISCUSSION AND CONCLUSIONS In summary, large-scale MD simulations revealed the mechanisms responsible for the inelastic deformation of MOFs; we find the process to be mediated by the nucleation and propagation of shear-collapse bands. These bands form by bending of flexible bonds between Zn−O clusters and BDC ligands, leading to the collapse of the extensive free volume in MOF-5 and an abrupt drop in stress of the system. As with any molecular dynamics simulation, our results depend on the description of atomic interactions. Thus, to confirm the mechanisms observed in the ReaxFF we repeated all simulations with a different, nonreactive force field also developed to study MOFs; the results, discussed in detail in section S2 of the SI, show identical processes. Additional simulations with larger system sizes and different strain rates (section S2 of the SI) establish the robustness of our conclusions with respect to simulation details. Finally, while higher accuracy ab initio calculations of such large systems are unfeasible, we confirmed the relaxation mechanisms at the unit cell level using density functional theory. We performed uniaxial compression simulation using ReaxFF on a single unit cell and took several snapshots from this process as input structures for density functional theory (DFT) calculations. We first relax the atomic structure via energy minimization maintaining the cell parameters fixed; these simulations show that the shear-collapse processes predictive by ReaxFF relaxed the uniaxial stress. Subsequent DFT simulations relaxed also the cell parameters under uniaxial and zero applied external stress; the cell parameters show little change with respect to the ReaxFF predictions and show that the sheared-collapsed structure are metastable and correspond to local minima of the enthalpy. These calculations are described in section S1 of the SI. Collapse of (001) planes involves shear along the two [100] directions contained on the plane and is driven both by compressive and resolved shear stresses. The interplay between shear and compressive driving forces explains the significant anisotropy in the yielding of these materials. Shear-collapse bands follow a nucleation and growth process. Nuclei are elongated along the slip direction, and growth on the collapsing plane occurs in the direction normal to the slip direction. Active planes during compression along [001] are parallel to each other; this results in relatively weak interaction between shearcollapse bands. However, intersection of active shear-collapse bands during compression along [101] and [111] results in the activation of a large number of slip-collapse systems and their interaction. The results show that partially collapsed planes could explain discrepancies between experimental and theoreti-

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by the US Department of Defense, Office of Naval Research, MURI Program Grant 2012 02341 01. Program Manager: Kenny Lipkowitz. Insightful discussions with K. Suslick are gratefully acknowledged and so are computational resources from nanoHUB.org and Purdue.

■

REFERENCES

(1) Yaghi, O. M.; Li, H.; Eddaoudi, M.; O’Keeffe, M. Design and Synthesis of an Exceptionally Stable and Highly Porous Metal-Organic Framework. Nature 1999, 402, 276−279. (2) Eddaoudi, M.; Kim, J.; Rosi, N.; Vodak, D. T.; Wachter, J.; O’Keeffe, M.; Yaghi, O. M. Systematic Design of Pore Size and Functionality in Isoreticular MOFs and Their Application in Methane Storage. Science 2002, 295, 469−472. (3) Rosi, N.; Eckert, J.; Eddaoudi, M.; Vodak, D. T.; Kim, J.; O’Keeffe, M.; Yaghi, O. M. Hydrogen Storage in Microporous MetalOrganic Frameworks. Science 2003, 300, 1127−1129. (4) Rowsell, J. L. C.; Spencer, E. C.; Eckert, J.; Howard, J. A. K.; Yaghi, O. M. Gas Adsorption Sites in a Large-Pore Metal-Organic Framework. Science 2005, 309, 1350−1354. (5) Lee, J. Y.; Farha, O. K.; Roberts, J.; Scheidt, K. A.; Nguyen, S. T.; Hupp, J. T. Metal−Organic Framework Materials as Catalysts. Chem. Soc. Rev. 2009, 38, 1450−1459. (6) McKinlay, A. C.; Morris, R. E.; Horcajada, P.; Ferey, G.; Gref, R.; Couvreur, P.; Serre, C. BioMOFs: Metal−Organic Frameworks for Biological and Medical Applications. Angew. Chem., Int. Ed. 2010, 49, 6260−6266. (7) Furukawa, H.; Cordova, K. E.; O’Keeffe, M.; Yaghi, O. M. The Chemistry and Applications of Metal-Organic Frameworks. Science 2013, 341, 1230444. (8) Czaja, A. U.; Trukhan, N.; Muller, U. Industrial Applications of Metal−Organic Frameworks. Chem. Soc. Rev. 2009, 38, 1284−1293. (9) Yot, P. G.; Boudene, Z.; Macia, J.; Granier, D.; Vanduyfhuys, L.; Verstraelen, T.; Speybroeck, V. V.; Devic, T.; Serre, C.; Ferey, G.; Stockd, N.; Maurin, G. Metal−Organic Frameworks as Potential Shock Absorbers: The Case of The Highly Flexible MIL-53(Al). Chem. Commun. 2014, 50, 9462−9464. (10) Su, Z.; Miao, Y. R.; Mao, S. M.; Zhang, G. H.; Dillon, S.; Miller, J. T.; Suslick, K. S. Compression-Induced Deformation of Individual G

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C Metal−Organic Framework Microcrystals. J. Am. Chem. Soc. 2015, 137, 1750−1753. (11) Bahr, D.; Reid, J.; Mook, W.; Bauer, C.; Stumpf, R.; Skulan, A.; Moody, N.; Simmons, B.; Shindel, M.; Allendorf, M. Mechanical Properties of Cubic Zinc Carboxylate IRMOF-1 Metal-Organic Framework Crystals. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 184106−1−7. (12) Biswas, M. M.; Cagin, T. High Pressure Structural Transformation of Selected Metal Organic Frameworks-A Theoretical Investigation. Mater. Chem. Phys. 2011, 131, 44−51. (13) Zhou, W.; Yildirim, T. Lattice Dynamics of Metal-Organic Frameworks: Neutron Inelastic Scattering and First-Principles Calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 180301. (14) Mattesini, M.; Soler, J. M.; Ynduráin, F. Ab initio Study of MetalOrganic Framework-5 Zn4O(1,4-benzenedicarboxylate)3: An Assessment of Mechanical and Spectroscopic Properties. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 094111. (15) Samanta, A.; Furuta, T.; Li, J. Theoretical Assessment of the Elastic Constants and Hydrogen Storage Capacity of Some MetalOrganic Framework Materials. J. Chem. Phys. 2006, 125, 084714−1−8. (16) Ortiz, A. U.; Boutin, A.; Fuchs, A. H.; Coudert, F. X. Anisotropic Elastic Properties of Flexible Metal-Organic Frameworks: How Soft are Soft Porous Crystals? Phys. Rev. Lett. 2012, 109, 195502. (17) Wu, H.; Yildirim, T.; Zhou, W. Exceptional Mechanical Stability of Highly Porous Zirconium Metal−Organic Framework UiO-66 and Its Important Implications. J. Phys. Chem. Lett. 2013, 4, 925−930. (18) Ortiz, A. U.; Boutin, A.; Fuchs, A. H.; Coudert, F.-X. Investigating the Pressure-Induced Amorphization of Zeolitic Imidazolate Framework ZIF-8: Mechanical Instability Due to Shear Mode Softening. J. Phys. Chem. Lett. 2013, 4, 1861−1865. (19) Hegde, V. I.; Tan, J.-C.; Waghmare, U. V.; Cheetham, A. K. Stacking Faults and Mechanical Behavior beyond the Elastic Limit of an Imidazole-Based Metal Organic Framework: ZIF-8. J. Phys. Chem. Lett. 2013, 4, 3377−3381. (20) Serre, C.; Millange, F.; Thouvenot, C.; Nogues, M.; Marsolier, G.; Louer, D.; Ferey, G. Very Large Breathing Effect in the First Nanoporous Chromium(III)-Based Solids: MIL-53 or CrIII(OH) •{O2C-C6H4-CO2}•{HO2C-C6H4-CO2H}x•H2Oy. J. Am. Chem. Soc. 2002, 124, 13519−13526. (21) Wang, X.; Eckert, J.; Liu, L.; Jacobson, A. J. Breathing and Twisting: An Investigation of Framework Deformation and Guest Packing in Single Crystals of a Microporous Vanadium Benzenedicarboxylate. Inorg. Chem. 2011, 50, 2028−2036. (22) Triguero, C.; Coudert, F. X.; Boutin, A.; Fuchs, A. H.; Neimark, A. V. Mechanism of Breathing Transitions in Metal-Organic Frameworks. J. Phys. Chem. Lett. 2011, 2, 2033−2037. (23) Coudert, F. X.; Boutin, A.; Fuchs, A. H.; Neimark, A. V. Adsorption Deformation and Structural Transitions in Metal-Organic Frameworks: From the Unit Cell to the Crystal. J. Phys. Chem. Lett. 2013, 4, 3198−3205. (24) Schrock, K.; Schroder, F.; Heyden, M.; Fischer, R. A.; Havenith, M. Characterization of Interfacial Water in MOF-5 (Zn4(O) (BDC)3)-A Combined Spectroscopic and Theoretical Study. Phys. Chem. Chem. Phys. 2008, 10, 4732−4739. (25) Hirth, J. P.; Lothe, J. Theory of Dislocations; McGraw-Hill: New York, 1968. (26) Lee, D. W.; Kim, H.; Strachan, A.; Koslowski, M. Effect of Core Energy on Mobility in a Continuum Dislocation Model. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 104101. (27) Hays, C. C.; Kim, C. P.; Johnson, W. L. Microstructure Controlled Shear Band Pattern Formation and Enhanced Plasticity of Bulk Metallic Glasses Containing in situ Formed Ductile Phase Dendrite Dispersions. Phys. Rev. Lett. 2000, 84, 2901−2904. (28) Jaramillo, E.; Sewell, T. D.; Strachan, A. Atomic-Level View of Inelastic Deformation in a Shock Loaded Molecular Crystal. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 064112. (29) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19.

(30) Han, S. S.; Choi, S. H.; van Duin, A. C. T. Molecular Dynamics Simulations of Stability of Metal−Organic Frameworks Against H2O Using the ReaxFF Reactive Force Field. Chem. Commun. 2010, 46, 5713−5715. (31) Raymand, D.; van Duin, A. C. T.; Spangberg, D.; Goddard, W. A., III; Hermansson, K. Water Adsorption on Stepped ZnO Surfaces From MD Simulation. Surf. Sci. 2010, 604, 741−752. (32) Wood, M. A.; van Duin, A. C. T.; Strachan, A. Coupled Thermal and Electromagnetic Induced Decomposition in the Molecular Explosive αHMX; A Reactive Molecular Dynamics Study. J. Phys. Chem. A 2014, 118, 885−895. (33) Han, S. P.; van Duin, A. C. T.; Goddard, W. A., III; Strachan, A. Thermal Decomposition of Condensed-Phase Nitromethane from Molecular Dynamics from ReaxFF Reactive Dynamics. J. Phys. Chem. B 2011, 115, 6534−6540. (34) Strachan, A.; Kober, E. M.; van Duin, A. C. T.; Oxgaard, J.; Goddard, W. A., III Thermal Decomposition of RDX From Reactive Molecular Dynamics. J. Chem. Phys. 2005, 122, 54502. (35) Greathouse, J. A.; Allendorf, M. D. Force Field Validation for Molecular Dynamics Simulations of IRMOF-1 and Other Isoreticular Zinc Carboxylate Coordination Polymers. J. Phys. Chem. C 2008, 112, 5795−5802. (36) Schultz, P. A. SeqQuest code, http:/dft.sandia.gov/Quest. (37) Haley B. P. LAMMPS 2015, https://nanohub.org/resources/ lammps. DOI: 10.4231/D3S17ST5Z. (38) Vedula, R. P. K.; Bechtol, G.; Haley, B. P.; Strachan A. nanoMATERIALS SeqQuest DFT 2014, https://nanohub.org/ resources/nmst_dft. DOI: 10.4231/D3FQ9Q61P. (39) Uchic, M. D.; Dimiduk, D. M.; Florando, J. N.; Nix, W. D. Sample Dimensions Influence Strength and Crystal Plasticity. Science 2004, 305, 986−989.

H

DOI: 10.1021/acs.jpcc.5b05446 J. Phys. Chem. C XXXX, XXX, XXX−XXX