Morphing graphene-based systems for applications: perspectives from simulations

Graphene, the one-atom-thick sp2 hybridized carbon crystal, displays unique electronic, structural and mechanical properties, which promise a large number of interesting applications in diverse high tech fields. Many of these applications require its functionalization, e.g. with substitution of carbon atoms or adhesion of chemical species, creation of defects, modification of structure of morphology, to open electronic band gap to use it in electronics, or to create 3D frameworks for volumetric applications. Understanding the morphology-properties relationship is the first step to efficiently functionalize graphene. Therefore a great theoretical effort has been recently devoted to model graphene in different conditions and with different approaches involving different level of accuracy and resolution. Here we review the modeling approaches to graphene systems, with a special focus on atomistic level methods, but extending our analysis onto coarser scales. We illustrate the methods by means of applications with possible potential impact.


Introduction
Since its first isolation 1 and especially after the assignment of the Nobel Prize for Physics in 2010 for its discovery, graphene has triggered great expectations. Most of its properties depend on the coincidence of an number of physical-chemical circumstances: three of the four valence electrons of carbon organize in three sp 2 hybridized orbitals producing an exactly planar geometry, while the forth contribute to a vast electronic delocalization over the plane, giving great stability to the electronic and geometric structure of the single sheet and very weak inter-layer interaction mediated only by van der Waals (vdW) forces and favoring bidimensionality. In addition, the specific geometry of in-plane bonds generates the highly symmetric honeycomb lattice, responsible for the peculiar band structure with conic geometry at the K points of the Brillouin Zone (BZ). This is responsible for the exotic electronic properties of graphene 2 , such as the ultra-relativistic-like high mobility conduction and the possible existence of Majorana fermions 3 .
Its -equally exceptional -mechanical properties, have received less attention. The peculiar chemistry and symmetry also brings an extremely large resistance to tensile strain especially considering its bidimensionality (Young modulus of the order of TPa 4 , five-fold that of steel), coupled to very low out of plane (bending) rigidity κ (1-2 eV 5 , 6 , 7 ) which is at least one order of magnitude smaller than that expected from materials with comparable in-plane rigidity 8 . This implies an extreme flexibility associated to strength suggesting a vast range of applications in a variety of high-tech fields 9 , included those involving exposure to extreme environmental conditions, such as aero-space technologies 10 . Flexibility and bidimensionality also imply the existence of low frequency quadratically dispersive acoustic phonon branches associate to out-of-plane displacement ("flexural" phonons), which can be described as travelling ripples 11 and have a main role in thermal behavior and related phenomena 12 .
Many of the applications, however, require graphene morphology manipulation of some kind (see Fig 1). Considering for instance nano-electronics, graphene has null density of carriers at the Fermi level, therefore needing doping to be used as conductor, or band gap opening to be used as semiconductor. These properties can be obtained by escaping from the perfect infinite 2D crystal case, e.g. in multilayers or breaking the honeycomb lattice symmetry. Graphene symmetry disruption is not difficult per se, but it is difficult to achieve in a controlled fashion.
Graphene is lightweight and with a huge surface to mass ratio, therefore is considered potentially interesting for gas storage applications. However, it is rather inert, little reactive and with weak vdW interactions. Therefore, this class of applications requires enhancing either physical or chemical interactions. Considering for instance hydrogen storage 13 for energy applications, accurate experimental evaluations of physisorption 14 show that a limiting value for gravimetric (i.e. the relative hydrogen mass stored with respect to the mass system) density in graphene-based systems is ~1% at room temperature (5%-6% at cryogenic temperatures). On the other hand, considering chemisorption one can in principle reach 6-8% gravimetric density at room temperature, but high chemi(de)sorption barriers makes the kinetic of the process very slow at room temperature 13 . Therefore, efforts are directed both to enhance physical-like interactions and to lower the kinetic barriers for chemisorption. Again, manipulation of morphology is required. As for electronics, achievement of interesting properties must proceed through rupture of the perfect graphene symmetry.
While for electronics 2-dimensionality is useful, for other applications the large surface to mass ratio must be declined in a 3D context. This is the case of energy storage applications, including also batteries and supercapacitors, besides hydrogen (or other gases) storage. Therefore great efforts are directed to building 3D graphene based frameworks 15 . This involves including spacers between sheets 16 with tailored size and mechanical properties, and control their amount and distribution 17 . Building 3D networks, controlled chemical functionalization is crucial in a number of other applications 18 , including bio-medical ones 19 . Figure 1. A schematic illustration of the relationship between graphene modification and properties in graphene. The bare single layer graphene is located in the upper-left corner. Background increasing level of shade is roughly proportional to the level of manipulation, from simpler ones (inclusion of defects or substitutions, rippling) to the creation of covalent networks. The effects on electronic and chemical properties and physical interactions are reported in white. The arrows indicate the inter-relationships: e.g. the substrate create rippling and in some cases doping and gap opening; chemical substitutions with N and B induce chemical doping and bad gap opening; defects opens the gap and induce curvature; curvature in turn induce local change in reactivity, and consequently controlled decoration and band gap opening. Finally, locally enhanced reactivity could be used to control decoration, chemical functionalization and finally to build 3D covalent network with enhanced physisorption properties. [The images of the chemical functionalization and 3D networks were adapted from ref [16]].
It appears that manipulation of graphene morphology, structure, electronic properties, chemistry and functionalization are interconnected tasks assuming a key role in all emerging graphene related technologies (see Fig 1). In this this context, theoretical investigations and computer modeling are tools of outmost importance to study and possibly design properties of these new materials. The socalled multi-scale approach, traditionally designed and used for biomolecular systems 20 , where the hierarchical organization in different length scales is more apparent, is currently being extended to materials 21 . In the cases under discussion, multi-scaling mainly involves the two basic levels of representations, namely the "ab initio" quantum chemistry (QM) level, with explicit electrons capable of describing the chemistry of the system, and the classical "molecular mechanics" (MM) representing the inter-atomic interactions with empirical force fields (FF).
Super-atomic descriptions are also "natural" in specific cases: For instance, fullerenes were sometimes treated with a Coarse Grained sphere-of-beads model 22 similar to that used for biomolecules 23 . For the graphene sheets, conversely, a natural low-resolution representation is a 2D membrane endowed with appropriate mechanical properties 24 , analogously to those used for biological membranes.
The scope of this report is to analyze the contribution of modeling to the morphing of graphene-based system for applications. Therefore, for each of the different types of morphing described, the focus is on theory and models able to better describe and predict the properties under examination. The well-known properties of perfect single sheet graphene are also summarized in the next section to have a base for comparison for changes due to manipulation. The subsequent section describes several ways aimed at controlling the electronic properties in graphene-derived systems. The manipulation of both physical and chemical interactions of graphene is addressed in section 4. Section 5 focuses on dynamical deformation of graphene. Section 6 reports a summary and conclusions.

Single sheet suspended graphene
The electronic and structural properties of graphene are very well known and described in a number of very good reviews 2 . The minimal unit cell of graphene (represented in blue in Figure 2a) includes two C atoms. The eight valence electrons organize in four filled bands whose structure was calculated with all available electronic structure methods 25 . The band structure evaluated within the Density Functional Theory (DFT) framework (PBE functional 26 ) is reported in Fig 2c (left plot). The π band crosses its empty counterpart (in red, the filled bands are in blue) at the K point, where the linear dispersion is described in terms of the Fermi velocity v F . Inclusion of many body effects show on average relatively small corrections to the DFT picture, resulting in a renormalization of v F , whose entity is still under debate 27 .
Clearly, when treating an isolated single sheet graphene with no defects or structural deformation, the unit cell is usually the preferred choice for the model system. However, there are infinite ways of defining a graphene supercell. Larger, rotated, or differently shaped cells might be necessary in the presence of substrates or other elements breaking the symmetry of the honeycomb lattice. A selection of them is reported in Figure 2a, while their corresponding Brillouin Zones (BZ) in Figure 2b. The size, shape and orientation of the BZs change accordingly. The definition of different BZs implies that the description of the band structure is celldependent. As an example, in Figure 2c we report the comparison of the graphene band structure evaluated in DFT-PBE in unit cell and in the 4√3×4√3R30 cell 26,28 (colored in green in Figure 2a,b), including 96 C atoms. Once the scale of energies of the two band structures are aligned, the bands of the 4√3×4√3R30 cell appear refolded. In order to understand the refolding, we follow for instance the band along the M-K line of the unit cell, which is fragmented in a sequence of Γ-M-Γ-M-Γ lines of the 4√3×4√3R30 cell, as it can be seen in Figure 2b. The fragmentation is reported in the central plot, and the refolded empty and filled π bands are highlighted in the right band structure. Interestingly, the crossing of bands at K point is remapped onto the Γ point in the 4√3×4√3R30 cell, as an effect of refolding/rotation. It is also interesting to remark that, while the band representation changes, the value of physical observables should not. As a matter of fact the density of electronic states (DoS) evaluated in the two cases coincide (within numerical error) once the energy scales are aligned. This is true in particular for the two peaks at ±1.5eV appearing where the π bands become flat (π edges).
Finally, we remark that the use of rectangular cells allows remapping the main symmetry directions along the two Cartesian directions x and y.

Modification of graphene electronic properties
Any breaking of the graphene perfect symmetry induces changes in the electronic or transport properties. In the following we review the effects of substrates, chemical substitution/adhesion and structural modifications on the electronic properties.

B and N substitutions and structural defects
Hexagonal Boron Nitride, BN, has a 2D honeycomb structure identical with that of graphene, with the two triangular sublattices occupied by B and N, respectively. The sublattices symmetry breaking induce the opening of a large band gap, estimated in ~4.5 eV within DFT scheme (generalized gradient approximation 29 ), making it the insulating material more similar to graphene. Graphene partially substituted with BN patches or strips of different configuration shows intermediate behavior and can be considered therefore semiconductor systems with band gap tunable as a function of the amount and size of BN areas with respect to graphene. Independent DFT estimates 30,29,26 indicate that the band gap opening is roughly linear at low B/N density, and can be estimated in ~0.03-0.06eV for 1% of substituted C atoms. In presence of excess of N or B, doping (of n-type or p-type respectively) is also present, estimated in ±0.3-0.4eV of Fermi level shift with respect to Dirac point for each % point of N(B) excess 31 . In this case, for given shapes of the substitution patches, mid-gap states and specific magnetic properties can also appear.
Purely structural defects in graphene such as isolated or aligned dislocations delimiting grain boundaries can also produce severe changes in the electronic structure 32 . Depending on the specific periodicity and topology, and on the relative orientation of domains, the system can display conductive behavior or band gaps with reflective transport behavior at the grain boundaries 33 . Vacancy type and Stone-Wales type defects can open band gaps on the scale of eV leading to semiconducting behavior 34 which opens to the possibility of using defects to engineer the electronic properties of graphene. Clearly, due to rupture of the regular π electronic system, structural defects of any kind also produce structural deformation and changes in local reactivity, which are discussed in Section 4.

Effects of the substrate
In supported graphene, the interaction with the substrate can induce electronic properties modifications directly, or indirectly (e.g. by substrate-induced defects) or influence electronic modification due to other effects. Considering graphene grown on SiC by Si evaporation, for instance, it was early inferred by DFT studies that -for grown on the C-rich face, the substrate can enhance the effect of Al, P, N and B substitutional doping 35 . Of all dopants, N is shown to prefer the substituting the graphene upper layer, while other dopants prefer interstitial or buffer layer locations. The enhancement effect seems due to an increased stability of dopants operated by interaction with the substrate. In a model without the buffer layer 36 (whose presence is still debated for grown on the C-exposed surface) a direct effect of the substrate producing n-type doping of ~0.3-0.4 eV, was reported, very sensitive to the intercalation of passivant H atoms. This work also reports a marked sensitivity to the use of different density functionals and treatments of vdW interactions, confirming a main role of the interaction with the substrate. In the Si exposed surface the presence of a covalently bond buffer layer was early established, and calculations with a minimalist cell shows noticeable levels of doping 37 .
Early calculations however, were performed in very small simulation cells, inducing strain on the graphene layers to achieve commensuration with the substrate. Subsequent calculations in larger and more relaxed supercells revealed a reduction of all the effects: on the C-rich surface, the first graphene layer becomes almost detached and recover graphene-type metallicity, with little doping 38 . In the case of the Si-rich surface the stress can be completely relaxed only using a very large supercell 39 , where doping seems to be negligible. The question is still controversial however, because experimental observations indicate a variable level of n-doping in monolayer graphene, which tends to decrease in quasi free standing graphene (QFSG) obtained intercalating H underneath the buffer layer, depending on the amount of H, and negligible interactions of QFSG with the substrate are confirmed by DFT calculation in the fully relaxed model 40 , reporting a ~5-8meV band gap and basically no doping 41 . A more advance theoretical analysis indicate that perturbation from the non-interacting situation and others effects (such as the appearance of specific non linearities in the bands) might be due to many body electron interaction and phonon-coupling effects not included in the standard DFT treatment 42 . Conversely, metallic substrates induce a genuine doping whose origin is ascribed by the difference in work function and to the more or less "physical" or "chemical" interactions between graphene-metal. A systematic DFT based study shows a n-type doping for Al, Ag, Cu, and p-type for Au and Pd 43 .
Effects of substrates and of defects can combine. For instance, vacancy type defects created on QFSG are shown to strongly interact with the H coverage and induce doping and/or localized states 44 and metals intercalation (e.g. Li) produces a strong n-type doping. In addition, the procedures to induce substitutional doping, e.g. with N, often induce also vacancy type nitrogenated defects 45 (e.g. with pyridinic or pyrrolic type reconstruction), which combines and often enhances the doping and band gap opening effects, allowing in principle a fine tuning of the electronic properties of the graphene sheet 46 .

Effects of adatoms
Chemisorption of adatoms (typically H and F) is another common way to change graphene properties. Because the fully hydrogenated sheet -graphane -is an insulator 47 , partially hydrogenated graphene display semiconducting properties, with band gap depending on the coverage and decoration. Decoration in strips produces the typical ribbon-like behavior 48 , with gap dependent on the edge type (zig-zag or harmchair) and decreasing with the width of the hydrogenated strips, and therefore increasing with the H coverage. Decoration in patches produces very variable situations 49,50 depending on the shape and connectivity of the patches, sometimes including midgap non dispersive states 51 . However, when experimental and DFT data are put together as a function of the H coverage, irrespective of the decoration type ( Figure 3), they seem to accumulate onto a line which was empirically fitted with an almost square root behavior (precisely, a 0.6 exponential 51 ). Best adherence to the empirical curve is observed for the uniform coverage 52 or in general in experimental cases, when the coverage is likely to be more random. Cases of regular or symmetric coverage seem to bring larger dispersion from the curve. Considering the dispersion of the curve, a measurement of the band gap e.g. with STS techniques could give an evaluation of H coverage with an error of 10-15%. Due to its symmetry, single layer graphene is not piezoelectric, although there are indications of its reactivity to electric fields by means of flexo-electric behavior 53,54 . Again, the sheets functionalization with specific elements (H, F, Li, K) was shown an effective way to induce piezo-electricity: fluorination and decoration with Lithium, in particular, are shown to bring piezoelectric coefficients of the order of magnitude of those of 3D matherials 55 . Piezo-electricity can also be induced creating holes with the right symmetry on the layer 56 , again breaking the graphene inversion symmetry. Chemical manipulation and substitution also influences the flexo-electric properties: curvature changes induced by an external electrostatic field are enhanced in the presence of N substitutions 26 (Figure 4).  Electronic and transport properties dependence on the topology of defects and relative orientation of grain boundaries [32]   The QFSG is weakly interacting, with very small band gap and basically no doping [40] (2011) Quasi free standing graphene on SiC 13×13 with buffer layer and substrate and complete intercaled H coverage; PBE+vdW The interaction between substrate and graphene is very weak, leading uniform electronic density on rather flat sheet [42]  In conclusion to this section, we remark that while standard DFT methods seems appropriate to evaluate the energetics and electronic properties for chemical manipulations, the effect of substrates, especially when the interaction is more physical in nature, calls into play the use of more advanced DFT schemes, involving the use of vdW corrections and evaluation of many body electronic effects.

Reactivity and interactions manipulation
From previous discussion the necessity of controlling adsorption or substitution of chemical species onto graphene emerges as a way to endow graphene with electronic and electro-mechanical properties of outmost interest in applications. Controlling interaction of graphene with atomic or molecular species has also a direct use in energy applications: graphene has a huge surface to mass ratio, therefore it would be in principle an ideal candidate not only for gas (mainly hydrogen) storage applications, but also for super-capacitors and batteries 57 . In addition, reactivity control is a key step for building 3D networks with organic inter-layer spacers. In the following a selection of theoretical and simulation works on adsorption of gases (mainly hydrogen) in graphene based system is reported. A summary of the literature is in Table 2.

Physisorption
Taking hydrogen as a paradigmatic case for molecular interactions with graphene, one can generally separate adsorption onto graphene in two main classes, chemisorption and physisorption 13 . Physisorption occurs barrierless via dispersive vdW interactions in case of neutral molecules. In the case of H 2 , the binding energy on bare graphene is estimated around 0.01eV 58 , bringing a precise linear dependence of the gravimetric density onto the specific surface area 14 and posing a strict upper limit to the gravimetric density for physisorption onto graphene at 1-2% at room temperature. Besides the obvious strategy of working at cryogenic temperature, in several theoretical works physical interactions were shown enhanced within nanocavities produced by rippling 58 , within multilayers with specific (nano-metric) interlayer spacing 59,60 or with nano-sized perforations of the sheet 61 .
From the methodological point of view, it is interesting to observe that these works address the problem of physisorption combining different methodologies besides DFT (see Table 2), including classical approaches based on empirical Force Fields, different types of dynamical treatment of hydrogen (QM dynamics, Monte Carlo and thermodynamic evaluation of the adsorption). This is basically due to two circumstances: first, gravimetric density evaluation must be treated at the statistical level, which requires large system and sampling methods proper to this aim. Second the dispersive interactions are very elusive and difficult to treat with standard DFT methods. Commonly used LDA and GGA Density Functionals are based on the (semi)local representation of exchange-correlation, insufficient to address instantaneous fluctuations of the electron density. As a consequence, these functionals fail in reproducing the medium-long range behavior of inter-atomic interactions 62,63 , which was corrected with a number of semiempirical schemes. In addition, hydrogen, the lightest element, calls a quantum treatment of its dynamics. Because all of these circumstances make the simulation unaffordable on large systems, recursion to combination of high accuracy/high cost with lower accuracy/low cost methods in multi-scale combined approach is compulsory. Overall, however, the lesson one learns from this analysis is that cavities, convexities, porosity or in general structures enhance the physisorption at best when they are at the nano-scale (specifically ~1nm). Therefore nano-scale structuring must be considered a guideline also in designing graphene-based 3D networks.

Manipulation of reactivity
In contrast to physisorption, chemisorption on graphene is barrier driven process, requiring a change in hybridization from sp 2 to sp 3 64, 65 . The barrier for atomic H is estimated around 0.3eV, but the barrier for H 2 , involving dissociative chemisorption, is ~1.5eV per atom. Once adsorbed H can hop from site to site with a barrier of ~0.75eV 66 . These processes display slow kinetic at room temperature.
Enhancement of reactivity can be obtained by perturbing the π delocalization in several ways, which mostly superimpose to those already mentioned to manipulate electronic structure. Structural and substitutional 67,68,69,70 defects generally constitutes hot-spots of reactivity. The action mechanism is usually ascribed to the presence of localized states or doping due to the defect, which are able to induce dissociation of the molecule and subsequent adsorption of the reactive radicals, resulting in a reduction of the dissociative chemisorption barrier. External electric fields were shown to act as catalysts for this process 71,72 , being able to deform the already perturbed orbitals. This effect is synergetic with that of substitutional doping 73,74,75 promising fast kinetics, large GD (6.73%) at room temperature and the possibility of using electric field for switch uptake/release of molecules.
The reactivity enhancement effect of purely structural defects is partially due to the local distortion from flatness, which induces a protruding and partial pyramidalization of a C site, favoring the sp 2 to sp 3 transition. In fact, significant lowering of the barriers for H 2 dissociation were previously observed on nanotubes 76 , while stabilization of the adsorbate was observed on corrugated graphene without defects 77,78, 58 . This indicate that control of local curvature can induce control of graphene decoration, upload and even release, as demonstrated by a simulation showing H 2 release by curvature inversion 58 . A side (but related) effect, is that the band gap also depends on mechanical deformation, therefore rippling 28 or stretching 79 could be used directly to manipulate electronic properties, and indirectly by curvature induced chemical functionalization. Because chemisorption involves much larger energies than physisorption, in a first approximation vdW correction to standard DFT schemes play a less important role. However they become crucial in the approaching phase of the molecule, when different orientations and locations with respect to C atoms determine the effective reaction path and barrier for the chemisorption reaction 80,81,82 .

Metal mediated adsorption
Metal mediated adsorption of hydrogen is often located 83 , in between chemisorption and physisorption. This interaction is specifically mediated by d orbitals in transition (e.g. Ti 84 ) or heavy metals (e.g. Pd 85 , 86 which also works as catalyst and dissociate the molecule producing a "spillover" effect of atomic H onto graphene), while it is often described by an "enhanced" vdW interaction for lighter metals such as Li and Ca 87 . In both cases the interaction energy is increased of at least one order of magnitude, which suggests useful applications in storage. However the problem is turned in how to control the distribution of metal onto graphene. In fact in most cases metals tend to form clusters, reducing the active surface for adsorption with respect to total mass. Again, combination of different kinds of graphene manipulation has a key role. For instance, it was shown that inducing N-defects onto graphene before exposure to Ti can reduce the cluster size, optimizing the active surface for adsorption 88 . Combination of Li decoration, N-doping and electric field predict reversibility of hydrogenation in a DFT study 89  Multilayer graphene-oxyde framework pillared with diboronic-derived molecules exposed to CO2 ~2nm×2nm×2nm DFT-PBE, DFTB Classical MMFF94 Using realistic framework models with the correct spacing and corrugation enhances selectivity and the adsorption coefficient.
[60] (2015) Multilayer perforated graphene exposed to H2  Morphology manipulations have diverse effect on selectivity of adsorption process of molecules on graphene, so that it can be used for the gas sensor applications. Doping graphene with heteroatoms can results in different amount of valence electrons, and usually protrude (i.e. dopants Al, Si, P, S, Mn, Cr) 90,91 from the graphene surface and create reactive binding sites for molecular adsorption. In general, as in the case of physisorption and more than for chemisorption, the vdW interaction have a main role, therefore the use of appropriate functional for their representation is crucial 92 .

Dynamical morphing
Reversibility of graphene manipulation is a key property in tasks such as hydrogen storage. Mechanical deformations are the most reversible. In fact the reversible storage using the reactivity dependence on local curvature was demonstrated in simulations 58 . Therefore, the problem is turned on how to control local curvature of the graphene sheet. One possibility is to exploit the natural curvature of supported graphene: on given substrates, such as SiC, graphene shows multi-stable patterns of curvature 93 , allowing the possibility of switching between each other by changing external environmental conditions (e.g. temperature, electric fields).
However, graphene offers the unique possibility of changing the local curvature dynamically, exploiting the naturally occurring out of plane wave-like deformations, i.e. the flexural phonons (or acoustic z polarized branches, i.e. ZA phonons). The capability of those vibrations to chemisorbed hydrogen and transport and pump it though a multilayer system was demonstrated by DFT 58 and classical MD simulations based on empirical FFs 11 . Clearly, for such applications, the key issue is how to create and maintain coherent ZA phonons of specific wavelength and amplitude.
Many efforts were devoted to the calculation the phonon dispersion curves of graphene. The reason is that due to the 2d nature of this material, the vibrational properties play a crucial role on its thermal and mechanical properties. Much attention has been reserved to the ZA mode, particularly important in 2d materials as graphene 94 . Its dispersion law is parabolic, implying that these modes are the lowest energies of the spectrum, and are the easiest to be excited. The very low phonon-phonon scattering rate and the large thermal population have as a consequence that the flexural phonons give a fundamental contribution to thermal conductivity both for graphene monolayer 95 and multilayers 96 . Moreover, they are responsible of the negative thermal expansion coefficient observed in a wide temperature interval, up to an inversion temperature value, which is still controversial.
Graphene monolayer and graphite Standard lattice dynamics with LCBOPII empiric interatomic potential Phonon dispersion and bending rigidity [104] Graphene monolayer Tersoff and Brenner potential Optimization of and modification of FF for the accurate reproduction of phonons [105] (2010) Graphene monolayer Tersoff potential, parameters in [11] Calculation of thermal conductivity [106] (2012) Graphene monolayer Monte Carlo with empiric interatomic potentials.
isothermalisobaric ensemble, periodic boundary conditions along x,y, reparameterized long-range carbon bondorder potential LBOP and LCBOP Thermodinamic properties; empirical potentials limited to nearest-neighbour interactions give rather dispersed results.

Graphene monolayer
Brenner (first generation) Analytic form for the bending modulus [108] (2009) Graphene monolayer and singlewall nanotube Analytic from interatomic potential tension and bending rigidity directly from the interatomic potential. [109] (2010) Graphene monolayer Continuum approach; Molecular mechanics, REBO interatomic empirical potential Nonlinear elastic properties under uniaxial stretch and tension [110] (2012) Graphene layer on Si terminated SIC environment-dependent interatomic empirical potential !(EDIP) Graphene layer basic cell of 338 atoms, Determination of the atomistic structure of the graphene buffer layer on Si-terminated SiC. The solution of minimal energy forms a hexagonal pattern composed of stuck regions separated by unbonded rods that release the misfit with the SiC surface [111] Graphene layer Molecular dynamics simulation, Tersoff empirical interatomic potentials for C-C, Si-Si, and Si-C interactions Two layers of carbon atoms of the 6H-SiC!(0001)! subsurface after sublimation of Si atoms undergo a transformation from a diamondlike phase to a graphenelike structure at annealing temperature above 1500 K [112] (2008) Table 3. Vibrational mechanical and thermal properties of graphene systems, theoretical approaches Among the most used computational approaches we can mention the elastic continuum model 97 , 98 , the first principles DFT also with Perdew-Burke-Ernzerhof generalized gradient approximation 99 , 100 , 101 , and the molecular dynamics simulations associated with the use of empirical interatomic potentials for the interactions between carbon atoms, which allows to treat larger systems 102 . The most effective have been demonstrated the Long-range Carbon Bond Order Potential (LCBOPII) 103 and the parameterization by Lindsay et al. 104 of the Tersoff potential 105 , which has been also used to calculate the thermal conductivity of graphene 106 , 107 . Empirical interatomic potentials have been widely used also for the investigation of the mechanical and elastic properties of graphene sheets, in particular for the elastic bending modulus 108 , the thickness and Young's modulus 109 , 110 , also on Si-terminated SiC 111 and 6H-SiC(0001) 112 .

Conclusions and Perspectives
Morphology manipulation of graphene requires flexible experimental and modeling approaches. In particular on the modeling side, both accuracy in the representation of the phenomena, and the use of model systems very large compared to the atomic scale are required. The former calls for the use of at least DFT based schemes, most often corrected to account for long range electron correlation effects to better account for specific electronic properties or for the dispersive part of vdW interaction; in specific cases the recursion to quantum mechanical treatment of nuclei is also required. On the other hand, large model systems are needed to analyze systems at the nano-micro scale with defects and functionalization, which have a low level of symmetry. The evaluation of thermodynamic properties is also often needed, requiring extensively long simulations. This is often incompatible with the high computational cost of accurate methods, therefore the recursion to empirical treatment of interactions is used. In some cases the two methods, ab initio and empirical, are mixed at several levels, intrinsically (such as in the empirical vdW corrections to DFT) or explicitly.
Ab initio level currently allows addressing the 10 5 atoms scale (corresponding approximately to the 100nm size in graphene) in short simulations. Considering the exponential increase of computer power ensured by Moore's law, the μm scale can be reached in less than 10 years, while it is already feasible, in principle, using empirical approaches on computing systems with extensive parallelization.
However, as the scales feasible in simulations increase, new questions emerge on the reliability of the theories underlying the simulation. This is especially true for more empirical approaches, because the parameterization of interactions are usually tested on smaller/shorter size/time scales, but also for more ab initio approaches such as DFT, which always include some hidden level of empiricism. Other sources of "systematic" errors in modeling might rise in the algorithms used for sampling of the conformations and phase space of the system and in the large scale average or thermodynamic properties.
In conclusion, as the size of the model system increases, the comparison with measurements becomes more and more important. The sizes of experiment and of simulations tend to meet -the first coming from bottom, the second from top -at the meso-scale. This allows a direct comparison of measured observable with the corresponding evaluated quantity. Clearly this brings an advantage in the interpretation of experiments, but also a feedback on the theory, allowing creating models more adherent to reality.