Van der Waals interlayer potential of graphitic structures: from Lennard-Jones to Kolmogorov-Crespy and Lebedeva models
aa r X i v : . [ c ond - m a t . m e s - h a ll ] M a r Van der Waals interlayer potential of graphitic structures: from Lennard-Jones toKolmogorov-Crespy and Lebedeva models.
Zbigniew Kozio l ∗ , Grzegorz Gawlik, and Jacek Jagielski
1, 2 National Center for Nuclear Research, Materials Research Laboratory,ul. Andrzeja So ltana 7, 05-400 Otwock-´Swierk, Poland Institute of Electronic Materials Technology, 133 W´olczy´nska Str.,01-919 Warszawa, Poland
The experimental knowledge on interlayer potential of graphenites is summarized and comparedwith computational results based on phenomenological models. Besides Lennard-Jones approxima-tion, the Mie potential is discussed, Kolmogorov-Crespy model and equation of Lebedeva et al. Anagreement is found between a set of reported physical properties of graphite (compressibility alongc-axis under broad pressure range, Raman frequencies for bulk shear and breathing modes underpressure, layer binding energies), when a proper choice of model parameters is made. It is arguedthat the Kolmogorov-Crespy potential is the preferable one for modelling. A simple method of fastnumerical modelling, convenient for accurate estimation of all these discussed physical properties isproposed. It is useful in studies of other van der Waals homo/heterostructures.
PACS numbers: 61.46.-w, 73.20.At, 73.22.-f 61.48.De, 68.35.Gy,73.22.Pr
INTRODUCTION.
There is recently a growing interest in artificial het-erostructures formed by 2-dimensional layers of grapheneand other newly discovered materials, like these of tran-sition metal dichalcogenides (TMDCs), and named ’vander Waals heterostructures’ [1]. The reason of that isreachness of physical phenomena observed and high po-tential of their possible applications [2]. Moir´e patternsresult when two layers are rotated [3], [4] and are relatedto van Hove singularities. Recently, in bilayer graphenetwisted at 1.1 degree superconductivity was found withcritical temperature of 1.7K [5]. By manipulating dop-ing levels these singularities are observed at angles upto 31 degrees [6]. It was shown that free-standing vdWheterostructure of graphene on hexagonal boron nitride(hBN), where a small lattice mismatch exists shows abuckled atomic structure formed as Moir´e pattern [7].There is already a number of laboratory prototypes ofoptoelectronic devices [8], like bolometers [9], photode-tectos [10], [11]. Their characteristics can be tuned byelemental doping, surface chemical doping, intercalationand electrostatic gating.It is believed that artificial van der Waals heterostruc-tures offer a broad field for research on novel materials forgraphene/silicene, graphene/
M oS , and silicene/ M oS systems [2], as well silicene [12] or M oS alone [13]. Ab-initio calculations suggest the existence of magnetism in ReS doped with Co, Fe, and Ni [14]. Formation ofquantum-well Type-I heterostructure in van der Waals- interacting monolayers of M oS and ReS has beendemonstrated [15] and Type-II band alignment was foundin vdW heterojunction diodes based on InSe and
GaS ∗ e-mail: [email protected] and graphene [16]. W Se and M oS monolayers havebeen shown useful to create charge density modulationin electronic band valleys, resulting in valley polariza-tion, and N a -doped
W S under strain is considered acandidate for spintronics [17]. After spin injection from aferromagnetic electrode transport of spin-polarized holeswithin the W Se layer was observed [18]. hBN andgraphene vdW heterojunctions indicate on large poten-tial of their use in spintronics [19]. Light-induced neg-ative differential transconductance phenomenon was re-alized on graphene/ W Se heterojunction transistor [20].DFT calculations show that vdW interactions dominatebetween antimonene and graphene layers and by apply-ing electric field between them one can tune the heightand type of Schottky contacts [21]. Lattice dynamicaltheory and ab initio calculations indicate on the exis-tence of piezoelectricity in 2D lattices comprising h-BN,2 H − M oS , and other transition-metal dichalcogenides[22].It is generally thought that van der Waals interaction,ubiquitous in nature, is caused by temporal fluctuationsof electronic charge that induces dipoles of random ampli-tude [23]. It plays a role in friction and adhesion betweenmaterials, absorption of molecules on surfaces. These rel-atively weak forces are responsive for binding betweenseparate sheets of graphene.While many electronic properties of graphene-basedsystems are known [24] still there are open questionsconcerning the detailed form and physical mechanismsleading to inter-plane potential in these systems. Thereare two types of approaches used to describe the van derWaals potentials of graphitic systems: an ab-initio one,based on density functional theory, and the other usesempirical potentials [25].The DFT calculations often give an underestimatedvalues of binding energy energy in case of van der Waalsinteractions [26], [27], [28], [29], [30] and energy differ-nce between bindings of (preferred) AB and AA typeof stacking, E AB and E AA , which seems too large, aswe will discuss. Recent diffusion quantum Monte Carlocalculations predict for instance E AA /E AB of about 0.65[31]. Results however depend strongly on the functionalsused to model vdW forces [32], [33].In this work we compare how well several properties ofgraphene/graphite can be described by using a few em-pirical potentials. Some of them are well based on DFTcalculations. An advantage of using empirical potentialsrelies on ease of their implementation in any numericalmethods. Moreover, we describe a scheme of comput-ing several basic properties based on a simple analyti-cal approach. The method may be easily extended toinvestigation of vdW homo- and heterostructures otherthan graphene/graphite. We demonstrate that by carefulchoosing of parameters of phenomenological potentialswe may reproduce self-consistently experimental valuesof compressibility of graphite as well Raman shifts ob-served in shear mode under pressure and predict changeof Raman shift under pressure for layer breathing mode(LBM). Lennard-Jones and Mie Potentials.
An often used approximation of vdW potential is the12-6 Lennard-Jones one, U ( r ) = 4 ǫ (cid:0) ( σ/r ) − ( σ/r ) (cid:1) ,with two only adjustable parameters, where r is distancebetween atoms.It reproduces well interlayer distance and elastic con-stant of graphite. It was proved useful in describingseveral basic properties of C molecules [34], with ǫ =2 . meV , σ = 3 . A . It was used for modeling carbonnanotubes [35], with ǫ = 4 . meV and σ = 3 . A . Heet al. [36] obtain an analytical approximation for vdW forces acting between nanotubes. Authors extend theircontinuum model [37] to study vibrations of multilay-ered graphene sheets. It is a convenient approximationin modeling reinforcement of composite materials withcarbon nanotubes [38].A more general form of vdW potential is oneintroduced by Mie [39], [40] and it often repro-duces well properties of carbon compounds, U ( r ) = αǫ (( σ/r ) m − ( σ/r ) n ), where α = ( m/ ( m − n )) · ( m/n ) n/ ( m − n ) . When m = 12 and n = 6, the form ofLennard-Jones potential, with α = 4 is assumed. Kolmogorov-Crespi and Lebedeva Potentials.
There are at least two deficiencies of LJ or Mie type ofpotentials. First is that they are isotropic, i.e. the poten-tial depends on distance between atoms, only. However,in graphenites, the binding between atoms on differentplanes is due to overlap between π electrons from adja-cent layers and as such it ought to depend on the anglebetween orbitals. The second problem, as it was noticedfirst by Kolmogorov and Crespi [41], [42] and will also beshown in this work, is that these potentials produce toosmall energy difference between bindings of (preferred)AB and AA type of stacking, E AB and E AA . Albeitthere is some controversy about how large that energydifference is. The KC potential has an r − two-body vander Waals-like attraction and an exponentially decayingrepulsion terms, very short ranged, falling essentially tozero at two transverse interatomic distances . The direc-tionality of the overlap is reflected by a function whichrapidly decays with the transverse distance ρ (Fig. 1).Most often it is used in the following form, for interactionbetween atoms m and l : U lm = − A (cid:18) z r lm (cid:19) + exp ( − λ ( r lm − z )) · [ C + f ( ρ lm ) + f ( ρ ml )] (1) ρ lm = r lm − ( n l r lm ) (2) ρ ml = r lm − ( n m r lm ) (3) f ( ρ lm ) = exp (cid:0) − ( ρ lm /δ ) (cid:1) · X n =0 C n ( ρ lm /δ ) n (4)where n k is the vector normal to the sp plane in thevicinity of the atom k , and z is close to the interlayerdistance at equilibrium. The summation over n in Eq. 4is usually limited from n = 0 to n = 2 [43].KC potential has been used broadly in numerical mod- eling (molecular dynamics) [3], [4], [44],[45] as well in de-scription of ballistic nanofriction [46], [47].Lebedeva et al. [48], [49] use another form ofanisotropic potential which was obtained by fitting DFTcalculations and was found to describe well the ex-2 mn m n l r lm r ml r lm FIG. 1: Schematics of Kolmogorov-Crespi interaction be-tween atoms l and m located on different flakes of graphene.The vector n k (k=l,m) is normal to the sp plane in the vicin-ity of atom k. perimental graphite compressibility and the corrugationagainst sliding, while Jiang [50] uses 5 for modeling ther-mal conductivity of FLG: U = − A (cid:16) z r (cid:17) + Bexp ( − α ( r − z )) ++ C (cid:0) D ρ + D ρ (cid:1) exp (cid:0) − λ ρ (cid:1) exp (cid:0) − λ (cid:0) z − z (cid:1)(cid:1) (5)where r is the interatomic distance and ρ = √ r − z isthe projection of the distance within the graphene plane. ANALYTICAL APPROXIMATION.
Lattice spacing in graphene in hexagonal c-direction,is known well from X-ray diffraction, it is 3.3538 ˚A atroom temperature, and it does not change strongly withtemperature [51]. The equilibrium interlayers spacing is3.34 ˚A in AB stacked [51] and 3.44 ˚A in turbostraticgraphene [52]. Neutron diffraction studies show that allin-plane C-C lengths are 1 . A ± . A ([53]). Distances between atoms in AB and AA stacking.
There are two types of ordering of atoms on the nearestplane with respect to atoms on another plane. Let us call them type I and type II ordering. The type I results whenatom is placed over another atom of the hexagon and typeII when atom is placed over the center of hexagon. TypeI ordering is the only one in case of AA stacking, and incase of AB stacking equal number of atoms is in orderingof type I and II.We can create a convenient scheme of calculating nu-merically potential energy for different stackings by com-puting it for each of these types of orderings. For that,we need distances projection (in plane) between an atomposition over the plane together with number of atoms atthese distances (rings of equi-distant atoms, as these inFig. 2). These numbers can be found quickly with cus-tom written scripts that compute distance between atomsin a graphene structure which was build also by scripts.The results are provided in Tables I and II, for both typesof ordering. Analytical expressions on these distances areavailable however we do not know a method to expressby formula the number of atoms in equi-distant rings.In numerical simulations, e.g. performed in LAMMPS[43], it is a standard procedure to introduce a cut-offdistance in Equations ?? – 5. Potential energy betweenparticles exceeding that distance is assumed to be equalzero, while the function given by these equations is ap-propriately smoothed out at that distance in order toavoid possible discontinuities in computational results.That cut-off is taken usually at around 12 ˚A in case ofgraphene modelling.Limiting any approximation to atoms listed in TablesI and II is equivalent of using cut-off distance in Eq. ?? of about 15.4 ˚A and then we expect to obtain the sameaccuracy of results as by using LAMMPS. FIG. 2: Rings of equi-distant atoms on neighbouring planesof graphene. On the left is type I, present for AA and ABstacking and type II (right) present in AB stacking, only. ABLE I: Distance between an atom on one graphene plane and nearest neighbors on next plane, for type I of neighbors, asin Fig. 2 (left), found in AA and AB stacking of layers. N is number of neighbors and d is distance (in units of graphene unitcell value of 1.42 ˚A). d N d N d N N is number of neighbors and d is distance (in units of graphene unit cellvalue of 1.42 ˚A). d N d N
12 6 12 12 12 12 24 12 6 12 12 12
Potential energy Φ I ( z ) for an atom in position AA ata distance z from the plane will be given by an infinitesum of contributions from equi-distant atoms of type I :Φ I ( z ) = ∞ X i =1 N ( i ) · U (cid:16)p ( z + ( a · d i ) ) (cid:17) , (6)where the number of neighbors N ( i ) at distance d ( i ) istaken from entries of Table I and U is any of potentialsdiscussed (LJ, Mie, KC or Lebedeva’s). In Eq. 6, a is theunit cell length for in-plane atoms ( a = 1 .
42 ˚A in case ofgraphene).In case of AB planes, there is equal number of atomsthat are in ordering of type I and type II . Equation onenergy for type II atoms, Φ II ( z ), is similar as 6, withsummation taken on the data in Table II. We must take an average of Φ I ( z ) and Φ II ( z ) as an average energy ofeach atom in case of AB ordering.Figure 3 illustrates quick convergence of potential en-ergy Φ as a function of the total number of atoms N in asymmetric ”molecule” for type I and type II ”molecules”computed for LJ potential. An approximately 1 /N scal-ing, when the rings of equi-distant neighbours are addedin Eq. 6.Hence, an average potential energy for an atom at dis-tance z from the surface of bulk graphite is given as asum, for AA, AB, and ABC stackings, respectively: E AA ( z ) = ∞ X i =1 Φ I ( z + ic ) . (7) E AB ( z ) = ∞ X i =1
12 (Φ I ( z + (2 i − c ) + Φ II ( z + (2 i − c )) + ∞ X i =1 Φ I ( z + 2 ic ) , (8) E ABC ( z ) = 12 ∞ X i =1 [Φ I ( z + 3 ic ) + Φ II ( z + 3 ic )]+ [Φ I ( x + (3 i − c ) + Φ II ( x + (3 i − c )]+ [Φ I ( z + (3 i − c ) + Φ II ( z + (3 i − c )] . (9)In our numerical computation we limit the number ofplanes in Eqs 7–9 to 4, since contribution to potential energy from next planes diminishes quickly: it is of theorder of 90%, 10%, 2% and 0.5% for the first, second,4
49 49.5 50 50.5 51 51.5 52 52.5 53 0 0.00005 0.0001 0.00015 0.0002type IItype I -F [ m e V / a t o m ] FIG. 3: Scaling of potential energy depth as a function of1 /N , the number of atoms in a symmetric ”molecule” fortype I and II molecules, computed for LJ potential with pa-rameters ǫ = 2 . meV and σ = 3 .
38 ˚ A . Asymptotic splittingfor N → ∞ between type I and type II energy minima,(Φ II − Φ I ) / Φ II , is 0 . third and forth plane, respectively (Fig. 4). Considering4 planes only is equivalent to assuming cut-off distanceof 5 c , that is 16.7 ˚A.The proposed method of computing interlayer poten-tial is convenient for quick testing of potentials other thanthat of Mie or Lennard-Jones as well, since that requiresonly replacement of the function U ( z ) in Eq. 6 by an-other one.In case of perfect AA or AB stacking, due to symmetry, π -orbitals must point out in directions normal to planes.In that case, in Kolmogorov-Crespi Equations 4, valuesof n l r and n m r reduce to z , where z is distance of anatom from the plane. That is, ρ lm and ρ ml become equalto d from Tables I and II. We may rewrite Equations 4in this form: U KC ( x ) = − A (cid:16) z r (cid:17) + exp ( − λ ( r − z )) · " C + 2 exp (cid:0) − ( d/δ ) (cid:1) · X n =0 C n ( d/δ ) n , (10)Similarly, we may rewrite Lebedeva’s version of potential: U Leb ( x ) = − A (cid:16) z r (cid:17) + Bexp ( − α ( r − z )) + C (cid:0) D d + D d (cid:1) exp (cid:0) − λ d (cid:1) exp (cid:0) − λ (cid:0) z − z (cid:1)(cid:1) . (11)In 10 and 11, r = √ z + d . TABLE III: Summary of parameters’ values used in Kolmogorov-Crespi equation 10.A [meV] C [meV] C [ meV ] C [ meV ] C [ meV ] z [˚A] λ [˚A − ] δ [˚A] Ref.10.238 3.030 15.71 12.29 4.933 3.34 3.629 0.578 [42],[47]9.89 0.7 5.255 2.336 1.168 3.449 3.2 2.1 this work, set A8.881 0.6286 4.719 2.098 1.049 3.487 3.4 1.7 this work, set BTABLE IV: Summary of parameters’ values used in Lebedeva equation 11.A [meV] B [meV] C [meV] α [˚ A − ] z [˚A] λ [˚ A − ] λ [˚A − ] D [˚A − ] D [˚A − ] Ref.10.510 11.652 35.883 4.16 3.34 0.48703 0.46445 -0.86232 0.10049 [48]14.558 21.204 1.8 4.16 3.198 0.6 0.4 -0.862 0.10049 this work In case of Lennard-Jones potential, results presented here are computed with ǫ = 2 . meV and σ = 3 . E p [ m e V / a t o m ] z [Å] FIG. 4: Contribution to potential energy of a single atom(represented by large red dot in the insert) as a functionof distance z from the surface of three neighbouring planes,separated also by distance z . The case of ABA and ABCordering of planes is considered. The atom is consideredin two positions with respect to the first plane, AA(1) orAB(1). The energy contribution from the third plane isexactly the same regardless whether it is A plane or Cplane. Energy from the second and third planes are rescaled5 and 10 times, respectively. Classical L-J potential 12-6with parameters ǫ = 2 . meV and σ = 3 .
38 ˚ A was used.The depth of potential well for the lowest curve denotedas ”AB(1)+AA(2)+AB(3)+AA(4)+AB(5)” (which shows thesum of energies from 5 consecutive planes in AB configura-tion) is 52 meV. For models of Kolmogorov-Crespi and Lebedeva et al.we used values of parameters as listed in tables III andIV, respectively[77]. One ought to be aware that usuallythere is a broad range of parameters values for any ofthe above models of potential that may provide a func-tionally nearly identical dependencies U ( z ) and reason-able agreement with experiments; compare for instance[42],[47] with [41] and [45]. The original parameters ofLebedeva equation were derived by fitting to results of ab − initio DFT modeling and with the purpose of de-scribing graphene corrugation experiments. While werecognize the usefulness of proposed new equations forproviding phenomenological description of certain mate-rials properties, we find no strong justification for adher-ing to original values of parameters. In particular, as wediscuss later, the value of C in Table IV proposed in theoriginal work is likely significantly too large.The meaning of some parameters is as follows. σ and z decide most about position of potential minimum (equi-librium interlayer spacing), which is also sensitive to D and λ , while λ , δ , λ and mainly C in Eq. 11 decidemost about AB-AA energy difference at potentials min-ima. Value of λ influences mostly the slope of dE/dz inFig. 4 on low-z side of potential curve and as a conse-quence it decides about the compressibility as a function of pressure. Compressibility.
Equation 8 may be used for computing compressibilityalong c-axis, since force and pressure acting on a particleis proportional to derivative of potential energy: P ( c ) = η · dE AB ( c ) dc , (12)where c is lattice constant perpendicular to the planesunder pressure P and η is the number of atoms per unitarea, η = 61 . GP a · eV − · ˚ A − for graphene. The first(Eq. 12) and second derivatives of E AB ( z ) could be com-puted by using derived analytical expressions. Numericalapproach of finding derivatives is however convenient toimplement and fast.Figure 5 shows the found ratio of c/c as a function ofpressure for LJ, KC and Lebedeva potentials, with pa-rameters as these in Tables III and IV. In case of LJmodel, the parameter ǫ used was 3 . meV , which leads tobinding energy of 64 meV /atom , larger than most com-monly accepted value of around 50 meV /atom but leadingto better agreement between results of other calculationsand measurements, as reported in the following sections.The straight solid line with a slope of − . /GP a at P = 0 for LJ potential corresponds to elasticity mod-ulus C = 38 . GP a , and for KC one the slope gives C = 43 . GP a , in a very good agreement with experi-ment. Ab − initio DFT computations result in a broaderspread of values, from 43 . . GP a [54].Pressure dependence of c can be found or deduced fromseveral measurements ([55], [56], [57], [58], [53]). Exper-imental results in Fig. 5 have a broad dispersion of datapoints, which is, in part, due to subtle differences in theused measurement techniques. Stacking order
There is much controversy around prevailing stackingorder in graphite and in a few layers graphene (FLG).It is known that there are three types of stacking,an AA one, which is simple hexagonal (it is not foundin natural graphite and exists only in intercalated com-pounds such as C Li and C K , an AB stacking, hexag-onal, known also as Bernal, and a rhombohedral ABCstacking. Additionally, a random stacking of these threetypes is called a turbostratic TS structure and often is ob-tained in laboratories. In bulk graphite, it was reportedthat the volume fraction AB:ABC:TS was about 80:14:6([60]).Lee et al. [61] were able to obtain AA structuredgraphene on diamond surface, with an interplanar spac-6 c / c P [GPa] Lennard−JonesKolmogorov−Crespi AKolmogorov−Crespi BLebedevaA)B)C)D)
FIG. 5: Compressibility of graphite along c-axis. The solidcurves are computed numerically by using Eq. 12 with L-J,K-C, and Lebedeva potentials. Datapoints are collected fromworks of Lynch and Drickamer ([56], A), Zhao and Spain ([57],B), Hanfland et al. ([59], C), Clark et al. ([58], D). ing of ∼ . . SiC substrate and interpret their results in terms of possible modification ofsecond-plane interactions by the substrate, as explainedby Slonczewski-Weiss-McClure model [64]. Yoshizawa etal. [65] propose that the AB stacking of layers in graphiteis a consequence of orbital interactions between layersrather than of generally accepted van der Waals forces.The energy difference between AA and AB is reported,to be 17.31 meV/atom ([31]), in favour of AB, while be-tween ABC and AB it was reported as 0.11 meV/atom([66]), which are both based however on DFT calcula-tions, only.It was observed that in the case of graphene, its bilayerexhibits AB stacking while trilayer prefers ABC stacking[67]. When the number of layers increases, again ABstacking is favored. The binding energies are found toincrease from 23 meV /atom in bilayer to 39 meV /atom inpentalayer, while the interlayer distance decreases from3 . A in bilayer to 3 . A in pentalayer. Our models donot reproduce so strong change of binding energy withnumber of layers. In 2-layer graphene we find around90% of binding energy per atom of that in bulk graphite,when LJ model is used (Fig. 4). However, the change ofinterlayer distance agrees with that reported [68], for allthree models discussed. TABLE V: Energies of FLG with number of layers from 2 to 6, for AA, AB and ABC stacking.ordering 2 layers 3 layers 4 layers 5 layers 6 layersAA AA(1) 2AA(1)+AA(2) 3AA(1)+2AA(2)+AA(3) 4AA(1)+3AA(2)+2AA(3) 5AA(1)+3AA(2)+3AA(3)+2AA(4)AB AB(1) 2AB(1)+AA(2) 3AB(1)+2AA(2)+AB(3) 4AB(1)+3AA(2)+2AB(3) 5AB(1)+4AA(2)+3AB(3)+2AA(4)ABC – 2AB(1)+AA(2) 3AB(1)+2AB(2)+AA(3) 4AB(1)+3AB(2)+2AA(3) 5AB(1)+4AB(2)+3AA(3)+2AB(4)
If to assume that interaction between next nearestplanes (NNP) is negligible than there is no reason fordifference in stability of ABC and AB structures. Hence(NNP) interaction must play a role. Table V summarizestotal potential energy of a FLG. The notation used thereis the same as in Fig. 4, and it should be understood asfollows (taking AB stacking for 3 layers, as an example):2AB(1)+AA(2) means that we have 2 pairs of layers inAB stacking that are apart for 1 c , and one pair of AAordering of layers apart for 2 c . One should notice thatenergy of AC stacking is the same as that of AB, sinceatoms configuration of neighboring planes are the samein both cases. We see from Table V that for 4 and morelayers ABC stacking is energetically favorable. Howeverthe energy difference between AB and ABC is caused byinteraction of layers that are 2 c apart and therefore it isvery small, of around 3 µeV /atom .The explanation to AB:ABC ratio observed, of about 80:14 ([60]), could be found by adding to already exist-ing isotropic attractive term of vdW potential a smallanisotropic contribution to equations like these of KC orLebedeva. Since interplane binding is caused by stronglyanisotropic interaction between π orbitals one would ex-pect that not only repulsive but also attractive compo-nent of that interaction is anisotropic. Raman shifts in LBM and shear modes ( C and C elastic constants) Early neutron scattering measurements confirmed theexistence of the low-frequency acoustic mode in bulk py-rolitic graphite [69], that is longitudinal waves in thedirection of the hexagonal axis, at 3 . ± . T Hz .The transverse waves (in shear mode) were less pro-nounced, at 1 . ± . T Hz . Based on these data elas-7ic constants have being deduced, C = 39 ± GP a (fordirection perpendicular to graphite planes; LBM), and C = 4 . ± GP a (shear mode parallel to planes). Thesevalues correspond to Raman shift energies of 130 and 43 cm − , in agreement with several DFT calculations [70].The experimental evidence of the longitudinal mode wasreported in ultrafast laser pump-probe spectroscopy [71]on FLG, at ∼ cm − , demonstrating also the presencethere of the phonon and electron coupling through theBreit-Wigner-Fano resonance, as at more pronounced G-mode [72].Position of peaks in both LBM and shear modes is ageometrical effect and depends strongly on the number oflayers in FLG. Moreover, unique multipeak features areobserved, characteristic for the number of layers investi-gated. Their frequencies in case of LBM [73] mode aredescribed well using a simple linear-chain model basedon nearest-neighbor couplings between the layers [74]: ω N ( n ) = ω sin [( N − n ) π/ N ] , (13)where ω is the frequency of the bulk mode, N it thenumber of layers in FLG and n enumerates observed Ra-man frequencies. It follows from measurements that inthe limit of large number of layers the bulk LBM Ramanpeak should have value of 264 . cm − [73].For shear mode, the nature of interactions betweenplanes, their collective behaviour, leads to quantitativelysimilar dependence of Raman frequencies on the numberof layers in FLG as in the case of LBM mode: the fre-quency in the bulk is √ σ = 3 . A and ǫ = 2 . meV . We find that the en-ergy difference between planes in AA and AB positionsis of around 1 meV /atom , only, at P = 0 GP a .By having the curvature of parabolic potential wellsat low-values of departure from optimal position of bothlayers (in AB stacking; we used data for x less then 0.05˚A) we are able to compute Raman frequencies in shearmode in a broad range of pressure. The curvature k ofparabolic fit of data in Figure 6, E p ( x ) − E AB = k/ · x ,with equation R [ cm − ] = √ · √ k · k is in eV / ˚ A and R is the Raman shift. The √ P = 0 GP a is about twice too low to explain the Ra-man shift observed in experiment . This means that theenergy difference between AA and AB arrangement ofplanes must be around 5 meV per atom at P = 0 GP a ,which is significantly smaller than reported in some DFTcalculations [31], [48]. E p − E AB [ m e V / a t o m ] x [Å]13.910.88.26.34.63.12.061.230.50 FIG. 6: Datapoints show change in potential energy of a planethat is moved from AB position over another large plane, fordistance equal to a=1.42 ˚A, which results in AA position, un-der several values of pressure [GPa], as listed in the legend.The classical 12-6 Lennard-Jones potential was used in mod-elling, with σ = 3 . A and ǫ = 2 . meV . Solid lines show atentative fit of parabolic dependencies. For LBM mode results on the Raman shift are knownat P = 0 GP a , only, and our model predicts that valueaccurately [73]. In this case there is no need to useLAMMPS in computation. The scheme described in thissection offers an alternate and convenient way of numer-ical computing frequency of oscillations from potentialenergy curvature by starting with Eq. 8. We need to usealso the data available in Fig. 5 and find out from therehow planes’ spacing c changes with pressure. We insert c values into Eq. 8 and numerically find out the secondderivative d E AB ( z ) /dz , which gives us curvature of po-tential minimum and, next, frequencies of Raman shift,as a function of pressure.Figure 8 shows pressure dependence of Raman shift inLBM mode computed for LJ, KC and Lebedeva mod-els, by using this method. If similar experimental datawere available we could have an additional indication onwhich of considered models of vdW interactions fits bestto reality.8 R a m an s h i ft [ c m − ] P [GPa] Lennard−JonesKolmogorov−Crespi AKolmogorov−Crespi BLebedevaHanfland et al.
FIG. 7: Raman shift energy as a function of pressure, forbulk shear mode. The lines are computed from curvature k ofparabolic fit of data in Figure 6. Blue circles are experimentaldata of Hanfland et al. [59]. R a m an s h i ft [ c m − ] P [GPa]Lennard−JonesKolmogorov−Crespi AKolmogorov−Crespi BLebedeva
FIG. 8: Predicted pressure dependence of Raman shift forLBM mode of bulk graphite.
SUMMARY AND CONCLUSIONS.
Interlayer interaction modeling in graphenites basedon phenomenological van der Waals potentials (Lennard-Jones, Kolmogorov-Crespi and Lebedeva’s was per-formed. Results have been compared self-consistentlywith ab-initio calculations and experimental data oncompressibility and Raman shifts in LBM and shearmodes under pressure, favoring strongly anisotropicKolmogorov-Crespi model. Computation was done byusing molecular dynamics package LAMMPS and a pro-posed convenient, extendable scheme of computationsuitable for fast numerical modeling of several physicalquantities. The method is useful for studying other 2-dimensional homo- and heterostructures with van derWaals type interaction between layers. It is argued that the value of the known Raman shift in shear mode is con-sistent with the difference in energy between AA and ABstacking of around 5 meV per atom. The models do notprovide explanation for the reported low content of ABCstacking in natural graphite. [1] A. K. Geim and I. V. Grigorieva. Vander waals heterostructures.
Nature , 499:419–425, 2013. doi: 10.1038/nature12385. URL http://dx.doi.org/10.1038/nature12385 .[2] Nam B. Le, Tran Doan Huan, and Lilia M.Woods. Interlayer interactions in van der waalsheterostructures: Electron and phonon proper-ties.
ACS Appl. Mater. Interfaces , 9(9):6286–6292, 2016. doi: 10.1021/acsami.6b00285. URL http://dx.doi.org/10.1021/acsami.6b00285 .[3] M M van Wijk, A Schuring, M I Katsnel-son, and A Fasolino. Relaxation of moir´e pat-terns for slightly misaligned identical lattices:graphene on graphite.
2D Materials , 2(3):034010,2015. doi: 10.1088/2053-1583/2/3/034010. URL http://doi.org/10.1088/2053-1583/2/3/034010 .[4] M M van Wijk, A Schuring, M I Katsnelson, and A Fa-solino. Moir´e patterns as a probe of interplanar interac-tions for graphene on h-bn.
Phys. Rev. Lett. , 113:135504,2014. doi: 10.1103/PhysRevLett.113.135504. URL http://doi.org/10.1103/PhysRevLett.113.135504 .[5] Yuan Cao, Valla Fatemi, Shiang Fang, Kenji Watan-abe, Takashi Taniguchi, Efthimios Kaxiras, and PabloJarillo-Herrero. Unconventional superconductivity inmagic-angle graphene superlattices.
Nature , 2018.doi: http://dx.doi.org/10.1038/nature26160. URL http://dx.doi.org/10.1038/nature26160 .[6] Han Peng, Niels B. M. Schr¨oter, Jianbo Yin, HuanWang, Ting-Fung Chung, Haifeng Yang, Sandy Eka-hana, Zhongkai Liu, Juan Jiang, Lexian Yang,Teng Zhang, Heng Ni Cheng Chen, Alexey Bari-nov, Yong P. Chen, Zhongfan Liu, Hailin Peng, andYulin Chen. Substrate doping effect and unusuallylarge angle van hove singularity evolution in twistedbi- and multilayer graphene.
Adv. Mater. , page1606741, 2017. doi: 10.1002/adma.201606741. URL http://doi.org/10.1002/adma.201606741 .[7] Giacomo Argentero, Andreas Mittelberger, MohammadReza Ahmadpour Monazam, Yang Cao, Timothy J. Pen-nycook, Clemens Mangler, Christian Kramberger, JaniKotakoski, A. K. Geim, , and Jannik C. Meyer. Unravel-ing the 3d atomic structure of a suspended graphene/hbnvan der waals heterostructure.
Nano Lett. , 17:1409–1406,2017. doi: http://doi.org/10.1021/acs.nanolett.6b04360.URL http://doi.org/10.1021/acs.nanolett.6b04360 .[8] Shisheng Lin, Yanghua Lu, Juan Xu, SiruiFeng, and Jianfeng Li. High performancegraphene/semiconductor van der waals heterostruc-ture optoelectronic devices.
Nano Energy , 40:122–148,2017. doi: 10.1016/j.nanoen.2017.07.036. URL https://doi.org/10.1016/j.nanoen.2017.07.036 .[9] Grigory Skoblin, Jie Sun, and August Yurgens. Graphenebolometer with thermoelectric readout and capacitivecoupling to an antenna.
Appl. Phys. Lett. , 112:063501, https://doi.org/10.1063/1.5009629 .[10] Xia Wei, Fa-Guang Yan, Chao Shen, Quan-Shan Lv,and Kai-You Wang. Photodetectors based on junc-tions of two-dimensional transition metal dichalco-genides. Chinese Phys. B 26 , 26 (3):038504,2017. doi: 10.1088/1674-1056/26/3/038504. URL https://doi.org/10.1088/1674-1056/26/3/038504 .[11] Mu Congpu, Xiang Jianyong, and Liu Zhongyuan. Pho-todetectors based on sensitized two-dimensional transi-tion metal dichalcogenides—a review.
Journal of Ma-terials Research , 32(22):4115–4131, 2017. doi: 10.1557/jmr.2017.402.[12] M Houssa, A Dimoulas, and A Molle. Sil-icene: a review of recent experimental andtheoretical investigations.
Journal of Physics:Condensed Matter , 27(25):253002, 2015. URL http://stacks.iop.org/0953-8984/27/i=25/a=253002 .[13] D. G. Kvashnin and L. A. Chernozatonskii. Elec-tronic and transport properties of heterophase com-pounds based on mos2.
Jetp Lett. , 105:250,2017. doi: 10.1134/S0021364017040117. URL http://doi.org/10.1134/S0021364017040117 .[14] M. Luo, Y. H. Shen, and T. L. Yin. Struc-tural, electronic, and magnetic properties of transi-tion metal doped res2 monolayer.
Jetp Lett. , 105:255, 2017. doi: 10.1134/S0021364017040038. URL http://doi.org/10.1134/S0021364017040038 .[15] Matthew Z. Bellus, Ming Li, Samuel D. Lane, FrankCeballos, Qiannan Cui, Xiao Cheng Zeng, and HuiZhao. Type-i van der waals heterostructure formedby mos2 and res2 monolayers.
Nanoscale Horiz. , 2:31–36, 2017. doi: 10.1039/C6NH00144K. URL http://doi.org/10.1039/C6NH00144K .[16] Faguang Yan, Lixia Zhao, Amalia Patan`e, PingAn Hu,Xia Wei, Wengang Luo, Dong Zhang, Quanshan Lv,Qi Feng, Chao Shen, Kai Chang, Laurence Eaves, andKaiyou Wang. Fast multicolor photodetectors basedon graphene-contacted p-gase/n-inse van der waalsheterostructures.
Nanotechnology , 28(7):27LT01, 2017.doi: http://doi.org/10.1088/1361-6528/aa749e. URL http://stacks.iop.org/0957-4484/28/i=27/a=27LT01 .[17] M. Luo, H. H. Yin, and J. H. Chu. Magnetic prop-erties of a na-doped ws2 monolayer in the presence ofan isotropic strain.
JETP Letters , 6(10):672–676, 2017.doi: http://doi.org/10.1134/S0021364017220039. URL http://doi.org/10.1134/S0021364017220039 .[18] Oriol Lopez Sanchez, Dmitry Ovchinnikov, ShikharMisra, Adrien Allain, and Andras Kis. Valley polariza-tion by spin injection in a light-emitting van der waalsheterojunction.
Nano Lett. , 16:5792–5797, 2016. doi:http://doi.org/10.1021/acs.nanolett.6b02527.[19] Mallikarjuna Gurram, Siddhartha Omar, and Bart J.van Wees. Bias induced up to 100ferromagnet/bilayer-hbn/graphene/hbn heterostructures.
Nature Commu-nications , 8:248, 2017. doi: http://doi.org/10.1038/s41467-017-00317-w.[20] Jaewoo Shim, Seo-Hyeon Jo, Minwoo Kim, Young JaeSong, Jeehwan Kim, and Jin-Hong Park. Light-triggeredternary device and inverter based on heterojunction ofvan der waals materials.
ACS Nano , 11(6):6319–6327,2017. doi: http://doi.org/10.1021/acsnano.7b02635.[21] Wei Li, Xinlian Wang, and Xianqi Dai. Tunable schot-tky contacts in the antimonene/graphene van der waals heterostructures.
Solid State Communications , 254:37–41, 2017. doi: 10.1016/j.ssc.2017.02.008. URL https://doi.org/10.1016/j.ssc.2017.02.008 .[22] K. H. Michel, D. C¸ akır, C. Sevik, and F. M. Peeters.Piezoelectricity in two-dimensional materials: Compar-ative study between lattice dynamics and ab initio cal-culations.
Phys. Rev. B , 95:125415, 2017. doi: https://doi.org/10.1103/PhysRevB.95.125415.[23] Shigeki Kawai, Adam S. Foster, Torbj¨orn Bj¨orkman, Syl-wia Nowakowska, Jonas Bj¨ork, Filippo Federici Canova,Lutz H. Gade, Thomas A. Jung, and Ernst Meyer.Van der waals interactions and the limits of isolatedatom models at interfaces.
Nature Communications ,7:11559, 2016. doi: 10.1038/ncomms11559. URL http://dx.doi.org/10.1038/ncomms11559 .[24] A.V. Rozhkov, A.O. Sboychakov, A.L. Rakhmanov,and Franco Nori. Electronic properties of graphene-based bilayer systems.
Physics Reports , 648:1–104,2016. doi: 10.1016/j.physrep.2016.07.003. URL http://dx.doi.org/10.1016/j.physrep.2016.07.003 .[25] L. A. Girifalco and Miroslav Hodak. Van der waals bind-ing energies in graphitic structures.
Phys. Rev. B , 65:125404, 2002. doi: 10.1103/PhysRevB.65.125404. URL http://doi.org/10.1103/PhysRevB.65.125404 .[26] J.-C. Charlier, X. Gonze, and J.-P. Michenaud.Graphite interplanar bonding: Electronic de-localization and van der waals interaction.
Europhysics Letters , 28:403, 1994. URL http://stacks.iop.org/0295-5075/28/i=6/a=005 .[27] S. B. Trickey, F. M¨uller-Plathe, G. H. F. Diercksen,, and J. C. Boettger. Interplanar binding and latticerelaxation in a graphite dilayer.
Physical Review B ,46:4460, 1992. doi: 10.1103/PhysRevB.45.4460. URL http://doi.org/10.1103/PhysRevB.45.4460 .[28] H. Rydberg, M. Dion, N. Jacobson, E. Schr¨oder,P. Hyldgaard, S. I. Simak, D. C. Langreth, and B. I.Lundqvist. Van der walls density functional for lay-ered structures.
Physical Review Letters , 91:126402,2003. doi: 10.1103/PhysRevLett.91.126402. URL https://doi.org/10.1103/PhysRevLett.91.126402 .[29] D. P. DiVincenzo, E. J. Mele, and N. A. W.Holzwarth. Density-functional study of interpla-nar binding in graphite.
Physical Review B , 27:2458, 1983. doi: 10.1103/PhysRevB.27.2458. URL https://doi.org/10.1103/PhysRevB.27.2458 .[30] Matthias C. Schabel and Jos´e Lu´ıs Martins. Energeticsof interplanar binding in graphite.
Physical Review B ,46:7185, 1992. doi: 10.1103/PhysRevB.46.7185. URL https://doi.org/10.1103/PhysRevB.46.7185 .[31] E. Mostaani, N. D. Drummond, and V. I. Fal’ko.Quantum monte carlo calculation of the binding en-ergy of bilayer graphene.
Phys. Rev. Lett. , 115:115501,2015. doi: 10.1103/PhysRevLett.115.115501. URL https://doi.org/10.1103/PhysRevLett.115.115501 .[32] M. Birowska, K. Milowska, and J.A. Ma-jewski. Van der waals density function-als for graphene layers and graphite.
ActaPhysica Polonica A , 120(5):845, 2011. URL http://przyrbwn.icm.edu.pl/APP/PDF/116/a116z520.pdf .[33] S. D. Chakarova-Kack, E. Schroder, B. I. Lundqvi-stand, and D. C. Langreth. Application of vander waals density functional to an extendedsystem: Adsorption of benzene and naphtha-lene on graphite.
Phys. Rev. Lett. , 96:146107, http://doi.org/10.1103/PhysRevLett.96.146107 .[34] Jian Ping Lu, X.-P. Li, , and Richard M. Martin. Groundstate and phase transitions in solid c60. Phys. Rev. Lett. ,68:1551, 1992. doi: 10.1103/PhysRevLett.68.1551. URL https://doi.org/10.1103/PhysRevLett.68.1551 .[35] L.Y. Jiang, Y. Huang, H. Jiang, G. Ravichan-dran, H. GaoandK.C. Hwang, and B. Liu. A co-hesive law for carbon nanotube/polymer interfacesbased on the van der waals force.
Journal ofthe Mechanics and Physics of Solids , 54:2436–2452,2006. doi: 10.1016/j.jmps.2006.04.009. URL http://doi.org/10.1016/j.jmps.2006.04.009 .[36] X.Q. He, S. Kitipornchai, and K.M. Liew. Bucklinganalysis of multi-walled carbon nanotubes: a contin-uum model accounting for van der waals interaction.
Journal of the Mechanics and Physics of Solids , 53:303–326, 2005. doi: 10.1016/j.jmps.2004.08.003. URL http://doi.org/10.1016/j.jmps.2004.08.003 .[37] S. Kitipornchai, X. Q. He, and K. M. Liew.Continuum model for the vibration of multilay-ered graphene sheets.
Phys. Rev. B , 72:075443,2005. doi: 10.1103/PhysRevB.72.075443. URL http://doi.org/10.1103/PhysRevB.72.075443 .[38] H. Tan, L.Y. Jiang, Y. Huang, B. Liu, and K.C. Hwang.The effect of van der waals-based interface cohesivelaw on carbon nanotube-reinforced composite materials.
Composites Science and Technology , 67(14):2941–2946,2007. doi: 10.1016/j.compscitech.2007.05.016. URL https://doi.org/10.1016/j.compscitech.2007.05.016 .[39] Gustav Mie. Zur kinetischen theorie dereinatomigen k¨orper.
Ann Phys , 316:657,1903. doi: 10.1002/andp.19033160802. URL http://dx.doi.org/10.1002/andp.19033160802 .[40] C. Avendano, T. Lafitte, A. Galindo, C. S. Adji-man, G. Jackson, and E. Muller. Saft- γ force fieldfor the simulation of molecular fluids. 1. a single-sitecoarse grained model of carbon dioxide. J Phys ChemB , 115:11154, 2011. doi: 10.1021/jp204908. URL http://dx.doi.org/10.1021/jp204908 .[41] Aleksey N. Kolmogorov and Vincent H. Crespi.Smoothest bearings: interlayer sliding in multiwalledcarbon nanotubes.
Phys. Rev. Lett. , 85(22):4727–4730, 2000. doi: 10.1103/PhysRevLett.85.4727. URL http://doi.org/10.1103/PhysRevLett.85.4727 .[42] Aleksey N. Kolmogorov and Vincent H. Crespi.Registry-dependent interlayer potential for graphiticsystems.
Physical Review B , 71(23):235415,2005. doi: 10.1103/PhysRevB.71.235415. URL https://doi.org/10.1103/PhysRevB.71.235415 .[43] S. Plimpton. Fast parallel algorithms for short-rangemolecular dynamics.
J Comp Phys, 117, 1-19 (1995) ,117:1–19, 1995. URL http://lammps.sandia.gov .[44] M. M. van Wijk, A. Schuring, M. I. Katsnel-son, and A. Fasolino. Moir´e patterns as aprobe of interplanar interactions for grapheneon h-bn.
Phys. Rev. Lett. , 113:135504, 2014.doi: 10.1103/PhysRevLett.113.135504. URL http://dx.doi.org/10.1103/PhysRevLett.113.135504 .[45] Annelot Schuring. Moir´e patterns in graphene in hexag-onal substrates.
Radbound University Niimegen , MasterThesis, 2014.[46] Tuck Wah Ng, Chun Yat Lau, Esteban Bernados-Chamagne, Jefferson Zhe Liu, John Sheridan, and Ne Tan. Graphite flake self-retraction response basedon potential seeking.
Nanoscale Res Lett. , 7(1):185, 2012. doi: 10.1103/PhysRevB.93.235414. URL http://doi.org/10.1103/PhysRevB.93.235414 .[47] Roberto Guerra, Ugo Tartaglino, Andrea Vanossi, andErio Tosatti. Ballistic nanofriction.
Nature Materi-als , 9:634–637, 2010. doi: 10.1038/nmat2798. URL http://doi.org/10.1038/nmat2798 .[48] Irina V. Lebedeva, Andrey A. Knizhnik, Andrey M.Popov, Yurii E. Lozovik, and Boris V. Potapkin. Mod-eling of graphene-based nems.
Physica E , 44(6):949–954, 2012. doi: 10.1016/j.physe.2011.07.018. URL https://doi.org/10.1016/j.physe.2011.07.018 .[49] Irina V. Lebedeva, Andrey A. Knizhnik, Andrey M.Popov, Yurii E. Lozovik, and Boris V. Potapkin.Molecular dynamics simulation of the self-retractingmotion of a graphene flake.
Phys. Rev. B , 84:245437, 2011. doi: 10.1103/PhysRevB.84.245437. URL https://doi.org/10.1103/PhysRevB.84.245437 .[50] Jin-Wu Jiang. Registry effect on the thermal con-ductivity of few-layer graphene.
Journal of AppliedPhysics , 116:164313, 2014. doi: 10.1063/1.4900526. URL http://dx.doi.org/10.1063/1.4900526 .[51] V. Baskin and L. Meyer. Lattice constants ofgraphite at low temperatures.
Phys. Rev. B , 100(2):544, 1955. doi: 10.1103/PhysRevB.100.544. URL http://doi.org/10.1103/PhysRevB.100.544 .[52] C. H. Kiang, M. Endo, P. M. Ajayan, G. Dres-selhaus, and M. S. Dresselhaus. Size effects incarbon nanotubes.
Phys. Rev. Lett. , 81(9):1869,1998. doi: 10.1103/PhysRevLett.81.1869. URL http://doi.org/10.1103/PhysRevLett.81.1869 .[53] P. Trucano and R. Chen. Structure of graphite by neu-tron diffraction.
Nature , 258:136, 1975.[54] Wei Gao and Rui Huang. Effect of surface rough-ness on adhesion of graphene membranes.
J.Phys. D: Appl. Phys. , 44:452001, 2011. URL http://stacks.iop.org/0022-3727/44/i=45/a=452001 .[55] Yuejian Wang, Joseph E. Panzik, Boris Kiefer, andKanani K. M. Lee. Crystal structure of graphite un-der room-temperature compression and decompression.
Scientific Reports , 2:520, 2012. doi: 10.1038/srep00520.URL http://doi.org/10.1038/srep00520 .[56] R. W. Lynch and H. G. Drickamer. Effect of high pres-sure on the lattice parameters of diamond, graphite,and hexagonal boron nitride.
The Journal of Chemi-cal Physics , 44:181, 1966. doi: http://doi.org/10.1063/1.1726442. URL http://doi.org/10.1063/1.1726442 .[57] You Xiang Zhao and Ian L. Spain. X-ray diffrac-tion data for graphite to 20 gpa.
Phys. Rev. B ,40(2):993, 1989. doi: 10.1103/PhysRevB.40.993. URL http://doi.org/10.1103/PhysRevB.40.993 .[58] S.M. Clark, Ki-Joon Jeon, Jing-Yin Chen, and Choong-Shik Yoo. Few-layer graphene under high pressure: Ra-man and x-ray diffraction studies.
Solid State Communi-cations , 154:15–18, 2013. doi: 10.1016/j.ssc.2012.10.002.URL http://doi.org/10.1016/j.ssc.2012.10.002 .[59] M. Hanfland, H. Beister, and K. Syassen. Graphiteunder pressure: Equation of state and first-order raman modes.
Phys Rev B , 39:12598,1989. doi: 10.1103/PhysRevB.39.12598. URL http://doi.org/10.1103/PhysRevB.39.12598 .[60] H. Lipson and A.R. Stokes. The structure of graphite.
Proc. Roy. Soc. A , 181:101–105, 1942.
61] Jae-Kap Lee, Seung-Cheol Lee, Jae-Pyoung Ahn, Soo-Chul Kim, John I. B. Wilson, and Phillip John. Thegrowth of aa graphite on (111) diamond.
J. Chem.Phys. , 129:234709, 2008. doi: 10.1063/1.2975333. URL http://dx.doi.org/10.1063/1.2975333 .[62] J. R. Dahn, Rosamaria Fong, and M. J. Spoon.Suppression of staging in lithium-intercalated car-bon by disorder in the host.
Phys. Rev. B , 42:6424, 1990. doi: 10.1103/PhysRevB.42.6424. URL https://doi.org/10.1103/PhysRevB.42.6424 .[63] Wataru Norimatsu and Michiko Kusunoki. Se-lective formation of abc-stacked graphene lay-ers on sic(0001).
Phys. Rev. B , 81:161410,2010. doi: 10.1103/PhysRevB.81.161410. URL http://doi.org/10.1103/PhysRevB.81.161410 .[64] J.-C. Charlier, X. Gonze, and J.-P. Michenaud.First-principles study of the electronic prop-erties of graphite.
Phys. Rev. B , 43:4579,1991. doi: 10.1103/PhysRevB.43.4579. URL http://doi.org/10.1103/PhysRevB.43.4579 .[65] Kazunari Yoshizawa, Takashi Kato, and Tokio Yam-abe. Interlayer interactions in two-dimensional sys-tems: Second-order effects causing abab stacking oflayers in graphite.
J. Chem. Phys. , 105(5):2099,1996. doi: http://dx.doi.org/10.1063/1.472076. URL http://dx.doi.org/10.1063/1.472076 .[66] J. C. Charlier, X. Gonze, and J. P. Michenaud.First-principles study of the stacking effecton the electronic properties of graphite(s).
Carbon , 32 (2):289–299, 1994. doi: https://doi.org/10.1016/0008-6223(94)90192-9. URL https://doi.org/10.1016/0008-6223(94)90192-9 .[67] W. Gao and A. Tkatchenko. Sliding mecha-nisms in multilayered hexagonal boron nitride andgraphene: The effects of directionality, thickness, andsliding constraints.
Phys. Rev. Lett. , 114:096101,2015. doi: 10.1103/PhysRevLett.114.096101. URL https://doi.org/10.1103/PhysRevLett.114.096101 .[68] Wang Tao, Guo Qing, Liu Yan, and ShengKuang. A comparative investigation of anab- and aa-stacked bilayer graphene sheet un-der an applied electric field: A density func-tional theory study.
Chin. Phys. B , 21:067301,2012. doi: 10.1088/1674-1056/21/6/067301. URL http://doi.org/10.1088/1674-1056/21/6/067301 . [69] G. Dolling and B. N. Brockhouse. Lattice vi-brations in pyrolitic graphite.
Phys. Rev. , 128:1120, 1962. doi: 10.1103/PhysRev.128.1120. URL https://doi.org/10.1103/PhysRev.128.1120 .[70] Alexander V. Lebedev, Irina V. Lebedeva, An-drey M. Popov, and Andrey A. Knizhnik. Stackingin incommensurate graphene/hexagonal-boron-nitrideheterostructures based on ab initio study of in-terlayer interaction.
Phys. Rev. B , 96:085432,2017. doi: 10.1103/PhysRevB.96.085432. URL https://doi.org/10.1103/PhysRevB.96.085432 .[71] Jingzhi Shang, Chunxiao Cong, Jun Zhang, Qihua Xiong,Gagik G. Gurzadyan, and Ting Yu. Observation oflow-wavenumber out-of-plane optical phonon in few-layergraphene.
J Raman Spec , 40:70–74, 2013. doi: 10.1002/jrs.4141. URL http://doi.org/10.1002/jrs.4141 .[72] J. M. Baranowski, M. Mozdzonek, P. Dabrowski,K. Grodecki, P. Osewski, W. Kozlowski, M. Kopciuszyn-ski, M. Jalochowski, and W. Strupinski. Observation ofelectron-phonon couplings and fano resonances in epitax-ial bilayer graphene.
Graphene , 2:115–120, 2013. doi:http://dx.doi.org/10.4236/graphene.2013.24017. URL http://dx.doi.org/10.4236/graphene.2013.24017 .[73] Chun Hung Lui and Tony F. Heinz. Mea-surement of layer breathing mode vibrations infew-layer graphene.
Phys. Rev. B , 87:121404,2013. doi: 10.1103/PhysRevB.87.121404. URL http://doi.org/10.1103/PhysRevB.87.121404 .[74] Steven T. Thornton and Jerry B. Marion. Classical dy-namics of particles and systems.
Fifth Edition , BrooksCole, Pacific Grove, CA, 2003.[75] P. H. Tan, W. P. Han, 1 Z. H. Wu W. J. Zhao, K. Chang,H. Wang, Y. F. Wang, N. Bonini, N. Marzari, G. Savini,A. Lombardo, and A. C. Ferrari. The shear mode ofmultilayer graphene.
Nat. Mater. , 11:294–300, 2012. doi:http://doi.org/10.1038/nmat3245.[76] Chunxiao Cong and Ting Yu. Enhanced ultra-low-frequency interlayer shear modes in foldedgraphene layers.
Nature Communications , 5:4709, 2014. doi: 10.1038/ncomms5709. URL http://doi.org/10.1038/ncomms5709 .[77] LAMMPS implementation of equation of Lebedeva et al.is available from the corresponding author..[77] LAMMPS implementation of equation of Lebedeva et al.is available from the corresponding author.