Bending energy of 2D materials: graphene, MoS₂ and imogolite

Abstract

The bending process of 2D materials, subject to an external force, is investigated, and applied to graphene, molybdenum disulphide (MoS₂), and imogolite. For graphene we obtained 3.43 eV Ų per atom for the bending modulus, which is in good agreement with the literature. We found that MoS₂ is ∼11 times harder to bend than graphene, and has a bandgap variation of ∼1 eV as a function of curvature. Finally, we also used this strategy to study aluminosilicate nanotubes (imogolite) which, in contrast to graphene and MoS₂, present an energy minimum for a finite curvature radius. Roof tile shaped imogolite precursors turn out to be stable, and thus are expected to be created during imogolite synthesis, as predicted to occur by self-assembly theory.

---

1 Introduction

One of the most exciting possibilities of nanotechnology is to control electronic properties at the atomic level, since new degrees of freedom become available, as compared to bulk materials. Among them, applying out-of-plane strain is a very versatile way to engineer several properties. In two-dimensional (2D) van der Waals materials, the strain is very anisotropic, and the folding and bending has recently attracted particular attention.^1–5 Even the possibility of nano-kirigami has been explored, and stretchable transistors have been fabricated.^6,7

Among the systems where the bending of a graphene sheet is relevant, an example is the rippling of suspended graphene with mesoscopic amplitudes and wavelengths.⁸ Also, supported graphene monolayers have been intensely studied.⁹ The effect of different metallic substrates has been investigated as well, showing that roughness can be tuned by varying the support.¹⁰ Recently, Zhou et al.¹¹ reported on the interaction with a substrate where graphene forms wrinkles and also adopts different conformations. Theoretical studies have shown that the bending stiffness of a graphene monolayer is critical to achieve structural stability, and it determines the morphology for both suspended and supported graphene sheets.¹² Moreover, it could have significant effects on their electronic properties.

Other 2D materials like h-BN, MoS₂ or WSe₂ have also been the focus of attention, since the control of the roughness, or the use of wrinkles, favors some applications. Moreover, new stable 2D materials and their hetero-structures¹³ are being synthesized or proposed at a fast pace. Therefore it is highly desirable to have simple models that are able to describe their basic properties. However, models need parameters in order to distinguish among the different materials.

The main objective of this work is to implement a strategy to calculate the energy that is needed to bend 2D materials (i.e. to calculate the bending modulus). This strategy does not rely on a specific theory level or code. It uses an external force, specifically a harmonic potential with the desired geometry, that imposes a specific energy minimum. While we use a cylindrical geometry, other configurations such as a horseshoe are also possible. Once the harmonic potential is turned on, the system is allowed to relax and to adopt the desired shape, but it remains free of any additional constraint or strain.

Other strategies to study the bending response of 2D and layered materials have been put forward. For example, the bending energy of graphene has been studied for carbon nanotubes of different radii.^14–17 Recently Xiong et al. computed the response of a MoS₂ single layer¹⁸ with a method that applied external forces to the edges of a MoS₂ sheet. Another approach uses a bent carbon template to control the curvature of aluminosilicates, such as a montmorillonite layer.¹⁹

We start with a simple case – graphene – which allows us to test our method against available data.^14–16 Next, we apply the method to a MoS₂ layer; that is, we generalize to a two atomic species system three atoms thick, coupling classical potentials with DFT, in order to study the effects of bending on the electronic structure, at a low computational cost. Finally, we study the bending of imogolite, a versatile aluminosilicate material that spontaneously adopts a nanotube structure. Imogolite has four atomic species, and is ≈0.5 nm thick. In contrast to graphene and MoS₂, the bending energy of imogolite has a minimum as a function of diameter, and is highly monodisperse (i.e. quite uniform in diameter).^20–24

Two hypotheses on the formation of imogolite have been put forward, and they are subtly different; however, in both the bending energy plays a key role. The hypothesis by Maillet et al.²⁵ models the growth of imogolite in two stages. During the first one, the transformation of roof tile shaped protoimogolite into short nanotubes takes place. In the second stage the tip–tip coupling of these short nanotubes allows for the nanotube growth. The other hypothesis was presented by Yucelen et al.²⁶ and suggests the self-assembly of protoimogolite to form the nanotubes. While both theories imply the occurrence of protoimogolite, only the former considers an intermediate short nanotube stage. Here we show that curved pieces of imogolite could exist and that they are energetically stable, supporting the plausibility of both hypotheses.

2 Methods and model

2.1 Bending model

A 2D sheet subject to a cylindrically symmetrical harmonic external potential adopts a curvature imposed by the external force, as illustrated in Fig. 1a. In this paper we study three different systems: graphene, MoS₂ and imogolite sheets. Also, to show that this strategy is general enough, we carried out density functional theory (DFT) calculations, as implemented in the VASP code,^27–30 for graphene and MoS₂.

Fig. 1 General procedure for our simulations. (a) Starting with a flat ribbon (e.g. graphene), an external potential is applied to smoothly generate the desired shape. The only strain source is the external force. (b) Strain energy ΔE(κ) for several curvature radii. The inset shows the abrupt energy change as the nanotube closes. (c) All of the carbon bonds are colored according to their length, showing that the 2D structure is maintained during the whole bending process.

Our initial graphene configuration is a flat ribbon, centered at the origin (which corresponds to the harmonic potential being set equal to zero, as illustrated in Fig. 1a); next we allow the system to relax. Thereafter, we adopt (in small steps) finite values for the curvature radius, and relax again; this procedure is iterated until the smallest possible R is reached (i.e. the nanotube geometry). Throughout this paper we use periodic boundary conditions along the long edge of the ribbon in order to simulate an infinitely long system. During this process the 2D structure of graphene is preserved without inducing strain on the system. In fact, in Fig. 1c we show that during the bending process the bond length of the carbon atoms is kept around the equilibrium distance.

The energy associated to the external force is

where K_ext = 250 eV Å⁻¹, R is the curvature radius, and R_i is the distance of each atom, subject to the external force, to the center of curvature. On the other hand, the strain energy of the system, as a function of the curvature R⁻¹ = κ, is given bywhere E is the total potential energy of the system (ignoring the external energy, E_ext), N is the number of atoms and κ is the curvature imposed by the external force (eqn (1), see Fig. 1b). The bending modulus¹⁶C_b is defined as the second derivative of the strain energy with respect to the curvature κ

2.2 Computational methods

In the study of the bending process we used classical molecular dynamics simulations, as implemented in the large-scale atomic/molecular massively parallel simulator (LAMMPS) code.³¹ The bending due to the external force was combined with energy optimization strategies. We used the fast inertial relaxation engine³² (FIRE) together with conjugate gradient minimization. We chose this combination because the FIRE is able to avoid getting stuck in local minimum energy configurations, while the conjugate gradient is faster.

The graphene C–C atomic interactions were modeled by the ReaxFF force field.³³ The ReaxFF parameters for the C–C interactions were first generated by Chenoweth et al.³⁴ and later improved by Srinivasan et al.³⁵ This force field has been tested for different carbon allotropes such as fullerenes, graphene, carbon nanotubes, amorphous carbon and diamond, with very good results.^35–38 The MoS₂ sheet was modeled with a second-generation reactive bond-order formalism by Liang et al.,^39,40 and has been used to investigate the mechanical properties of MoS₂ nanotubes,⁴¹ nano-indentation of crystalline MoS₂,⁴² for a few layers of MoS₂ under tension,⁴³ and for the properties of a MoS₂ layer.^18,44–47 Finally, in order to model imogolite atoms we used the CLAYFF potential,⁴⁸ since it has been proven to be adequate to model this aluminosilicate nanotube.^20,21,49,50 The CLAYFF potential, developed by Cygan et al.,⁴⁸ incorporates the charges of every single atom, the van der Waals interactions, a harmonic potential for the O–H group stretching, and harmonic stretching for the angles between Al–O–H and Si–O–H.

The band structure calculations were obtained with the VASP package.^27–30 PAW pseudopotentials^51,52 were used and the PBE exchange correlation⁵³ was employed. The energy cutoff was set to 400 eV for graphene and to 300 eV for MoS₂. We used 20 K-points along the only periodic direction. At least 15 Å of the vacuum were included in the non-periodic directions. The open visualization tool (OVITO) was used⁵⁴ for the graphics and for the post-processing of our simulations, and PyProcar was employed to perform the band structure analysis.⁵⁵

3 Results and discussion

3.1 Graphene

The initial configuration for our simulations is a graphene ribbon with periodic boundary conditions along the bending axis (i.e. perpendicular to the external forces) which is 8 nm wide. In Fig. 1a we illustrate the general procedure for our simulations of the ribbon bending, with dangling bonds and zig-zag termination. We represent the minimum of the external potential imposed for each curvature as a red dashed line. As shown in Fig. 1b the energy increases as a function of curvature κ, as expected. This behavior of the graphene strain energy is well known. Since a graphene sheet has dangling bonds at the edges there is an abrupt decrease in the total energy when the system finally adopts the nanotube geometry (κ ≈ 0.78 nm⁻¹). To quantify the effect of the external force on the ribbon geometry in Fig. 1c we plot the pair correlation function g(r) for the three stages of the bending process. These stages are: the beginning, where the sheet is completely flat; an intermediate curvature; and finally as a CNT, when we turn off the external force and allow the system to relax. During the whole process the g(r) peaks associated to the first, second and third neighbors are located at 1.42, 2.46 and 2.84 Å, respectively. This constitutes evidence that the external force preserves the 2-D nature of the system during the whole process. Some smaller peaks, observed during the bending stage, are associated with the relaxation of the free edges, and disappear when the CNT geometry is achieved, when all the C atoms are sp² coordinated.

In Fig. 1b we fitted a quadratic curve, of the form , for the 8 nm wide ribbon, obtaining a value for the bending modulus of 34.2 meV nm² per atom = 3.42 eV Ų per atom. However, for a narrow width ribbon finite size effects can be important, and have to be included. In fact, in Fig. 2 we observe that for a 0.7 nm ribbon the C_b is around 65% of the 8 nm value; but, for a 4 nm ribbon the C_b is already ≈95% of the 8 nm value. In other words, ribbons less than ∼2.5 nm wide are substantially easier to fold.

Fig. 2 Bending modulus of a graphene ribbon as a function of its width d. While the usual methods, such as indentation or CNT bending energy, ignore finite size and edge effects, our method captures them. On the other hand, in the large diameter limit, where edges are negligible, these effects can be ignored. In the inset, we show a linear regression of the bending modulus as a function of 1/d (solid red line). From this linear regression we estimate a value of 3.43 eV Ų per atom for graphene (1/d → 0 limit).

To obtain the C_b value for graphene we adopt the result obtained by the linear regression shown as an inset in Fig. 2, which corresponds to 3.43 eV Ų per atom. This value is closer to the DFT values^14,15 of 3.9 and 4.3 eV Ų per atom than the classically calculated ones of 1.8 and 2.2 eV Ų per atom.¹⁶

3.1.1 Changes in electronic properties due to bending. To study the effects of bending acting upon a graphene nanoribbon (GNR), we employed DFT. We found an almost undistorted electronic structure of a bent zig-zag GNR, even for a very small curvature radius. In Fig. 3 we calculate the band structure of a zig-zag GNR with a thickness of 2.4 nm. The pristine and bent cases, Fig. 3a and b, look almost identical: a textbook projection of the graphene’s Brillouin zone on the GNR reciprocal space (the sub-bands already resemble graphene’s Dirac cone). The valence and conduction bands are degenerated from ,⁵⁶ due to the breaking of the sub-lattice symmetry at the boundaries. Additionally, there are flat bands, typical of the dangling bonds of terminal C atoms (they would disappear after hydrogenation). These dangling bonds are somewhat hybridized with bulk states for the bent GNR, which is the only noticeable difference compared to a flat GNR.

Fig. 3 Band structure of a flat GNR (a) and a bent GNR, with the curvature radius R = 0.5 nm (b). The reddish bands are derived from the edges of the ribbon, the blue bands do not have a projection on the edges and intermediate colors apply for mixed states.

3.2 MoS₂

Increasing the complexity of the 2D material under study, we now investigate the bending of a MoS₂ ribbon. In this case we apply the bending force on the central Mo layer, while we allow the S atoms to move freely. It is worth remarking that the Mo atoms are able to move orthogonally to the constraining force, thus avoiding unwanted strain.

To compare the bending response of a MoS₂ monolayer with the graphene case, we compute the bending stiffness⁴⁷D, which is related to the bending modulus by D = NC_b/A, where A is the transverse area of the ribbon. This way, the bending stiffness D of an 8 nm MoS₂ ribbon turns out to be more than 11 times larger than that of a graphene ribbon of the same width. Additionally, we estimate the Föppl–von-Kármán number⁶ of a MoS₂ monolayer to be about 40 times smaller than that of graphene. The Föppl–von-Kármán number is defined as E/D, where E corresponds to the Young modulus.^18,38 Inspection of Fig. 4b of the paper by Blees et al.⁶ suggests that displacement, due to a laser, is negligible (less than the measurement noise). This argument does not rule out the feasibility of building nano-kirigami or foldable transistors^6,7 but is an indication that the activation process would have to be different.

In a recent article Xiong and Cao compared DFT and molecular mechanics results for some MoS₂ force fields.⁴⁷ For the molecular mechanics they reported a bending stiffness of 8.7–13.4 eV depending on the force field used. On the other hand, they compared this with previous DFT results from Xiao et al. based on MoS₂ single wall nanotube calculations.⁵⁷ They have shown that the force field used in the present work slightly overestimates the bending stiffness with respect to the DFT results. Finally, using the values shown in Fig. 4, with a linear regression for 1/d → 0 limit, we obtain a bending stiffness for the MoS₂ bulk layer of D = 16 eV. Casillas et al.⁵⁸ reported a bending stiffness value of 4.1–13.24 eV, based on continuum theory. Our results for a MoS₂ single layer are slightly over the values reported by other computational calculations and experiments.

Fig. 4 Bending modulus of a MoS₂ ribbon as a function of its width d. MO: blue; S: red.

In relation to the MoS₂ finite width effects, in contrast with graphene ribbons, the difference between a narrow (0.95 nm) and a wide (8 nm) MoS₂ monolayer is more pronounced. As can be seen in Fig. 4 the bending coefficient of our narrowest ribbon is only 42% of the value obtained for large sheets (d ≥ 8 nm).

3.2.1 Electronic property changes due to bending. To obtain some insight into the changes in the electronic properties of a MoS₂ ribbon due to curvature we employ first principles calculations. To alleviate the computational cost we used structures optimized classically (see previous section), but allowed the S atoms to relax; also, the strain in the periodic direction was minimized.

The (armchair) ribbon edges were not passivated, leaving dangling bonds. Eigenvalue analysis shows that these states remain highly localized at the edge of the ribbon – as expected of an insulating material – leaving bulk-like states in the inner region of the ribbon; therefore the exact details at the ribbon edges are not an issue. We will elaborate on this topic later on.

The effect of bending on the MoS₂ band structure is shown in Fig. 5 for a flat and finite curvature radius system. It is evident that the inner, bulk-like bands (the red lines in Fig. 5) are qualitatively similar for R = 8.12 Å and R → ∞, but the gap between the valence and conduction bands is ∼0.8 eV smaller. Instead, the electronic states at the edge of the nanoribbon (blue in Fig. 5) suffer larger changes. This finding suggests that the bending effects on the electronic structure are rather modest; in fact, a bent MoS₂ layer still behaves like flat MoS₂ but with a smaller band gap, even for large curvatures. The surface states are likely to be more affected by the bending, but they remain localized at the edges. To quantify the extent of the change of the band-structure due to bending, Fig. 6 shows the band gap at two representative points (Γ and M) as a function of the curvature radius. Moderate or small bending results in small gap changes, while a tiny curvature radius (7 Å) can reduce the gap by ∼60%. We do not expect a metallic state to develop, or band inversion to occur, even for smaller radii. However, a more accurate calculation is required to obtain a precise quantitative picture.

Fig. 5 Band structure of an armchair MoS₂ nanoribbon. The upper and lower panels correspond to a flat (R → ∞) and a bent nanoribbon, respectively. The color code shows the (real-space) localization of the eigenvalues: red is used for the inner subbands (i.e. bulk-like), while the edge subbands are blue; intermediate colors correspond to the hybridization of bulk and edge states.

Fig. 6 Electronic band gap at the Γ and M points as a function of curvature. The system band gap is always located at Γ; the band gap at M is just a parameter to quantify bending effects. Surface states were ignored in the gap measurement.

3.3 Imogolite

Finally we report our results for the strain energy of a curved imogolite sheet. In Fig. 7 we provide on the left hand side an axial view of the imogolite tube, and on the right a view of its basic unit [(OH)₃Al₂O₃SiOH]₂. The single walled imogolite nanotube has the chemical composition [(OH)₃Al₂O₃SiOH]_2Nθ; this unit is repeated angularly, around the longitudinal axis, N_θ times, and N_Z times along the longitudinal direction.⁵⁹ In Fig. 8 we show the energy of an imogolite sheet with this composition and N_θ = 5, N_Z = 6, assuming periodic boundary conditions along the longitudinal direction. The bending force is applied on the Al atoms. In Fig. 8 we include three insets of the sheet configuration: the initial planar sheet; the minimum energy configuration (half cylinder); and the structure obtained for the smallest curvature radius we calculated. Previously^21,49 we showed that the energy minimum for an imogolite nanotube, using the CLAYFF force field, occurs for N_θ = 10. In the present case, a ribbon of half the width (N_θ = 5) reaches its minimum energy for a half-cylinder shape. We noticed that our results for the semi-cylinder indicate that it is energetically favorable as compared to forcing it to close upon itself to form a full cylinder. This reinforces the growth model proposed by Maillet et al.,²⁵ suggesting that during the growth of imogolite, curved protoimogolite appear and assemble. The shallow minimum in Fig. 8 constitutes an indication that imogolite nanotubes of N_θ < 8 are energetically very unfavorable, due to the excess bending energy that is required. Moreover, this shallowness of the minimum also suggests that during the formation process, depending on the synthesis conditions,²⁶ nanotubes 9 ≤ N_θ ≤ 15 can form, since the energy required is quite small.

Fig. 7 Structure of an imogolite nanotube. (a) An axial view of imogolite for N_θ = 10. The basic unit is marked by the continuous black line. (b) A side view of the basic imogolite unit [(OH)₃Al₂O₃SiOH]₂, that is repeated N_θ times to form the nanotube. Green Al; red O; yellow Si; and light gray H.

Fig. 8 Energy vs. curvature of an imogolite sheet. The force is applied on the Al atoms. A minimum of the bending energy occurs for 9 < N_θ < 10. The imogolite structure inset corresponds to N_θ = 5 angular repetitions. The dashed vertical lines, and the numbers that label them, correspond to the curvature radius of a completely closed nanotube with that N_θ value. Similar data, for imogolite sheets for N_θ = 4 and 6, are included as ESI.† Green Al; red O; yellow Si; and light gray H.

4 Summary and conclusions

We studied the bending process of nano-sheets, by modeling and quantitatively treating the rolling of 2D materials into tubular conformations. To display the versatility and the wide range of applicability of the procedure we studied three different systems: graphene, MoS₂ and imogolite. Moreover, our model is straightforward to generalize to other types of nano-sheets,⁶⁰ like h-BN, MoSe₂ and TiO₂. It can also be applied to nano-sheets that roll around different kinds of nano-cylinder^61,62 such as TiO₂ nanotubes, metallic nanowires, CNTs (or MWCNTs) and imogolite-like nanotubes.^63,64

For graphene (κ = 0) we obtained a bending modulus of 3.42 eV Ų per atom, which is in agreement with the values in the literature.¹⁶ We also found that for thin graphene ribbons, edge effects can significantly decrease the magnitude of the bending modulus. This facilitates the graphene folding. We found that MoS₂ is ∼11 times harder to bend than graphene. We also showed that it is possible to control the bandgap of a nano-sheet by bending it.

In conclusion, we have shown a strategy to compute the bending energy of 2D materials, which can be readily generalized to other 2D systems. It can also be used to predict the mechanical properties of novel heterostructures.¹³ In addition, we calculated how finite size effects, due to the presence of edge atoms, modify electronic and mechanical properties.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the Fondo Nacional de Investigaciones Científicas y Tecnológicas (FONDECYT, Chile) under grants #1150806 (FM), #1150718 (JAV), #1160639 and #1130272 (MK and JR), and Financiamiento Basal para Centros Científicos y Tecnológicos de Excelencia CEDENNA (JR, JAV, FM, FV, RG and MK). FV was supported by the CONICYT Doctoral Fellowship grant #21140948.

References

1. S. Xu, Z. Yan, K.-I. Jang, W. Huang, H. Fu, J. Kim, Z. Wei, M. Flavin, J. McCracken, R. Wang, A. Badea, Y. Liu, D. Xiao, G. Zhou, J. Lee, H. U. Chung, H. Cheng, W. Ren, A. Banks, X. Li, U. Paik, R. G. Nuzzo, Y. Huang, Y. Zhang and J. A. Rogers, Science, 2015, 347, 154–159.
2. C. Ye and V. V. Tsukruk, Science, 2015, 347, 130–131.
3. A. van der Zande and J. Hone, Nat. Photonics, 2015, 6, 804–806.
4. N. Horiuchi, Nat. Photonics, 2015, 6, 804.
5. J. Feng, X. Qian, C.-W. Huang and J. Li, Nat. Photonics, 2012, 6, 866–872.
6. M. K. Blees, A. W. Barnard, P. A. Rose, S. P. Roberts, K. L. McGill, P. Y. Huang, A. R. Ruyack, J. W. Kevek, B. Kobrin, D. A. Muller and P. L. McEuen, Nature, 2015, 524, 204–207.
7. B. F. Grosso and E. Mele, Phys. Rev. Lett., 2015, 115, 195501.
8. R. Breitwieser, Y.-C. Hu, Y. C. Chao, R.-J. Li, Y. R. Tzeng, L.-J. Li, S.-C. Liou, K. C. Lin, C. W. Chen and W. W. Pai, Carbon, 2014, 77, 236–243.
9. D. Shikai and B. Vikas, Mater. Today, 2016, 19, 197–212.
10. W. Bao, F. Miao, Z. Chen, H. Zhang, W. Jang, C. Dames and C. N. Lau, Nat. Nanotechnol., 2009, 4, 562–566.
11. Y. Zhou, Y. Chen, B. Liu, S. Wang, Z. Yang and M. Hu, Carbon, 2015, 84, 263–271.
12. A. Fasolino, J. Los and M. Katsnelson, Nat. Mater., 2007, 6, 858–861.
13. G. Gibney, Nature, 2015, 522, 274.
14. K. N. Kudin, G. E. Scuseria and B. I. Yakobson, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 235406.
15. J. Mintmire and C. White, Carbon, 1995, 33, 893–902.
16. M. Arroyo and T. Belytschko, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 69, 115415.
17. Q. Lu, M. Arroyo and R. Huang, J. Phys. D: Appl. Phys., 2009, 42, 102002.
18. S. Xiong and G. Cao, Nanotechnology, 2015, 26, 185705.
19. Y.-T. Fu, G. D. Zartman, M. Yoonessi, L. F. Drummy and H. Heinz, J. Phys. Chem. C, 2011, 115, 22292–22300.
20. S. Konduri, S. Mukherjee and S. Nair, ACS Nano, 2007, 1, 393.
21. R. I. González, R. Ramirez, J. Rogan, J. A. Valdivia, F. Munoz, F. Valencia, M. Ramirez and M. Kiwi, J. Phys. Chem. C, 2014, 118, 28227–28233.
22. A. Thill, B. Guiose, M. Bacia-Verloop, V. Geertsen and L. Belloni, J. Phys. Chem. C, 2012, 116, 26841.
23. S. U. Lee, C. C. Young, S. G. Youm and D. Sohn, J. Phys. Chem. C, 2011, 115, 5226.
24. L. Guimaraes, A. N. Enyashin, J. Frenzel, T. Heine, H. A. Duarte and G. Seifert, ACS Nano, 2007, 1, 362.
25. P. Maillet, C. Levard, O. Spalla, A. Madison, J. Rose and A. Thill, Phys. Chem. Chem. Phys., 2011, 13, 2682.
26. G. I. Yucelen, R. P. Choudhury, A. Vyalikh, U. Scheler, H. W. Beckham and S. Nair, J. Am. Chem. Soc., 2011, 133, 5397.
27. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
28. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
29. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15.
30. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
31. S. J. Plimpton, J. Comp. Physiol., 1995, 117, 1.
32. E. Bitzek, P. Koskinen, F. Gahler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
33. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
34. K. Chenoweth, A. C. Van Duin and W. A. Goddard, J. Phys. Chem. A, 2008, 112, 1040–1053.
35. S. G. Srinivasan, A. C. van Duin and P. Ganesh, J. Phys. Chem. A, 2015, 119, 571–580.
36. B. D. Jensen, K. E. Wise and G. M. Odegard, J. Phys. Chem. A, 2015, 119, 9710–9721.
37. K. Yoon, A. Ostadhossein and A. C. van Duin, Carbon, 2016, 99, 58–64.
38. B. D. Jensen, K. E. Wise and G. M. Odegard, Modell. Simul. Mater. Sci. Eng., 2016, 24, 025012.
39. T. Liang, S. R. Phillpot and S. B. Sinnott, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 245110.
40. T. Liang, S. R. Phillpot and S. B. Sinnott, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 199903.
41. E. W. Bucholz and S. B. Sinnott, J. Appl. Phys., 2012, 112, 123510.
42. J. A. Stewart and D. E. Spearot, Modell. Simul. Mater. Sci. Eng., 2013, 21, 045003.
43. A. McCreary, R. Ghosh, M. Amani, J. Wang, K.-A. N. Duerloo, A. Sharma, K. Jarvis, E. J. Reed, A. M. Dongare, S. K. Banerjee, M. Terrones, R. R. Namburu and M. Dubey, ACS Nano, 2016, 10, 3186–3197.
44. K. Q. Dang, J. P. Simpson and D. E. Spearot, Scr. Mater., 2014, 76, 41–44.
45. S. K. Singh, M. Neek-Amal, S. Costamagna and F. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 014101.
46. X. Wang, A. Tabarraei and D. E. Spearot, Nanotechnology, 2015, 26, 175703.
47. S. Xiong and G. Cao, Nanotechnology, 2016, 27, 105701.
48. R. T. Cygan, J. J. Liang and A. G. Kalinichev, J. Phys. Chem. B, 2004, 108, 1255.
49. R. I. González, R. Ramírez, J. Rogan, J. A. Valdivia, F. Munoz, F. Valencia, M. Ramírez and M. Kiwi, AIP Conf. Proc., 2015, 1702, 050004.
50. R. I. González, J. Rogan, E. M. Bringa and J. A. Valdivia, J. Phys. Chem. C, 2016, 120, 14428–14434.
51. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
52. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
53. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
54. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
55. F. Munoz, A. Romero, J. Mejia-Lopez, I. V. Roshchin, R. González and M. Kiwi, J. Magn. Magn. Mater., 2015, 393, 226–232.
56. K. Nakada, M. Fujita, G. Dresselhaus and M. S. Dresselhaus, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 17954–17961.
57. J. Xiao, M. Long, X. Li, H. Xu, H. Huang and Y. Gao, Sci. Rep., 2014, 4, 4327.
58. G. Casillas, U. Santiago, H. BarroÌĄn, D. Alducin, A. Ponce and M. JoseÌĄ-YacamaÌĄn, J. Phys. Chem. C, 2014, 119, 710–715.
59. P. D. G. Cradwick, V. C. Farmer, J. D. Russell, C. R. Masson, K. Wada and N. Yoshinaga, Nat. Phys. Sci., 1972, 240, 187.
60. X. Zhuang, Y. Mai, D. Wu, F. Zhang and X. Feng, Adv. Mater., 2015, 27, 403–427.
61. Y. Wang, H. Zhan, C. Yang, Y. Xiang and Y. Zhang, Comput. Mater. Sci., 2015, 96, 300–305.
62. D. Xia, Q. Xue, K. Yan and C. Lv, Appl. Surf. Sci., 2012, 258, 1964–1970.
63. A. Kuc and T. Heine, Adv. Mater., 2009, 21, 4353.
64. M. Ramírez, R. I. González, F. Munoz, J. A. Valdivia, J. Rogan and M. Kiwi, AIP Conf. Proc., 2015, 1702, 050005.

---

Footnote

† Electronic supplementary information (ESI) available: Energy vs. curvature curves for imogolite sheets with N_θ = 4 and 6 are included. Also, an illustration of the bending process is included. See DOI: 10.1039/c7ra10983k

---