Bending energy of 2D materials: graphene, MoS2 and imogolite

The bending process of 2D materials, subject to an external force, is investigated, and applied to graphene, molybdenum disulphide (MoS2), and imogolite. For graphene we obtained 3.43 eV Å2 per atom for the bending modulus, which is in good agreement with the literature. We found that MoS2 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 MoS2, 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.


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][2][3][4][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. 8 Also, supported graphene monolayers have been intensely studied. 9 The effect of different metallic substrates has been investigated as well, showing that roughness can be tuned by varying the support. 10 Recently, Zhou et al. 11 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. 12 Moreover, it could have signicant effects on their electronic properties.
Other 2D materials like h-BN, MoS 2 or WSe 2 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 13 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 specic theory level or code. It uses an external force, specifically a harmonic potential with the desired geometry, that imposes a specic energy minimum. While we use a cylindrical geometry, other congurations 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 2 single layer 18 with a method that applied external forces to the edges of a MoS 2 sheet. Another approach uses a bent carbon template to control the curvature of aluminosilicates, such as a montmorillonite layer. 19 We start with a simple casegraphenewhich allows us to test our method against available data. 14-16 Next, we apply the method to a MoS 2 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 z0.5 nm thick. In contrast to graphene and MoS 2 , the bending energy of imogolite has a minimum as a function of diameter, and is highly monodisperse (i.e. quite uniform in diameter). [20][21][22][23][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. 25 models the growth of imogolite in two stages. During the rst 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. 26 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.

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 2 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 2 .
Our initial graphene conguration is a at 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. Thereaer, we adopt (in small steps) nite 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 innitely 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Å À1 , 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 À1 ¼ k, is given by where E is the total potential energy of the system (ignoring the external energy, E ext ), N is the number of atoms and k is the curvature imposed by the external force (eqn (1), see Fig. 1b). The bending modulus 16 C b is dened as the second derivative of the strain energy with respect to the curvature k

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. 31 The bending due to the external force was combined with energy optimization strategies. We used the fast inertial relaxation engine 32 (FIRE) together with conjugate gradient minimization. We chose this combination because the FIRE is able to avoid getting stuck in local minimum energy congurations, while the conjugate gradient is faster. The graphene C-C atomic interactions were modeled by the ReaxFF force eld. 33 The ReaxFF parameters for the C-C interactions were rst generated by Chenoweth et al. 34 and later improved by Srinivasan et al. 35 This force eld has been tested for different carbon allotropes such as fullerenes, graphene, carbon nanotubes, amorphous carbon and diamond, with very good results. [35][36][37][38] The MoS 2 sheet was modeled with a secondgeneration reactive bond-order formalism by Liang et al., 39,40 and has been used to investigate the mechanical properties of MoS 2 nanotubes, 41 nano-indentation of crystalline MoS 2 , 42 for a few layers of MoS 2 under tension, 43 and for the properties of a MoS 2 layer. 18,[44][45][46][47] Finally, in order to model imogolite atoms we used the CLAYFF potential, 48 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., 48 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 53 was employed. The energy cutoff was set to 400 eV for graphene and to 300 eV for MoS 2 . We used 20 Kpoints 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 54 for the graphics and for the post-processing of our simulations, and PyProcar was employed to perform the band structure analysis. 55 3 Results and discussion

Graphene
The initial conguration 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 k, 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 nally adopts the nanotube geometry (k z 0.78 nm À1 ). 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 at; an intermediate curvature; and nally 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 rst, 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 2 coordinated.
In Fig. 1b we tted a quadratic 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 z95% of the 8 nm value. In other words, ribbons less than $2.5 nm wide are substantially easier to fold.
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Å 2 per atom. This value is closer to the DFT values 14,15 of 3.9 and 4.3 eVÅ 2 per atom than the classically calculated ones of 1.8 and 2.2 eVÅ 2 per atom. 16 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 k $ 2p 3a to p a , 56 due to the breaking of the sub-lattice symmetry at the boundaries. Additionally, there are at bands, typical of the dangling bonds of terminal C atoms (they would disappear aer hydrogenation). These dangling bonds are somewhat hybridized with bulk states for the bent GNR, which is the only noticeable difference compared to a at GNR.

MoS 2
Increasing the complexity of the 2D material under study, we now investigate the bending of a MoS 2 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 2 monolayer with the graphene case, we compute the bending stiffness 47 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 2 ribbon turns out to be more than 11 times 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Å 2 per atom for graphene (1/d / 0 limit).
larger than that of a graphene ribbon of the same width. Additionally, we estimate the Föppl-von-Kármán number 6 of a MoS 2 monolayer to be about 40 times smaller than that of graphene. The Föppl-von-Kármán number is dened as E/D, where E corresponds to the Young modulus. 18,38 Inspection of Fig. 4b of the paper by Blees et al. 6 suggests that displacement, due to a laser, is negligible (less than the measurement noise). This argument does not rule out the feasibility of building nanokirigami 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 2 force elds. 47 For the molecular mechanics they reported a bending stiffness of 8.7-13.4 eV depending on the force eld used. On the other hand, they compared this with previous DFT results from Xiao et al. based on MoS 2 single wall nanotube calculations. 57 They have shown that the force eld 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 2 bulk layer of D ¼ 16 eV. Casillas et al. 58 reported a bending stiffness value of 4.1-13.24 eV, based on continuum theory. Our results for a MoS 2 single layer are slightly over the values reported by other computational calculations and experiments.
In relation to the MoS 2 nite width effects, in contrast with graphene ribbons, the difference between a narrow (0.95 nm) and a wide (8 nm) MoS 2 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 2 ribbon due to curvature we employ rst 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 ribbonas expected of an insulating materialleaving 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 2 band structure is shown in Fig. 5 for a at and nite 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 / N, 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 nding suggests that the bending effects on the electronic structure are rather modest; in fact, a bent MoS 2 layer still behaves like at MoS 2 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 (G 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.

Imogolite
Finally we report our results for the strain energy of a curved imogolite sheet. In Fig. 7 we provide on the le hand side an axial view of the imogolite tube, and on the right a view of its basic unit [(OH) 3 Al 2 O 3 SiOH] 2 . The single walled imogolite nanotube has the chemical composition [(OH) 3 Al 2 O 3 SiOH] 2N q ; this unit is repeated angularly, around the longitudinal axis, N q times, and N Z times along the longitudinal direction. 59 In Fig. 8 we show the energy of an imogolite sheet with this composition and N q ¼ 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 conguration: the initial planar sheet; the minimum energy conguration (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 eld, occurs for N q ¼ 10. In the present case, a ribbon of half the width (N q ¼ 5) reaches its   6 Electronic band gap at the G and M points as a function of curvature. The system band gap is always located at G; the band gap at M is just a parameter to quantify bending effects. Surface states were ignored in the gap measurement.  The force is applied on the Al atoms. A minimum of the bending energy occurs for 9 < N q < 10. The imogolite structure inset corresponds to N q ¼ 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 q value. Similar data, for imogolite sheets for N q ¼ 4 and 6, are included as ESI. † Green Al; red O; yellow Si; and light gray H. 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.,25 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 q < 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, 26 nanotubes 9 # N q # 15 can form, since the energy required is quite small.

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 2 and imogolite. Moreover, our model is straightforward to generalize to other types of nanosheets, 60 like h-BN, MoSe 2 and TiO 2 . It can also be applied to nano-sheets that roll around different kinds of nanocylinder 61,62 such as TiO 2 nanotubes, metallic nanowires, CNTs (or MWCNTs) and imogolite-like nanotubes. 63,64 For graphene (k ¼ 0) we obtained a bending modulus of 3.42 eVÅ 2 per atom, which is in agreement with the values in the literature. 16 We also found that for thin graphene ribbons, edge effects can signicantly decrease the magnitude of the bending modulus. This facilitates the graphene folding. We found that MoS 2 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. 13 In addition, we calculated how nite size effects, due to the presence of edge atoms, modify electronic and mechanical properties.

Conflicts of interest
There are no conicts to declare.