Probing the Soft and Nanoductile Mechanical Nature of Single and

Dec 2, 2016 - Department of Mechanical Engineering, Materials Science and Engineering Program, FAMU-FSU College of Engineering, Florida State Universi...
2 downloads 11 Views 13MB Size
Probing the Soft and Nanoductile Mechanical Nature of Single and Polycrystalline Organic− Inorganic Hybrid Perovskites for Flexible Functional Devices Jingui Yu,† Mingchao Wang,† and Shangchao Lin* Department of Mechanical Engineering, Materials Science and Engineering Program, FAMU-FSU College of Engineering, Florida State University, Tallahassee, Florida 32310, United States S Supporting Information *

ABSTRACT: Although organic−inorganic hybrid perovskites have been extensively investigated for promising applications in energy-related devices, their mechanical properties, which restrict their practical deployment as flexible and wearable devices, have been largely unexplored at the atomistic level. Toward this level of understanding, we predict the elastic constant matrix and various elastic properties of CH3NH3PbI3 (MAPbI3) using atomistic simulations. We find that singlecrystalline MAPbI3 is much stiffer and exhibits higher ultimate tensile strength than polycrystalline samples, but the latter exhibit unexpected, greatly enhanced nanoductility and fracture toughness, resulting from the extensive amorphization during the yielding process. More interestingly, polycrystalline MAPbI3 exhibits inverse Hall−Petch grain-boundary strengthening effect, in which the yield stress is reduced when decreasing the grain size, due to amorphous grain boundaries. By monitoring the centro-symmetry parameter and local stress evolution, we confirm the soft and nanoductile nature of defective MAPbI3 with a crack. By conducting atomic stress decomposition, we attribute such fracture toughness primarily to the strong electrostatic interactions between the ionic components. The observed limited brittle fracture behavior is attributed to the transformation of partial edge dislocations to disordered atoms (nanovoid formation). A significant plastic deformation region is observed when nanovoids enlarge and coalesce with adjacent ones, which ultimately leads to crack propagations via ionic-chain breaking. After comparing with traditional inorganic energy-related materials, we find that hybrid perovskites are more compressible and can absorb more strain energy before fracture, which makes them well suited for wearable functional devices with high mechanical flexibility and robustness. KEYWORDS: hybrid perovskite, polycrystalline thin films, fracture mechanics, crack propagation, microstructural evolution, molecular dynamics

V

organic−inorganic hybrid perovskites for PV and LED applications, it is essential to advance the existing understanding of their degradation pathways and mechanical stabilities, in particular, fracture mechanics considering cracks or defects, beyond the ideal picture of elasticity and plasticity for pristine samples. It is well-known that as-synthesized crystalline materials always possess invisible cracks and defects at the microscopic level. The initiation and propagation of these cracks or defects dictate the material’s failure and lifetime,

ery recently, organic−inorganic hybrid perovskites, such as the prototypical methylammonium lead iodide (CH3NH3PbI3 or MAPbI3), have emerged as promising energy-related materials, such as light absorbers in photovoltaic (PV) cells, emitters in light-emitting diodes (LEDs),1 or potential thermoelectric materials with high figure of merit.2,3 They provide the possibility of integrating useful organic and inorganic characteristics within a single crystalline molecularscale hybrid, enabling unique electronic,4 optical,5 thermoelectric,2,3 and surface6,7properties. These ambipolar semiconductors have attracted increasing attention due to their scalable, low-temperature and low-cost synthesis processes, small and direct electron bandgaps, high extinction coefficients, and high carrier mobility.8,9 With the growing interest in using © 2016 American Chemical Society

Received: September 1, 2016 Accepted: December 2, 2016 Published: December 2, 2016 11044

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

www.acsnano.org

Article

ACS Nano

Table 1. Comparison between Elastic Properties of Pristine (except for the Last Row for the Defective MD Model) MAPbI3 Obtained Using MM (from This Work), MD (from This Work), DFT19 Calculations, and Available Experimental NanoIndentation Results:17,18 Bulk Modulus (K), Shear Modulus (G), Young’s Modulus (E), and Poisson’s Ratio (ν)a MAPbI3 phases (methods) 19

orthorhombic (DFT) orthorhombic (MM)b tetragonal (DFT)19 tetragonal (MM)b tetragonal (experiment) cubic (DFT)19 cubic (MM)b single-crystalline cubic (MD)b polycrystalline cubic (MD)b defective single-crystalline cubic (MD)b

KVRH (GPa)

GVRH (GPa)

EVRH (GPa)

νVRH

18.1 17.80 12.2 20.25 13.918 16.4 18.08 13.11 2.51−4.87 12.32

3.6 9.27 3.7 10.51 5.418 8.7 11.91 5.70 2.15−3.21 5.36

15.0 23.68 12.8 26.87 10.4−10.7,17 14.0−14.318 22.2 29.29 14.94 5.08−7.89 14.04

0.36 0.28 0.33 0.28 0.3319 0.28 0.23 0.31 0.16−0.23 0.31

The subscript “VRH” denotes the Voigt−Reuss−Hill approximation, which is applied for MM and DFT calculations only. For MD (predicted both E E E and ν) and experimental (measured E and used ν from DFT) results, K and G were estimated based on K = 3(1 − 2ν) and G = 2(1 + ν) , assuming

a

that MAPbI3 is isotropic and homogeneous. bThis work.

their fracture mechanical properties and failure mechanism are still missing. Current experimental techniques have difficulty in obtaining high-quality bulk perovskite samples to probe crack propagation, while DFT calculations cannot handle length scales larger than 10 nm and time scales longer than 10 ns, both are relevant for probing fracture mechanics. Therefore, largescale MD simulations are well-suited for studying the fracture mechanics of hybrid perovskites to understand their failure mechanisms and predict the fracture toughness as measured by the critical strain energy release rate and the critical stress intensity factor for crack propagation. With the above in mind, we first computed the elastic constant matrix of pristine MAPbI3 using force field-based molecular mechanics (MM) calculations, for the cubic, tetragonal, and orthorhombic phases. We used the firstprinciple derived classical model potential for hybrid perovskites (MYP), recently developed by Mattoni et al.20,21 We utilized these elastic constants to determine the elastic properties of MAPbI3 and compared with our MD simulation results for pristine single and polycrystalline MAPbI 3, considering the finite temperature effect, together with previous DFT and experimental results. Then we carefully explored the fracture mechanical behavior and microstructural evolution of defective single-crystalline MAPbI3 (with a crack) under tensile loading using MD simulations (see Methods for details). The elastic, plastic, and fracture mechanical properties of MAPbI3 were characterized and found to be very soft and nanoductile, with major contributions from interatomic Coulombic interactions (determined through atomic stress decomposition, see eq 8 in Methods), an important feature for ionic crystals. We paid special attention to the underlying mechanism of the nanoductility observed during the microstructural evolution. A local plastic deformation region was found at the nanoscale, while only localized partial edge dislocations were observed near the crack tip. The movement of these dislocations forms nanovoids, which enlarge and coalesce with adjacent nanovoids to achieve crack propagations. Finally, we compared the mechanical properties, in particular, stiffness, fracture toughness, and flexibility, of hybrid perovskites against existing energy-related inorganic PV, LED, and thermoelectric materials.

rather than the ultimate stress−strain predicted for pristine, crack-free materials. There has been rapid progress in studies of the electronic band structures and exciton diffusion pathways of hybrid perovskites,9,10 but very few experimental or theoretical investigations have been reported on their mechanical properties. Similar to silicon-based11,12 solar cells and inorganic LED devices (such as those using InGaN, ZnSe, GaP, and etc.), there are great concerns that applications of hybrid perovskites in energy-related devices may be limited by their brittleness and lack of mechanical flexibility compared to organic PV and OLED (organic LED) materials based on conductive polymers or small molecules. In particular, the stress state (residue or thermal stresses, cyclic bending, and surface scratching) and the crystallinity of the hybrid perovskite layer have strong impacts on the solar spectrum absorption performance, optical bandgap, and the resulting energy conversion efficiency.13 It also has large impacts on the structural phase transition,14 material degradation and mismatch at interfaces, and the resulting device lifetime, as studied in recent experiments on making flexible and fatigue-resistant hybrid perovskite solar cells.15,16 Experimental methods, while certainly of great value, have limitations in determining the long-term stability of hybrid perovskites, measuring their fracture mechanical properties, and decoupling the various atomistic interaction contributions to the above properties. Indeed, direct measurements of the elastic, plastic, and fracture responses of hybrid perovskites through tensile testing are challenging, since high sample quality (with large crystal sizes, or ideally, be single crystalline rather than polycrystalline) and large samples are required. For example, so far only nano-indentation tests have been carried out very recently to measure the hardness and deduce the Young’s modulus (10−20 GPa) of MAPbX3 (X = Br, I, Cl) at the room temperature.17,18 On the other hand, theoretical and computational methods, from first-principle calculations to force field-based molecular dynamics (MD) simulations, can predict material properties that are difficult to access experimentally. Recently, the elastic constants and elastic mechanical properties of pristine hybrid perovskites MABX3 (B = Sn, Pb; X = Br, I) have been determined using the density functional theory (DFT) for their cubic, tetragonal, and orthorhombic phases.19 Despite the above DFT19 and experimental17,18 studies on the elastic properties of hybrid perovskites, fundamental understanding and quantification of 11045

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

Figure 1. (a−c) Crystal structures of MAPbI3 under the orthorhombic, tetragonal, and cubic phases and the corresponding 3D illustrations of the surface contours of the MM-predicted Young’s modulus E. The x, y, and z axes used here reflect the [100], [010], and [001] crystallographic directions. (d) Stress−strain curve of the pristine cubic single and polycrystalline MAPbI3 (consisting of 6, 10, and 14 grains) under tensile loading. σ*zz is the effective tensile stress along the z axis considering the Poisson’s effect, and ε is the engineering strain along the z axis. The simulation setup and the chemical structure of MAPbI3 are shown in the top right corner.

RESULTS AND DISCUSSION We first validated the capability of the current MYP force field (see Methods for details) to reproduce the elastic properties of hybrid perovskites. The MYP potential has been shown to capture well the structural (temperature-dependent lattice constants for different phases), dynamical (rotations of MA+ cations and distortions of the PbI6 octahedra), and some mechanical properties (bulk modulus and cohesive energy) of MAPbI3.20,21 Force field-based MM simulations are conducted using the volume perturbation method in GULP22,23 to calculate the first-order elastic constant matrix [cij] and compliance constant matrix [sij] of the orthorhombic, tetragonal, and cubic phases of MAPbI3. The calculated cij values are summarized in Table S1 in the Supporting Information, and they are consistent with the DFT results,19 although small deviations were observed for c44. Overall, MM calculations predict slightly higher elastic constants than DFT. Based on the Voigt−Reuss−Hill (VRH) approximation,24

elastic constants cij can be utilized to evaluate various mechanical properties of MAPbI3, including Young’s modulus E, bulk modulus K, shear modulus G, and Poisson’s ratio ν (see Table 1). The VRH approach combines the upper (Voigt) and the lower (Reuss) bounds of the mechanical properties as the macroscopic effective mechanical properties. The mechanical anisotropy is another important factor in the fracture mechanics of MAPbI3, since the formation and propagation of microcracks are often related to the elastic anisotropy. To better describe the anisotropic elastic behavior of each phase, we constructed the three-dimensional (3D) surface contour plots for E with respect to different crystallographic directions, as depicted in Figure 1a−c (see Methods for details). In principle, for an isotropic crystal, the 3D surface contour should be spherical. As shown in Figure 1a,b, the orthorhombic and tetragonal phases show strong elastic anisotropy. The E value in the [010] direction for the orthorhombic phase and [001] direction for the tetragonal 11046

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

Figure 2. Microstructural evolution of pristine cubic-phase single and polycrystalline MAPbI3 (6-, 10-, and 14-grains) during the tensile loading process. Plots (a−c) are simulation snapshots colored based on the centro-symmetry parameter (CSP, defined by eq 11 in Methods). For single crystals, the corresponding strain values are 5.0%, 15.4%, and 35.0% for (a1), (b1) and (c1), respectively. For polycrystals, the corresponding strain values are 4.7%, 32.2%, and 58.2% for (a2−a4), (b2−b4), and (c2−c4), respectively.

in experimental nano-indentations and various theoretical methods (DFT and our MM and MD simulations) suggest that hybrid perovskites could be ductile for extensive deformation, primarily attributed to the low shear modulus G. Compared to MM results, we found that the temperature effect considered in MD simulations has a great influence on the predicted elastic properties and improves the predicted elastic properties compared to experiments.17,18 With the above validations, we are confident about the capability of the MYP force field utilized here to capture the mechanical properties of MAPbI3. Since MAPbI3 is deployed as a thin film composed of small grains in hybrid perovskite solar cells, here we also evaluate the mechanical properties of pristine polycrystalline MAPbI3. We studied three polycrystalline MAPbI3 models possessing different averaged grain sizes of 8.25−12.60 nm (or the inverse of the number of grains from 14 to 6, see Figure S1 in the Supporting Information). Based on the stress−strain curves in Figure 1d, their Young’s moduli E and Poisson’s ratio ν are in the ranges of 5.08−7.89 GPa and 0.16−0.23, respectively. These values are lower than those of single-crystalline MAPbI3, demonstrating the mechanical deterioration induced by grain boundaries. Remarkably, there exists obvious yielding regions in the stress−strain curves. The corresponding yield stresses and yield strains for all three polycrystalline MAPbI3 are in the ranges of 0.133−0.181 GPa and 3.4−4.6%, respectively, which are much lower than the ultimate stress and strain of singlecrystalline MAPbI3. It is also unveiled that the yield stress is reduced when decreasing the averaged grain size. Such a trend in hybrid perovskite (see Figure S1) is in good agreement with the recently observed inverse Hall−Petch (grain boundary strengthening) effect in polycrystalline grapheme.26 Reductions of the yield stress here can be attributed to the amorphous structures at the grain boundaries. The amorphization process normally weakens the fracture strength of solids,27 and therefore, polycrystalline MAPbI3 with smaller grains, which

phase have the lowest values. This is consistent with the fact that both directions have the largest stacking spaces and thus lowest packing densities. While for the cubic phase (Figure 1c), higher symmetry contributes to its weaker elastic anisotropy: the values of E along x, y, and z axes are all ∼25.5 GPa. The MM-predicted Young’s modulus E for MAPbI3 is in the range of 14−35 GPa for the different phases. These values are higher than those measured by recent nano-indentation tests17,18 (see Table 1) due to the neglect of temperature. MAPbI3 exhibits three types of temperature-dependent crystal structures, that is, the orthorhombic, tetragonal and cubic phases. For many MAPbI3-based functional devices (such as solar cells and LED), their working temperatures are higher than the tetragonal-to-cubic phase transition temperature (327 K).25 It is well-known that MM calculations usually overpredict the elastic constants due to the neglect of temperature. Therefore, we carried out mechanical tensile test on pristine cubic single and polycrystalline MAPbI3 using MD simulations at 350 K (see Figure 1d) to directly determine the various elastic properties of MAPbI3, as listed in Table 1. The effective tensile stress σ*zz (considering the Poisson’s effect for the fixed simulation box dimensions, see eq 9 in Methods) is computed during the tensile loading process and shown in Figure 1d as a function of the engineering strain ε. For pristine singlecrystalline MAPbI3, the resulting ultimate tensile stress and yield strain are 0.449 GPa and 5.0%, respectively. From the linear elastic region of the stress−strain curve and separate MD simulations with varying normal pressure tensors under the NPT ensemble, we found that the Young’s modulus E = 14.94 GPa and the Poisson’s ratio ν = 0.31, respectively, which were then used to estimate the bulk and shear modulus K and G in Table 1 for the cubic phase (MD). The K/G ratio has been previously identified as a measure of ductility of hybrid perovskites,19 while lower ν values were connected with larger compressibility or deformability of solid materials. The high K/ G ratios of 1.14−5.03 and low ν values of 0.16−0.36 obtained 11047

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

Figure 3. Local von Mises stress distribution of pristine single and polycrystalline MAPbI3 (with 6-, 10-, and 14-grains) during tensile loading. Plots (a−c) are colored based on the values of local von Mises stresses. For single crystals, the corresponding strain values are 5.0%, 15.4%, and 35.0% for (a1), (b1), and (c1), respectively. For polycrystals, the corresponding strain values are 4.7%, 32.2%, and 58.2% for (a2−a4), (b2− b4), and (c2−c4), respectively.

Figure 4. (a) Stress−strain curve during the tensile fracture and crack propagation of defective single-crystalline MAPbI3. The defective model of MAPbI3 (cubic phase structure) with a crack is shown in the plot. The stress−strain curve is colored to denote the various deformation and fracture stages during the tensile test. (b) Decomposed contributions from various interatomic interactions to the total stress−strain curve. The total effective stress was decomposed into contributions from short- and long-range Coulombic, vdW, bond, angle, and dihedral interactions. The nonbonded pairwise contribution to the stress is the sum of the vdW and Coulombic (short-range) contributions.

transformations in the yielding region. The crystalline structure is gradually absorbed by the amorphous grain boundaries and becomes highly amorphous without forming any nanovoids. Similar phase transformations were also confirmed to be a dominant yielding mechanism in polycrystalline aragonites (carbonate minerals).28 Here the large fraction of amorphous structures in polycrystalline MAPbI3 is believed to greatly facilitate such an amorphization process and hence results in the large yielding region. Interestingly, the stress enhancement in the yielding region for the 10- and 14-grain cases (compared to the 6-grain case) also reveals the strain hardening behavior of polycrystalline MAPbI3 with smaller grain sizes. This may be due to the grain rotation and grain-boundary sliding29,30 phenomena observed here (see Movies S1, S2, and S3 in the Supporting Information). The von Mises stress distribution is relatively uniform in single-crystalline MAPbI3, as shown in Figure 3a1. Local stress

correspondingly contain larger fractions of amorphous structures at the grain boundaries, will exhibit lower yield stress. We now discuss the surprising nanoductility of polycrystalline MAPbI3. In Figure 1d, the yielding regions of the three polycrystalline MAPbI3 are more extensive than that of the single-crystalline MAPbI3. To explore the underlying yielding mechanism, the simulated deformation configurations of both single and polycrystalline MAPbI3 at different tensile strains are shown in Figure 2. For single-crystalline MAPbI3, it deforms homogeneously during tensile loading (Figure 2a1). When reaching the ultimate tensile stress (∼0.449 GPa), a phase transformation (amorphization) from crystalline to amorphous begins, but is mitigated and suppressed by the formation of nanovoids (Figure 2b1). The adjacent nanovoids gradually merge (Figure 2c1), which is also accompanied by significant stress drops. For polycrystalline MAPbI3, Figure 2a2−c2, a3−c3, and a4−c4 shows that there are extensive and continuous phase 11048

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

Figure 5. Crack propagation at different stages. Plots a−c are colored based on the local von Mises stress; plots d−f are colored based on the CSP. The strains are 4.4%, 11.9%, and 44.7% in (a−c), respectively.

concentrations lead to lots of nanovoids (see Figures 2b1 and 3b1, which then emerge into a crack, releasing some stress (see Figure 3c1). There are some stress concentrations at the grain boundaries, as shown in Figure 3a2−a4, for polycrystalline MAPbI3. Smaller grain-size models possess more grain boundaries and show lots of stress concentrations at the grain boundaries, which lead to lower yield stresses and the inverse Hall-Petch effect. As the transformation from crystalline to amorphous structure proceeds, the stress concentrations will move from the grain boundaries to other regions, as shown in Figure 3b2−b4. After yielding, the system becomes unstressed, and the local von Mises stresses reduce to small values, as shown in Figure 3c2−c4. The single crystalline MAPbI3 is then selected in MD simulations to investigate its microstructural evolution during fracture and the resulting fracture toughness. We constructed a defective model for the single-crystalline MAPbI3 with a rectangular crack at the center and applied the tensile test (see Figure 4a). The central crack is terminated with the more stable charge-neutral (001) surfaces (i.e., [MAI]0 or [PbI2]0). The * is computed during the tensile loading effective tensile stress σzz process and shown in Figure 4a as a function of the engineering strain ε. From the linear elastic region of the stress−strain curve, we estimated that the Young’s modulus E is ∼14.04 GPa (see Table 1, lower than that of the pristine case due to the presence of crack), in perfect agreement with the experimental nano-indentation results in the range of 10.4−14.3 GPa.17,18 We also estimated the corresponding bulk and shear modulus, K and G, for the defective model in Table 1, in very good agreement with the experimental results. The predicted ultimate tensile stress and yield strain are 0.319 GPa and 4.4%, respectively, determined from the linear elastic region. The sharp drop of σ*zz after the linear elastic region reflects the brittle fracture nature, but it is then followed by a very slow decay of stress (stringing effect) until complete failure (σ*zz

approaches zero) at a large strain of 44.7%. Interestingly, MAPbI3 does not behave like traditional ceramics or inorganic oxide-based perovskites (such as SrTiO331) that are classical brittle materials. A significant plastic deformation region is observed with a constant stress level (half of the ultimate tensile stress), which occurs at a strain range of 8.1−24.4%. These observations imply strong electrostatic attractions between MA+ cations and octahedral PbI6 cages (discussed later). The colored mapping of the local von Mises stress (defined by eq 10 in Methods) of MAPbI3 during the tensile test is shown in Figure 5a−c. The crack does not propagate under a strain of 4.4% as shown in Figure 5a, but some stress concentration is found near the crack tip with local disordering of the crystal structure. The crack first gets blunt when the strain is 8.1%, then a number of nanovoids (small cavitations) are formed during the initial stage of plastic deformation in Figure 5b. However, the crack blunting and generation of nanovoids generally reflect the plastic deformation feature which is rarely observed in brittle ceramic materials.32 The number of nanovoids gradually increases when increasing the strain, which was also observed in defective amorphous/glassy ceramic and polymeric materials under tensile loading.33,34 When the strain reaches 12.7%, the nanovoids in the vicinity of the crack tip coalesce quickly during the last stage of plastic deformation. The crack gradually propagates when the strain is increased up to 24.4%. Complete fracture finally occurs at the strain of 44.7% in Figure 5c, and the system becomes unstressed and the local von Mises stress reduces to a small value, which also indicates that the strain energy is gradually released in the process of crack propagation. To decouple the effects from various interatomic interactions to the total stress, Figure 4b illustrates contributions, including van der Waals (vdW, modeled by the Buckingham and Lennard-Jones potentials here), short- and long-range (via Ewald summation) Coulombic, and bonded (i.e., covalent 11049

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano bonding, angular bending, and dihedral rotation within MA+) interactions, to the total effective stress. We found that the contributions from both long-range Coulombic and bonded interactions have negligible impact on the effective stress. Meanwhile, the contribution from nonbonded pairwise interactions (sum of the vdW and short-range Coulombic interactions) plays a leading role. The contributions from the short-range Coulombic interactions to the effective stress are positive, as shown in the inset in Figure 4b, suggesting that net electrostatic attractions hold the ionic components together and fight against the external tensile stress. On the other hand, vdW interactions generate net repulsive (negative) stress and follow the direction of the external tensile stress. The magnitude of the short-range Coulombic contribution is larger than that of the vdW contribution, suggesting that the former is the dominant factor in determining the mechanical properties of MAPbI3. By probing the microstructural evolution at the atomistic level, we captured and examined the various brittle and ductile responses, including the competing processes of brittle fracture (nanovoids formation) and amorphization35 (PbI6 octahedra disordering) in the single-crystalline MAPbI3. The colored CSP mapping illustrates the degrees of local defects and partial dislocations during the tensile test, confirming the nanoductility (featured by nanovoid formation and evolution) of the singlecrystalline MAPbI3 within a glassy region33 around the crack surface, as shown in Figure 5d−f. Large structural distortion and stress concentration exist near the crack tip with high CSP values in Figure 5d, indicating that its structure is rearranged and suggesting that there might be some dislocations near the crack tip in response to strain. The region with high CSP values becomes larger when increasing the strain in Figure 5e, accompanied by the enhanced disordering of the singlecrystalline MAPbI3, and finally, forms a plastic deformation zone (glassy state of MAPbI3). Interestingly, the von Mises stress distribution remains almost unchanged for a large range of strain values without crack propagation. In Figure 5f, the path of crack propagation is located within the region of high CSP values (see Movie S4 in the Supporting Information), indicating that the crack propagates through local defects. The existence of localized partial edge dislocations at the crack tip within MAPbI3 under small tensile strains is shown in Figure S2. This suggests that MAPbI3 has a certain degree of shear-induced ductility at the nanoscale. It is known that the polar character of perovskite oxides and their associated charge effects remarkably affect the formation of dislocations in perovskite oxides, such as SrTiO3.36 Similar charge effect in MAPbI3 may also play a notable role in their formation of dislocations, although dislocations in MAPbI3 are significantly more difficult to capture than those in metals. Moreover, atoms in MAPbI3 are quickly rearranged under a small strain of 4.4% and easily transit from ordered to disordered states via the slipping of dislocations. The formation of nanovoids near those disordered atoms leads to the observed brittle-fracture feature in MAPbI3. Although dislocation is a common phenomenon observed during plastic deformation, we found very few dislocations in MAPbI3, while the stress−strain curve suggests a wide plastic deformation region. This indicates that the dislocation observed here is not the primary reason for the plastic deformation in MAPbI3. The underlying reason lies in the existence of strong Coulombic interactions and nanovoid expansion during plastic deformation. In order to testify this

hypothesis, some simulation snapshots (Figure S3) are shown at various stages of plastic deformation. We found that the ductile behavior of MAPbI3 is mainly attributed to the enlarging of individual nanovoids and the coalescence of adjacent nanovoids. Figure S3 shows the process of nanovoid expansion and coalescence at different strain values. Ionic components in MAPbI3 are rearranged to form the nanovoids after the initial brittle fracture process. The two red ovals in Figure S3a highlight two representative nanovoids adjacent to each other at a strain of 18.4%. The two adjacent nanovoids are gradually enlarged (growing) near the crack tip when the strain is increased. Then, only ionic chains are present to connect the two sides of the material (necking) and separate the two enlarged nanovoids when the strain reaches 22.8% in Figure S3b. The adjacent nanovoids finally merge under higher strain values, with the ionic chains completely broken (see Movie S5 in the Supporting Information). The complete breakage of the ionic chains near the crack tip will initiate crack propagation at a strain of 24.4% in Figure S3c. From a lattice structure point of view, the plastic deformation in MAPbI3 is associated with the distortion and splitting of the building-block PbI6 octahedra, as shown in Figure S3d,e. The splitting of PbI6 octahedra to PbI4 tetrahedra contributes to the irreversible crack propagation during tensile loading in Figure S3f,g (for more details, see Figure S4). It has been reported that the initial failure in a perfect crystal is induced by atomistic disorders,37,38 which promotes the formation of nanovoids, similar to those observed in Figure S3a. The populating, growing, and coalescing of nanovoids with similar sizes observed here finally lead to ionic chain breakage in MAPbI3, similar to the fracture behavior observed in MD simulations for skutterudite CoSb3.39 To quantify nanovoid migration or population, the lengths of the simulated region populated with nanovoids as a function of strain are shown in Figure S5. These results indicate that the length of the region populated with nanovoids increases linearly with strain during the brittle fracture process, which is also directly correlated with the linear decrease in the effective stress in Figure 4a (see Movie S6 in the Supporting Information). Based on the method developed by Brochard et al.40 and Bauchy et al.,33,34 we determined the fracture mechanical properties of MAPbI3 at the nanoscale, without any specific assumptions about the mechanical behavior of the material during fracture. Common methods, such as the numerical and analytical J-integral41,42 and fictitious crack model of Hillerborg,43 have been used to determine the fracture energy release rate for either brittle or ductile materials under the assumption that the material is homogeneous and continuous in the macroscale.34 Therefore, these methods are not applicable for studying MAPbI3 because it is mechanically anisotropic and possesses unique fracture properties at the nanoscale. The method we used is based on the energy theory of fracture mechanics and uses thermodynamic integration (over the stress−strain curve during crack propagation) to determine the fracture toughness.44 This method has been successfully used to predict the fracture toughness of materials possessing significant plasticity due to nanovoids, such as silicate glasses33 and amorphous calcium silicate hydrate.34 During tensile loading, the internal stresses of a system increase from zero, together with the associated mechanical (strain) energy P. When the crack starts to propagate, the mechanical energy stored is released to create two new surfaces. Based on the energy theory of fracture mechanics, the Helmholtz free energy F can take the form of the mechanical 11050

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

Figure 6. (a) Total potential energy (Epot) profile of the simulated MAPbI3 system as a function of strain. The stress−strain curve in Figure 4a is also shown here for comparison. (b−d) Decomposed Epot into various contributions from interatomic interactions: (b) nonboned pairwise (Coulombic + vdW) interactions, in which both short- and long-ranged Coulombic energies are included, (c) separated Coulombic and vdW interactions, and (d) bonded molecular interactions.

energy P under the NVT ensemble. Therefore, the energy released per unit area (the xy plane in Figure 4) of the propagated crack defines the critical strain energy release rate Gc, an important measure of fracture toughness based on the Griffith’s theory of fracture:40,45 Gc =

LxLy ΔA∞

∫L

Lzmax

0 z

*dLz σzz

MAPbI3 is quite low compared to traditional inorganic perovskites such as LaCoO347 (0.73−0.98 MPa·m1/2) and BaTiO348,49 (0.68−1.0 MPa·m1/2). Combined with its low elastic modulus, MAPbI3 is both soft and easier to fracture (from a stress point of view) but exhibits plasticity and ductility (flexible from a strain point of view) at the nanoscale, which is absent in most inorganic perovskites. To quantify the degree of ductility of MAPbI3, we decomposed Gc into33

(1)

where ΔA∞ is the total area of the crack created after complete fracture (under a strain of 44.7% here), σzz * is the effective stress along the z direction, Lx (fixed), Ly (fixed), and Lz (varies from L0z to Lmax as the strain is increased from 0 to 44.7%) are the z sizes of the simulation box in the x, y, and z direction, respectively. Another important measure of fracture toughness is KIc, the critical stress intensity factor for the Mode I fracture that is commonly used in engineering applications. KIc denotes the ratio between the critical stress at the crack tip right before crack propagation and the external loading stress far away from the crack tip. KIc can be estimated from Gc based on the Irwin formula for the plane strain loading condition considered here and under the assumption of mechanical isotropy:46 KIc =

⎛ GcE ⎞1/2 ⎜ ⎟ ⎝ 1 − ν2 ⎠

Gc = Ge + Gd

(3)

where Ge is the reversible, elastic contribution to the fracture energy release rate in the linear elastic stage, and Gd is the irreversible, dissipated energy release rate during plastic deformation and crack propagation. The value of Gd for a perfectly brittle material should be equal to zero. We used the brittleness parameter, B = Ge/Gc, to quantify ductility, which is larger for more brittle materials. The computed value of B for MAPbI3 is 0.20, suggesting that it is even more ductile than amorphous glasses, such as sodium silicate (B = 0.73) and calcium aluminosilicate (B = 0.38).33 To further understand the unexpected ductility of MAPbI3, we evaluated the total surface energy of typical charge-neutral [MAI]0 and [PbI2]0 free surfaces along the (001) plane, which gives γ = γMAI + γPbI2 = 0.19 J/m2 (see Methods). The total surface energy, γ, should be very close to Ge = 0.188 J/m2 (determined from the linear elastic deformation stage as in eq 3) in order to initiate crack propagation, which is indeed the case here, which again validates our predicted fracture properties of MAPbI3. However, such brittle fracture through the creation of free surfaces does not proceed for long (for 4.4− 8.1% strains) before the plastic deformation begins. Here the much smaller value of Ge compared to that of Gc (i.e., very small B value) indicates more energy dissipation during the

(2)

where E = 14.04 GPa is the Young’s modulus predicted by MD for the defective MAPbI3 model, and ν = 0.31 is the Poisson’s ratio for the cubic-phase MAPbI3 predicted by MD and given in Table 1. The critical strain energy release rate Gc is calculated based on eq 1 by integrating the entire stress−strain curve in Figure 4a, which leads to Gc = 0.94 J/m2, and subsequently based on eq 2, KIc = 0.12 MPa·m1/2. The fracture toughness observed in 11051

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

Table 2. Comparison between the Mechanical Properties (Young’s Modulus E, Critical Stress Intensity Factor KIc, Critical Strain Energy Release Rate Gc, and Poisson’s Ratio ν) of MAPbI3 Predicted in This Work and Those of Traditional Inorganic PV (e.g., Silicon and GaAs), Thermoelectric (e.g., Bi2Te3 and PbTe), and LED (e.g., ZnSe, GaP, GaN, and InGaN) Materials Available from Various Experimental Studiesa

a

materials

E (GPa)

Si GaAs Bi2Te3 PbTe ZnSe GaP GaN InGaN pristine MAPbI3 defective MAPbI3 polycrystal MAPbI3

53

KIc (MPa·m1/2) 54

174.8, 150 121.5,57 85.5,58 11859 39.9,63 31−3464 58.0566 78.667 144.657 32072 ∼37675 14.94b 14.04b 5.08−7.89b

55

56

0.90, 0.67−0.95 0.46,60 0.44,58,61 0.5562 1.12,64 1.13−1.1564 0.3566 0.968 1.270 0.7973 − 0.15b 0.12b 0.10−0.13b

Gc (J/m2)

ν

1.253 0.8757 − − 0.667 0.9457 − − 1.45b 0.94b 1.93−2.04b

0.2754 0.3158 0.2565 0.2666 0.2869 0.3171 0.2274 0.3076 0.31b 0.31b 0.16−0.23b

Could be either single or polycrystalline depending on the synthesis method used. bThis work.

CONCLUSIONS Using atomistic simulations and solid mechanics theories, we investigated the elastic, plastic, and fracture mechanics of pristine and defective MAPbI3. Force field-based MM simulations were first conduced to evaluate the elastic constants of pristine orthorhombic, tetragonal, and cubic single-crystalline MAPbI3. Their elastic mechanical properties (i.e., Young’s modulus, bulk modulus, shear modulus, and Poisson’s ratio) based on the Voigt−Reuss−Hill approximation were predicted as well. According to the 3D surface contours of Young’s moduli of MAPbI3, its orthorhombic and tetragonal phases exhibit strong elastic anisotropy. While for its cubic phase, higher symmetry contributes to its weaker elastic anisotropy. The results calculated by the DFT, MM (in this work), and MD (in this work) simulations on pristine MAPbI3 are in good agreement with the experimental nano-indentation results. With the above validations, we are confident about the capability of the force field utilized here to capture the mechanical properties of MAPbI3. The mechanical properties of polycrystalline MAPbI3 possessing different grain sizes were evaluated as well. MD simulation results revealed that polycrystalline MAPbI3 is softer (lower Young’s modulus) and possesses much lower ultimate strength than the singlecrystalline one. Meanwhile, they exhibit a larger yielding region during tensile loading, and such notable nanoductility can be attributed to the extensive and continuous amorphization. An inverse Hall−Petch relation, in which the yield stress rises with the increase of gain size, was also observed in polycrystalline MAPbI3. Then we investigated the fracture toughness and microstructural evolution of the defective cubic MAPbI3. A significant plastic deformation (amorphous/glassy) region was observed after the initial elastic deformation and the limited brittle fracture process. Even after breakage, the stress value remains nonzero up to a large strain of 44.7%, suggesting strong Coulombic attractions between the separated ionic components. The predicted Young’s modulus and fracture toughness suggest that hybrid perovskites are relatively soft and easy to fracture, mechanically similar to reinforced polymers or glassy ceramics. Only localized edge dislocations were found in MAPbI3 near the crack tip, and the observed limited brittle and significant plastic behavior originate from the nanoductility nature of ionic MAPbI3 crystals. Specifically, the slipping of dislocations leads to nanovoid formation, and the limited

plastic deformation and fracture stage than that during the elastic stage, which leads to the ultimate large ductility. Notably, the underlying mechanism behind such trend (Ge ≪ Gc) in MAPbI3 is different from that for some brittle materials, such as Si50 and SiC,51 where the lattice trapping effect dissipates a portion of the strain energy and leads to the inequality between Ge and Gc. Finally, we analyzed the total potential energy (as an estimate of the Helmholtz free energy) as a function of strain to understand the microstructural evolution of MAPbI3 during the fracture process. The contributions from distinct interactions to the total potential energy were also considered in Figure 6b−d, with trends consistent with those observed from contributions to total stress in Figure 4b. It was noteworthy that the kinetic energy was kept constant here via temperature control. The total potential energy increases (or decreases in magnitude) when increasing strain in the linear elastic and the first half of the brittle fracture stage in Figure 6a. Due to the existence of partial dislocations, the potential energy continued to rise in the first half of the brittle fracture stage. Then the potential energy dropped slightly (stress released) in the second half of the brittle fracture stage, as a result of nanovoid formation. The primary contributor to the reduction of the total potential energy during nanovoid spreading (brittle fracture region) is the increase of the Coulombic potential energy (in magnitude), as shown in Figure 6c. As the distances between atoms are increased due to nanovoid formation, Coulombic interactions are decreased (with the energy increased in magnitude), which leads to stress release. When individual nanovoids start to enlarge, the total potential energy increases again, reaching the maximum value at the strain of 24.4%, primarily due to the increase in the vdW potential energy. Once the crack begins to propagate through nanovoid coalescence and ionic chain breakage, a large amount of potential energy is released, further reducing the Coulombic and bonded potential energies. When crack propagation is completed, the total potential energy drops to a lower value, but is still higher than that before tensile loading due to the formation of two new surfaces. Moreover, the total potential energy does not fall sharply during the crack propagation stage, suggesting a low crack propagation speed. This is similar to the fracture behavior of amorphous calciumsilicate-hydrate, again confirming the soft and ductile-fracture nature of MAPbI3. 11052

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

charge-neutral surfaces [MAI]0 and [PbI2]0 on the top and bottom of the crack. It is well-known that the reliability of classical force fields is crucial for reproducing the mechanical and fracture properties of materials.78,79 To predict the mechanical properties of single and polycrystalline MAPbI3, we used a recently developed classical force field for hybrid perovskites (MYP) to simulate MAPbI3.20 It has been confirmed that the MYP force field can nicely reproduce the temperature-dependent structural, dynamical, and physical properties of MAPbI3, including: (i) the phase diagram and rotation and angular distribution of MA+ cations,20 (ii) point-defect diffusivity that could be involved in extended defect propagation,80 (iii) phonon thermal transport that is important for energy release rate,81,82 and (iv) electroand mechanocaloric effects related to the interplay of molecular ordering and entropy.83 Since these properties are also highly related to the elastic and plastic responses, energy release rates, and crack/ defect propagation in MAPbI3, we expect that the MYP force field is also capable of predicting the mechanical and fracture properties of MAPbI3, as verified in Table 1 and Figure 1 after comparing to firstprinciple DFT results. However, we need to recognize the limitations of the present MYP force field that some phenomena are still beyond its capability and remain to be further validated. For example, degradation and decomposition mechanisms of MAPbI3 involving chemical reactions are still beyond the capability of the present MYP force field, such as the deprotonation of MA + into metylamine and molecular fragmentation.84 Even the most likely decomposition process85 (i.e., molecule fragmentation, which remains to be tested in our future work) of MAPbI3 → PbI2 + MAI might be beyond the capability of the present MYP force field, since it was originally not designed to model individual PbI2 and MAI crystals. In the present MYP force field, the cross-term interactions between inorganic and organic components are described by a combination of Buckingham, Lennard-Jones (LJ), and electrostatic (Coulombic) terms. In summary, the total nonbonded potential energy of MAPbI3 can be calculated by the expression:20

brittleness is due to the increasing number of nanovoids and their migrations. The significant plastic behavior is mainly attributed to the expansion of individual nanovoids and the coalescence of adjacent ones. The adjacent nanovoids gradually merge (necking), with the two sides of the crystalline material (separated by the glassy region) bridged or held only by ionic chains. Finally, the crack starts to propagate at a strain of 24.4% when some ionic chains start to break near the crack tip. Our theoretical predictions and understanding of fracture mechanics of MAPbI3 will enable the rational design of hybrid perovskite-based functional devices for energy-related applications with high stability and mechanical flexibility. One possible strategy could involve the design of highly charged (or higher valence) or polarized ionic species that hold each together more tightly via electrostatic interactions. In addition, our predicted fracture energy release rates in solid state (correlated to the vacuum permittivity) establish an upper bound for degradations of hybrid perovskites in liquid state or as composites (correlated to the dielectric constant of the liquid or matrix-phase), in which the electrostatic attractions among ionic species are reduced due to dielectric constants. Finally, in Table 2, we compared our predicted fracture mechanical properties of MAPbI3 with those of existing inorganic PV (e.g., silicon and GaAs), thermoelectric (e.g., Bi2Te3 and PbTe), and LED (e.g., ZnSe, GaP, GaN, and InGaN) materials reported from experimental studies. As a mixed organic−inorganic material, hybrid perovskites feature electrostatic or weak ionic bonding that dominates their mechanical properties and grain boundary structures, which is difficult to predict based on the simple combination of bonding features of soft organic (vdW interactions dominating) and hard inorganic (covalent or strong ionic interactions dominating) components. We found that MAPbI3 possesses the lowest Young’s modulus (being softest) and the critical stress intensity factor among all these materials due to its ionic bonding nature, suggesting that care must be taken in device fabrication with respect to its poor wear resistance and vulnerability against high local stress or impacts. On the other hand, we found that common polycrystalline MAPbI3 possesses the highest fracture energy and a very low Poisson’s ratio due to its remarkable ductility against extensive deformation, suggesting that it is well suited as a thin film in flexible solar cells,15,16 LEDs,52 and thermoelectric devices that can sustain many cycles of bending and distortion.



Uij =

∑ ⎢Aije−rij/ ρij − i,j

⎢⎣

Cij ⎤ ⎥+ rij6 ⎥⎦

∑ i≠j

1 qiqj 4πε0 rij

⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij + ∑ εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥(rij ≤ rc) r ⎝ rij ⎠ ⎦ i,j ⎣⎝ ij ⎠

(4)

where rij is the interatomic distance between atoms i and j; Aij, ρij, and Cij are the Buckingham parameters; qi and qj are the partial charges of atoms i and j; εij and σij are the LJ parameters; and rc is the cutoff distance for nonbonded interactions. Specifically, the pairwise interactions between Pb or I and heavy atoms (C and N) in MA+ are described by the Buckingham and Coulombic terms to stabilize the perovskite crystalline structure. The pairwise interactions between Pb or I and H atoms in MA+, as well as the pairwise interactions within MA+, are described by the LJ and Coulombic terms. The cutoff distances for both Buckingham/LJ and short-range Coulombic interactions are selected as rc = 10 Å. Molecular interactions of organic MA+ are described by the AMBER force field. The Ewald summation technique is also adopted to deal with long-range Coulombic interactions.86 The directional dependences of the reciprocal values of the Young’s modulus E (shown in Figure 1a−c) for the various phases of MAPbI3 can be defined as87 Orthorhombic:

METHODS Pristine cubic-phase single and polycrystalline MAPbI3 models with dimensions of 1.887 nm (x-axis) × 49.830 nm (y-axis) × 19.122 nm (z-axis) (as shown in Figure 1d) were built separately to evaluate their corresponding elastic and fracture properties, with periodic boundary conditions applied in all three directions. Three polycrystalline MAPbI3 models, containing 6-, 10-, and 14-grains (or possessing averaged grain sizes of 8.25, 9.76, and 12.60 nm, respectively), were generated using the Voronoi construction method.77 To avoid artificial atom overlaps in the polycrystalline MAPbI3 models, atoms protruding beyond the grain boundaries were deleted by setting the gap distance to be 0.8 Å. Then we prepared the defective single-crystalline defective MAPbI3 models to investigate the fracture performance of MAPbI3. A central crack was then introduced by removing certain central ions from the pristine cubic single-crystalline model. In order to obtain a strong stress concentration at the crack tip, we built the rectangular-shape crack with the dimensions of 10.164 nm (y-axis) × 1.747 nm (z-axis), as shown in Figure 4. The atoms on the crack surface were also carefully deleted to keep the model charge neutral, by rendering two

1 = l14s11 + 2l12l 22s12 + 2l12l32s13 + l 24s22 + 2l 22l32s23 + l34s33 E + l 22l32s44 + l12l32s55 + l12l 22s66

(5)

Tetragonal: 11053

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano 1 = (l14 + l 24)s11 + l34s33 + l12l 22(2s12 + s66) + l32(1 − l32) E × (2s13 + s44) + [2l1l 2(l12 − l 22)s16]

N /2

CSP =

where the N nearest neighbors or each atom are identified, and R⃗ i and R⃗ i+N/2 are vectors from the central atom to a particular pair of nearest neighbors. Here CSP is only evaluated for I atoms with N = 12 and used to characterize the local disorder of octahedra under tensile loading. Bulk and surface-exposed MAPbI3 models (cubic-phase MAPbI3 with dimensions of 3.758 nm (x-axis) × 3.701 nm (y-axis) × 3.833 nm (z-axis)) were generated separately to calculate the total surface energy of the typical charge-neutral [MAI]0 and [PbI2]0 free surfaces along the (001) plane (see Figure S6). We set the same vacuum boxes with a length of 3.928 nm along the z-axis on both the top and bottom of the surface-exposed model in order to avoid any short-range vdW and electrostatic interactions and minimize long-range electrostatic interactions between the top and bottom surfaces due to the periodic boundary condition. The total surface energy, γ = γMAI + γPbI2, can be defined as the following:

(7)

where l1, l2, and l3 are the directional cosines in the transformed Cartesian coordinate system with respect to the original coordinates, and [sij] is the compliance matrix such that [cij][sij] = 1. In order to obtain the initially equilibrium state, the pristine and defective models were first relaxed for 100 ps in the NVT ensemble after energy minimization and then run for another 100 ps in the NPT ensemble88 to reduce all the possible residual stresses inside the model. During the tensile test, the simulation model was stretched along z-axis (under Mode I loading and which is normal to the plane of the crack) and at a constant loading rate of 4 × 10−7 fs−1. The NVT ensemble was applied to maintain the system temperature.89,90 A simulation time step of 1 fs was chosen for all the simulations here. All MD simulations were carried out using large-scale atomic/molecular massively parallel simulator (LAMMPS)91 at 350 K. Visualization package OVITO was used to analyze and illustrate simulation results.92 The virial theorem93 was used to evaluate atomic stresses σxx, σyy, and σzz along corresponding directions. Normally, virial stress is equivalent to continuum Cauchy stress94 and can be utilized to represent the atomic stress at nanoscale.95 Virial stress consists of kinetic energy and potential energy components, and it can be defined as96

⎡ 1 1 σαβ = ⎢− ∑ mivαi vβi + Ω ⎢⎣ i 2

γMAI + γPbI = (Esurface − E bulk )/A

i

where Esurface and Ebulk denotes, respectively, the total potential energies of surface-exposed and bulk MAPbI3 models, A is the crosssectional surface area, and γMAI and γPbI2 are the surface energies of [MAI]0 and [PbI2]0 free surfaces, respectively.

ASSOCIATED CONTENT S Supporting Information *

Fαijrβij ⎥ ⎥⎦

j≠i

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsnano.6b05913. Supporting figures and table (PDF) Movie S1: grain rotation and grain-boundary sliding, which are colored by displacement vectors using program OVITO. The displacement vectors in the 6-grain model are denoted by orange arrows (AVI) Movie S2: grain rotation and grain-boundary sliding, which are colored by displacement vectors using program OVITO. The displacement vectors in the 10-grain model are denoted by orange arrows (AVI) Movie S3: grain rotation and grain-boundary sliding, which are colored by displacement vectors using program OVITO. The displacement vectors in the 14-grain model are denoted by orange arrows (AVI) Movie S4: crack propagation, colored by the local centrosymmetry parameter (CSP) map (AVI) Movie S5: nanovoid formation and migration, colored by the local centro-symmetry parameter (CSP) map using OVITO (AVI) Movie S6: von Mises stress distribution, colored using LAMMPS-computed local von Mises stress and our custom MATLAB program (AVI)

(8)

where α and β are the coordinate components (x, y, or z); Ω is the volume of the system; νiα and νiβ are the velocities of atom i along α and β directions, respectively; Fijα is the force between atoms i and j; and rijβ is the β component of the position of the atom i and j. Based on this virial stress definition, it is straightforward to carry out atomic stress decomposition by including only certain type of interatomic force (e.g., LJ, Buckingham, Coulombic, covalent bonding, and etc.) in eq 8 to obtain that specific stress contribution. In order to consider the Poisson’s effect (finite compressibility of a material) due to the NVT ensemble utilized here, the effective tensile stress along the z axis is computed as the following:97 σzz* = σzz − ν(σxx + σyy)

(9)

where ν = 0.23, 0.16, 0.18, and 0.31 are the Poisson’s ratios of the 6-, 10-, and 14-grain and defective single-crystalline MAPbI3, respectively, predicted using MD simulations as shown in Table 1. The von Mises stress is a popular engineering measure of the total effective stress for a 3D material system (considering all the normal and shear stresses) for deciding the yielding of a ductile material. It is defined based on the second deviatoric stress invariant, J2:98,99 σvM =

(12)

2



∑∑

(11)

i=1

(6)

Cubic:

⎛ 1 1 ⎞ = s11 − 2⎜s11 − s12 − s44⎟(l12l 2 2 + l 2 2l32 + l32l12) ⎝ E 2 ⎠

∑ |R⃗ i + R⃗ i + N /2|2

(σxx − σyy)2 + (σyy − σzz)2 + (σxx − σzz)2 + 6(σxy2 + σyz2 + σzx2 ) 2

(10)

AUTHOR INFORMATION

where σxx, σyy, and σzz are the normal stresses along x, y, and z-axes, respectively, and σxy, σyz, and σzx are the shear stresses. The mapping of local von Mises stress (with a grid size of 1.5 nm along the y-axis and 1 nm along the z-axis) distribution inside of a material can provide valuable information on local stress concentrations and help one understand the fracture process. The centro-symmetry parameter (CSP) was utilized to identify the local lattice disorder around an atom and characterize for most crystal defects (e.g., a dislocation or stacking fault). The CSP is close to zero for regular sites of a centro-symmetric crystal and becomes nonzero for defect atoms,100 and it can be defined as101

Corresponding Author

*E-mail: [email protected]. ORCID

Shangchao Lin: 0000-0002-6810-1380 Author Contributions †

These authors contributed equally to this work.

Notes

The authors declare no competing financial interest. 11054

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

(19) Feng, J. Mechanical Properties of Hybrid Organic-Inorganic CH3NH3BX3 (B= Sn, Pb; X= Br, I) Perovskites for Solar Cell Absorbers. APL Mater. 2014, 2, 081801. (20) Mattoni, A.; Filippetti, A.; Saba, M.; Delugas, P. Methylammonium Rotational Dynamics in Lead Halide Perovskite by Classical Molecular Dynamics: The Role of Temperature. J. Phys. Chem. C 2015, 119, 17421−17428. (21) Mattoni, A.; Filippetti, A.; Saba, M. I.; Caddeo, C.; Delugas, P. Temperature Evolution of Methylammonium Trihalide Vibrations at the Atomic Scale. J. Phys. Chem. Lett. 2016, 7, 529−535. (22) Gale, J. D. GULP: A Computer Program for the SymmetryAdapted Simulation of Solids. J. Chem. Soc., Faraday Trans. 1997, 93, 629−637. (23) Gale, J. D.; Rohl, A. L. The General Utility Lattice Program (GULP). Mol. Simul. 2003, 29, 291−341. (24) Hill, R. The Elastic Behaviour of a Crystalline Aggregate. Proc. Phys. Soc., London, Sect. A 1952, 65, 349−354. (25) Baikie, T.; Fang, Y.; Kadro, J. M.; Schreyer, M.; Wei, F.; Mhaisalkar, S. G.; Graetzel, M.; White, T. J. Synthesis and Crystal Chemistry of the Hybrid Perovskite (CH3NH3)PbI3 for Solid-State Sensitised Solar Cell Applications. J. Mater. Chem. A 2013, 1, 5628− 5641. (26) Sha, Z.; Quek, S.; Pei, Q.; Liu, Z.; Wang, T.; Shenoy, V.; Zhang, Y. Inverse Pseudo Hall-Petch Relation in Polycrystalline Graphene. Sci. Rep. 2014, 4, 5991. (27) Li, G.; An, Q.; Goddard, W. A.; Hanus, R.; Zhai, P.; Zhang, Q.; Snyder, G. J. Atomistic Explanation of Brittle Failure of Thermoelectric Skutterudite CoSb3. Acta Mater. 2016, 103, 775−780. (28) Shin, Y. A.; Yin, S.; Li, X.; Lee, S.; Moon, S.; Jeong, J.; Kwon, M.; Yoo, S. J.; Kim, Y.-M.; Zhang, T. Nanotwin-Governed Toughening Mechanism in Hierarchically Structured Biological Materials. Nat. Commun. 2016, 7, 10772. (29) Gunkelmann, N.; Bringa, E. M.; Kang, K.; Ackland, G. J.; Ruestes, C. J.; Urbassek, H. M. Polycrystalline Iron under Compression: Plasticity and Phase Transitions. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 144111. (30) Bringa, E. M.; Caro, A.; Wang, Y.; Victoria, M.; McNaney, J. M.; Remington, B. A.; Smith, R. F.; Torralva, B. R.; Van Swygenhoven, H. Ultrahigh Strength in Nanocrystalline Materials under Shock Loading. Science 2005, 309, 1838−1841. (31) Gumbsch, P.; Taeri-Baghbadrani, S.; Brunner, D.; Sigle, W.; Rühle, M. Plasticity and an Inverse Brittle-to-Ductile Transition in Strontium Titanate. Phys. Rev. Lett. 2001, 87, 085505. (32) Hossain, M.; Ali, M.; Parvin, F.; Islam, A. Mechanical and Optical Properties of Inverse-Perovskites Sc3InX (X= B, C, N). Comput. Mater. Sci. 2013, 73, 1−8. (33) Wang, B.; Yu, Y.; Lee, Y. J.; Bauchy, M. Intrinsic Nano-Ductility of Glasses: The Critical Role of Composition. Front. Mater. 2015, 2, 11. (34) Bauchy, M.; Laubie, H.; Qomi, M. A.; Hoover, C.; Ulm, F.-J.; Pellenq, R.-M. Fracture Toughness of Calcium−Silicate−Hydrate from Molecular Dynamics Simulations. J. Non-Cryst. Solids 2015, 419, 58− 64. (35) Huang, S.; Zhang, S.; Belytschko, T.; Terdalkar, S. S.; Zhu, T. Mechanics of Nanocrack: Fracture, Dislocation Emission, and Amorphization. J. Mech. Phys. Solids 2009, 57, 840−850. (36) Hirel, P.; Mrovec, M.; Elsässer, C. Atomistic Simulation Study of < 110> Dislocations in Strontium Titanate. Acta Mater. 2012, 60, 329−338. (37) Lynden-Bell, R. A Simulation Study of Induced Disorder, Failure and Fracture of Perfect Metal Crystals under Uniaxial Tension. J. Phys.: Condens. Matter 1995, 7, 4603−4624. (38) Romano, A.; Li, J.; Yip, S. Atomistic Simulation of Matter under Stress: Crossover from Hard to Soft Materials. Phys. A 2002, 304, 11− 22. (39) Li, G.; An, Q.; Goddard, W. A.; Hanus, R.; Zhai, P.; Zhang, Q.; Snyder, G. J. Atomistic Explanation of Brittle Failure of Thermoelectric Skutterudite CoSb3. Acta Mater. 2016, 103, 775−780.

ACKNOWLEDGMENTS The authors would like to acknowledge the startup funding from the Energy and Materials Initiative from the Florida State University. The authors would like to thank Professor Z. Yu from the Florida State University for helpful discussions. S. Lin also would like to acknowledge the First Year Assistant Professor (FYAP) award from the Florida State University. REFERENCES (1) Stranks, S. D.; Snaith, H. J. Metal-halide Perovskites for Photovoltaic and Light-Emitting Devices. Nat. Nanotechnol. 2015, 10, 391−402. (2) He, Y.; Galli, G. Perovskites for Solar Thermoelectric Applications: A First Principle Study of CH3NH3AI3 (A = Pb and Sn). Chem. Mater. 2014, 26, 5394−5400. (3) Mettan, X.; Pisoni, R.; Matus, P.; Pisoni, A.; Jaćimović, J.; Náfrádi, B.; Spina, M.; Pavuna, D.; Forró, L.; Horváth, E. Tuning of the Thermoelectric Figure of Merit of CH3NH3MI3 (M = Pb,Sn) Photovoltaic Perovskites. J. Phys. Chem. C 2015, 119, 11506−11510. (4) Yin, W.-J.; Shi, T.; Yan, Y. Unusual Defect Physics in CH3NH3PbI3 Perovskite Solar Cell Absorber. Appl. Phys. Lett. 2014, 104, 063903. (5) Mosconi, E.; Umari, P.; De Angelis, F. Electronic and Optical Properties of Mixed Sn−Pb Organohalide Perovskites: A First Principles Investigation. J. Mater. Chem. A 2015, 3, 9208−9215. (6) Haruyama, J.; Sodeyama, K.; Han, L.; Tateyama, Y. Termination Dependence of Tetragonal CH3NH3PbI3 Surfaces for Perovskite Solar Cells. J. Phys. Chem. Lett. 2014, 5, 2903−2909. (7) Mosconi, E.; Azpiroz, J. M.; De Angelis, F. Ab Initio Molecular Dynamics Simulations of Methylammonium Lead Iodide Perovskite Degradation by Water. Chem. Mater. 2015, 27, 4885−4892. (8) Grätzel, M. The Light and Shade of Perovskite Solar Cells. Nat. Mater. 2014, 13, 838−842. (9) Xing, G.; Mathews, N.; Sun, S.; Lim, S. S.; Lam, Y. M.; Grätzel, M.; Mhaisalkar, S.; Sum, T. C. Long-range Balanced Electron-and Hole-Transport Lengths in Organic-Inorganic CH3NH3PbI3. Science 2013, 342, 344−347. (10) Liu, M.; Johnston, M. B.; Snaith, H. J. Efficient Planar Heterojunction Perovskite Solar Cells by Vapour Deposition. Nature 2013, 501, 395−398. (11) Rupnowski, P.; Sopori, B. Strength of Silicon Wafers: Fracture Mechanics Approach. Int. J. Fract. 2009, 155, 67−74. (12) Buehler, M. J.; van Duin, A. C. T.; Goddard, W. A. Multiparadigm Modeling of Dynamical Crack Propagation in Silicon Using a Reactive Force Field. Phys. Rev. Lett. 2006, 96, 095505. (13) Bretschneider, S. A.; Weickert, J.; Dorman, J. A.; SchmidtMende, L. Research Update: Physical and Electrical Characteristics of Lead Halide Perovskites for Solar Cell Applications. APL Mater. 2014, 2, 040701. (14) Ong, K. P.; Goh, T. W.; Xu, Q.; Huan, A. Mechanical Origin of the Structural Phase Transition in Methylammonium Lead Iodide CH3NH3PbI3. J. Phys. Chem. Lett. 2015, 6, 681−685. (15) Poorkazem, K.; Liu, D.; Kelly, T. L. Fatigue Resistance of a Flexible, Efficient, and Metal Oxide-Free Perovskite Solar Cell. J. Mater. Chem. A 2015, 3, 9241−9248. (16) Park, M.; Kim, H. J.; Jeong, I.; Lee, J.; Lee, H.; Son, H. J.; Kim, D. E.; Ko, M. J. Mechanically Recoverable and Highly Efficient Perovskite Solar Cells: Investigation of Intrinsic Flexibility of Organic−Inorganic Perovskite. Adv. Energy Mater. 2015, 5, 1501406. (17) Sun, S.; Fang, Y.; Kieslich, G.; White, T. J.; Cheetham, A. K. Mechanical Properties of Organic−Inorganic Halide Perovskites, CH3NH3PbX3 (X= I, Br and Cl), by Nanoindentation. J. Mater. Chem. A 2015, 3, 18450−18455. (18) Rakita, Y.; Cohen, S. R.; Kedem, N. K.; Hodes, G.; Cahen, D. Mechanical Properties of APbX3 (A= Cs or CH3NH3; X= I or Br) Perovskite Single Crystals. MRS Commun. 2015, 5, 623−629. 11055

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano

(64) Zhao, L.-D.; Zhang, B.-P.; Li, J.-F.; Zhou, M.; Liu, W.-S.; Liu, J. Thermoelectric and Mechanical Properties of Nano-SiC-Dispersed Bi2Te3 Fabricated by Mechanical Alloying and Spark Plasma Sintering. J. Alloys Compd. 2008, 455, 259−264. (65) Tasi, C.-H.; Tseng, Y.-C.; Jian, S.-R.; Liao, Y.-Y.; Lin, C.-M.; Yang, P.-F.; Chen, D.-L.; Chen, H.-J.; Luo, C.-W.; Juang, J.-Y. Nanomechanical Properties of Bi2Te3 Thin Films by Nanoindentation. J. Alloys Compd. 2015, 619, 834−838. (66) Ueno, K.; Yamamoto, A.; Noguchi, T.; Inoue, T.; Sodeoka, S.; Obara, H. Optimization of Hot-Press Conditions of Zn4Sb3 for High Thermoelectric Performance. J. Alloys Compd. 2005, 392, 295−299. (67) Freiman, S.; Mecholsky, J., Jr The Fracture Energy of Brittle Crystals. J. Mater. Sci. 2010, 45, 4063−4066. (68) Evans, A.; Johnson, H. A Fracture Mechanics Study of ZnSe for Laser Window Applications. J. Am. Ceram. Soc. 1975, 58, 244−249. (69) Bilal, M.; Shafiq, M.; Ahmad, I.; Khan, I. First Principle Studies of Structural, Elastic, Electronic and Optical Properties of ZnChalcogenides under Pressure. J. Semicond. 2014, 35, 072001. (70) Brinkevich, D.; Vabishchevich, N.; Prosolovich, V. Micromechanical Properties of GaP Dy Epilayers. Inorg. Mater. 2012, 48, 768−772. (71) Jochum, M.; Werner-Zwanziger, U.; Zwanziger, J. Observable Effects of Mechanical Stress Induced by Sample Spinning in Solid State Nuclear Magnetic Resonance. J. Chem. Phys. 2008, 128, 052304. (72) Fujikane, M.; Leszczyński, M.; Nagao, S.; Nakayama, T.; Yamanaka, S.; Niihara, K.; Nowak, R. Elastic−Plastic Transition during Nanoindentation in Bulk GaN Crystal. J. Alloys Compd. 2008, 450, 405−411. (73) Ratschinski, I.; Leipner, H.; Heyroth, F.; Fränzel, W.; Hammer, R.; Jurisch, M. Dislocations and Cracks at Vickers Indentations in (0 0 0 1) GaN Single Crystals. Philos. Mag. Lett. 2010, 90, 565−571. (74) Deger, C.; Born, E.; Angerer, H.; Ambacher, O.; Stutzmann, M.; Hornsteiner, J.; Riha, E.; Fischerauer, G. Sound Velocity of AlxGa1‑xN Thin Films Obtained by Surface Acoustic-Wave Measurements. Appl. Phys. Lett. 1998, 72, 2400. (75) Jian, S.-R.; Fang, T.-H.; Chuu, D.-S. Nanomechanical Characterizations of InGaN Thin Films. Appl. Surf. Sci. 2006, 252, 3033−3042. (76) Jian, S.-R.; Jang, J. S.-C.; Lai, Y.-S.; Yang, P.-F.; Yang, C.-S.; Wen, H.-C.; Tsai, C.-H. Mechanical Properties of InGaN Thin Films Deposited by Metal-Organic Chemical Vapor Deposition. Mater. Chem. Phys. 2008, 109, 360−364. (77) Brostow, W.; Dussault, J.-P.; Fox, B. L. Construction of Voronoi Polyhedra. J. Comput. Phys. 1978, 29, 81−92. (78) Mattoni, A.; Ippolito, M.; Colombo, L. Atomistic Modeling of Brittleness in Covalent Materials. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 224103. (79) Bitzek, E.; Kermode, J. R.; Gumbsch, P. Atomistic Aspects of Fracture. Int. J. Fract. 2015, 191, 13−30. (80) Delugas, P.; Caddeo, C.; Filippetti, A.; Mattoni, A. Thermally Activated Point-Defects Diffusion in Methylammonium Lead Trihalide: Anisotropic and Ultra-High Mobility of Iodine. J. Phys. Chem. Lett. 2016, 7, 2356−2361. (81) Wang, M.; Lin, S. Anisotropic and Ultralow Phonon Thermal Transport in Organic−Inorganic Hybrid Perovskites: Atomistic Insights into Solar Cell Thermal Management and Thermoelectric Energy Conversion Efficiency. Adv. Funct. Mater. 2016, 26, 5297− 5306. (82) Caddeo, C.; Melis, C.; Saba, M. I.; Filippetti, A.; Colombo, L.; Mattoni, A. Tuning the Thermal Conductivity of Methylammonium Lead Halide by the Molecular Substructure. Phys. Chem. Chem. Phys. 2016, 18, 24318−24324. (83) Liu, S.; Cohen, R. E. Response of Methylammonium Lead Iodide to External Stimuli and Caloric Effects from Molecular Dynamics Simulations. J. Phys. Chem. C 2016, 120, 17274−17281. (84) Delugas, P.; Filippetti, A.; Mattoni, A. Methylammonium Fragmentation in Amines as Source of Localized Trap Levels and the Healing Role of Cl in Hybrid Lead-Iodide Perovskites. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 92, 045301.

(40) Brochard, L.; Hantal, G.; Laubie, H.; Ulm, F.-J.; Pellenq, R. J. Capturing Material Toughness by Molecular Simulation: Accounting for Large Yielding Effects and Limits. Int. J. Fract. 2015, 194, 149−167. (41) Rice, J. R. Mathematical Analysis in the Mechanics of Fracture. In Fracture: An Advanced Treatise; Liebowitz, H., Ed.; Academic Press: New York, 1968; pp 191−311. (42) Šmilauer, V.; Hoover, C. G.; Bažant, Z. P.; Caner, F. C.; Waas, A. M.; Shahwan, K. W. Multiscale Simulation of Fracture of Braided Composites via Repetitive Unit Cells. Eng. Fract. Mech. 2011, 78, 901− 918. (43) Hillerborg, A. The Theoretical Basis of a Method to Determine the Fracture Energy GF of Concrete. Mater. Constr. 1985, 18, 291− 296. (44) Griffith, A. A. The Phenomena of Rupture and Flow in Solids. Philos. Trans. R. Soc., A 1921, 221, 163−198. (45) Bocquet, M.-L.; Coasne, B.; Pellenq, R.; Ulm, F.-J. Organic− Clay Interfacial Chemical Bonds Probed by ab Initio Calculations. J. Phys. Chem. C 2015, 119, 6511−6517. (46) Irwin, G. R. Analysis of Stresses and Strains Near the End of a Crack Traversing a Plate. J. Appl. Mech. 1997, 24, 361−364. (47) Orlovskaya, N.; Kleveland, K.; Grande, T.; Einarsrud, M.-A. Mechanical Properties of LaCoO3 Based Ceramics. J. Eur. Ceram. Soc. 2000, 20, 51−56. (48) Gao, L.; Jiang, L.; Sun, J. Carbon Nanotube-Ceramic Composites. J. Electroceram. 2006, 17, 51−55. (49) Blamey, J. M.; Parry, T. V. The Effect of Processing Variables on the Mechanical and Electrical Properties of Barium Titanate PositiveTemperature-Coefficient-of-Resistance Ceramics. J. Mater. Sci. 1993, 28, 4317−4324. (50) Bernstein, N.; Hess, D. Lattice Trapping Barriers to Brittle Fracture. Phys. Rev. Lett. 2003, 91, 025501. (51) Mattoni, A.; Colombo, L.; Cleri, F. Atomic Scale Origin of Crack Resistance in Brittle Fracture. Phys. Rev. Lett. 2005, 95, 115501. (52) Bade, S. G. R.; Li, J.; Shan, X.; Ling, Y.; Tian, Y.; Dilbeck, T.; Besara, T.; Geske, T.; Gao, H.; Ma, B.; et al. Fully Printed Halide Perovskite Light-Emitting Diodes with Silver Nanowire Electrodes. ACS Nano 2016, 10, 1795−1801. (53) Chen, C.; Leipold, M. Fracture Toughness of Silicon. Am. Ceram. Soc. Bull. 1980, 59, 469−472. (54) Wortman, J.; Evans, R. Young’s Modulus, Shear Modulus, and Poisson’s Ratio in Silicon and Germanium. J. Appl. Phys. 1965, 36, 153−156. (55) Ericson, F.; Schweitz, J. Å. Micromechanical Fracture Strength of Silicon. J. Appl. Phys. 1990, 68, 5840−5844. (56) Ericson, F.; Johansson, S.; Schweitz, J.-Å. Hardness and Fracture Toughness of Semiconducting Materials Studied by Indentation and Erosion Techniques. Mater. Sci. Eng., A 1988, 105-106, 131−141. (57) White, G.; Freiman, S.; Fuller, E.; Baker, T. Effects of Crystal Bonding on Brittle Fracture. J. Mater. Res. 1988, 3, 491−497. (58) Hjort, K.; Soderkvist, J.; Schweitz, J.-A. Gallium Arsenide as a Mechanical Material. J. Micromech. Microeng. 1994, 4, 1−13. (59) Wang, Y. B.; Wang, L. F.; Joyce, H. J.; Gao, Q.; Liao, X. Z.; Mai, Y. W.; Tan, H. H.; Zou, J.; Ringer, S. P.; Gao, H. J. Super Deformability and Young’s Modulus of GaAs Nanowires. Adv. Mater. 2011, 23, 1356−1360. (60) Margevicius, R.; Gumbsch, P. Influence of Crack Propagation Direction on {110} Fracture Toughness of Gallium Arsenide. Philos. Mag. A 1998, 78, 567−581. (61) Adachi, S. GaAs and Related Materials: Bulk Semiconducting and Superlattice Properties; World Scientific: Singapore, 1994; pp 59−60. (62) Sweet, C. A.; Schulte, K. L.; Simon, J. D.; Steiner, M. A.; Jain, N.; Young, D. L.; Ptak, A. J.; Packard, C. E. Controlled Exfoliation of (100) Gaas-Based Devices by Spalling Fracture. Appl. Phys. Lett. 2016, 108, 011906. (63) Lamuta, C.; Cupolillo, A.; Politano, A.; Aliev, Z. S.; Babanly, M. B.; Chulkov, E. V.; Pagnotta, L. Indentation Fracture Toughness of Single-Crystal Bi2Te3 Topological Insulators. Nano Res. 2016, 9, 1032−1042. 11056

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057

Article

ACS Nano (85) Dualeh, A.; Gao, P.; Seok, S. I.; Nazeeruddin, M. K.; Grätzel, M. Thermal Behavior of Methylammonium Lead-Trihalide Perovskite Photovoltaic Light Harvesters. Chem. Mater. 2014, 26, 6160−6164. (86) Bresme, F.; Hafskjold, B.; Wold, I. Nonequilibrium Molecular Dynamics Study of Heat Conduction in Ionic Systems. J. Phys. Chem. 1996, 100, 1879−1888. (87) Nye, J. F. Physical Properties of Crystals: Their Representation by Tensors and Matrices; Oxford University Press Inc.: New York, 1985; pp 142−145. (88) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (89) Nosé, S. A Unified Formulation of The Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511−519. (90) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. (91) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (92) Stukowski, A. Visualization and Analysis of Atomistic Simulation Data with OVITO−The Open Visualization Tool. Modell. Simul. Mater. Sci. Eng. 2010, 18, 015012. (93) Tsai, D. The Virial Theorem and Stress Calculation in Molecular Dynamics. J. Chem. Phys. 1979, 70, 1375−1382. (94) Cui, C.; Beom, H. Molecular Dynamics Simulations of Edge Cracks in Copper and Aluminum Single Crystals. Mater. Sci. Eng., A 2014, 609, 102−109. (95) Diao, J.; Gall, K.; Dunn, M. L. Atomistic Simulation of the Structure and Elastic Properties of Gold Nanowires. J. Mech. Phys. Solids 2004, 52, 1935−1962. (96) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989; pp 60−62. (97) Andreev, G. E. Brittle Failure of Rock Materials; Balkema: Rotterdam, 1995; pp 112−115. (98) Wang, D.; Lee, J.; Holland, K.; Bibby, T.; Beaudoin, S.; Cale, T. Von Mises Stress in Chemical-Mechanical Polishing Processes. J. Electrochem. Soc. 1997, 144, 1121−1127. (99) Botari, T.; Perim, E.; Autreto, P.; van Duin, A.; Paupitz, R.; Galvao, D. Mechanical Properties and Fracture Dynamics of Silicene Membranes. Phys. Chem. Chem. Phys. 2014, 16, 19417−19423. (100) Stukowski, A. Structure Identification Methods for Atomistic Simulations of Crystalline Materials. Modell. Simul. Mater. Sci. Eng. 2012, 20, 045021. (101) Kelchner, C. L.; Plimpton, S.; Hamilton, J. Dislocation Nucleation and Defect Structure during Surface Indentation. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 11085−11088.

11057

DOI: 10.1021/acsnano.6b05913 ACS Nano 2016, 10, 11044−11057