A Torsional potential for graphene derived from fitting to DFT results
Georgios D. Chatzidakis George Kalosakas Zacharias G. Fthenakis Nektarios N. Lathiotakis
cond-mat.mtrl-sci



  

  
  Department of Physics, National Technical University of Athens, GR-15780 Athens, Greece 
  
  Materials Science Department, University of Patras, Rio GR-26504, Greece
  Crete Center for Quantum Complexity and Nanotechnology (CCQCN), Physics Department, University of Crete GR-71003 Heraklion, Greece
  
  Institute of Electronic Structure and Laser, FORTH, Heraklion, Greece
  Department of Physics, University of South Florida, Tampa, Florida 33620, USA
  
  Theoretical and Physical Chemistry Institute, National Hellenic
  Research Foundation, Vass. Constantinou 48, GR-11635 Athens, Greece  

  
  
  

We present a simple torsional potential for graphene to accurately describe its out-of-plane deformations.
The parameters of the potential are derived through appropriate fitting with suitable DFT calculations
regarding the deformation energy of graphene sheets folded around two different folding axes, along
an armchair or along a zig-zag direction. Removing the energetic contribution of bending angles,
using a previously introduced angle bending potential, we isolate the purely torsional deformation
energy, which is then fitted to simple torsional force fields. The presented out-of-plane torsional potential
can accurately fit the deformation energy for relatively large torsional angles up to 0.5 rad.
To test our proposed potential, we apply it to the problem of the vertical displacement of a single
carbon atom out of the graphene plane and compare the obtained deformation energy with
corresponding DFT calculations. The dependence of the deformation energy on the vertical displacement
of the pulled carbon atom is indistinguishable in these two cases, for displacements up to about 0.5 Å.
The presented potential is applicable to other sp^2 carbon structures.

A Torsional potential for graphene derived from fitting to DFT results
    Nektarios N. Lathiotakis
    December 30, 2023
======================================================================









§ INTRODUCTION



Following the isolation of single layer graphene <cit.> an enormous research effort has been
devoted to the study of this two-dimensional material and its properties <cit.>.
Potential applications have been explored in electronics <cit.>, opto-electronics <cit.>,
gas filtering <cit.>, energy storage <cit.>, uses related to its unique mechanical 
properties <cit.>, etc.

Many empirical force fields have been used in atomistic simulations, calculating various
structural, mechanical or phonon properties of graphene
<cit.>.
Besides the older, well known Tersoff <cit.> and Brenner <cit.> potentials, more
accurate force fields have been introduced the last two decades. For example, optimized parameter sets for the
latter potentials, providing better description of structural and phonon properties of graphene are presented
in Ref. <cit.>. LCBOP <cit.> and AIREBO <cit.> are efficient potentials
that have been widely applied in many calculations. Other potentials leading to good predictions of elastic
and thermal properties of graphene have been also discussed <cit.>.

More recently, we have presented simple analytical expressions for the accurate description of bond stretching and
angle bending potentials of graphene <cit.>. These potentials are derived by fitting analytical functions
to the deformation energy of proper distortions of graphene, obtained through accurate calculations
from first principles' methods (DFT). The presented force field is applicable only to distortions restricted
within the plane of graphene. These in-plane potentials can accurately describe elastic properties and the
mechanical response of graphene in various extensional loads <cit.>.
In this work, using similar ideas and methods, we extend this force field with torsional
energy terms, in order to be able to describe out-of-plane distortions in graphene. 
The basic motivation is to provide a simple and computational efficient classical potential which can be used
for accurate large-scale atomistic calculations. The torsional potential
presented here is also capable to describe other non-planar sp^2 carbon systems, like fullerenes
and carbon nanotubes <cit.>.

In the present work, we describe in detail the procedure followed and the necessary analytical calculations in order to fit the proposed
torsional potential to ab-initio data. 
The full potential is then tested in the case of the deformation energy due to the vertical displacement of a C atom outside graphene's plane. A more comprehensive benchmark study for fullerenes, nanotubes and graphene's phonons is presented elsewhere <cit.>.
We have considered two types of folding of graphene sheets around
different axes (either an armchair or a zig-zag one). The corresponding deformation energies are
calculated using DFT methods. Following the removal of the contribution of angle bending terms in the
total deformation energy, we isolated the pure torsional energy. Then, the analytic modeling of this energy in terms of 
individual torsional contributions, leads to a fitting procedure providing the optimal parameters of the
out-of-plane torsional energy.

This paper is organized as follows: In Sect. <ref> we describe the structures and methodology adopted for the
DFT simulations that was used to obtain the deformation energies. Then in Sect. <ref>, we
present the analytic work for removing the angle bending contributions from the deformation energies. 
The analytic expressions of the torsional energy terms in terms of the folding angles are provided in
Sec. <ref>. The fitting of the torsional terms of both 
Model 1 and 2 is described in Sec. <ref> completing the presentation of the derivation of the
new torsional force fields. Then, in Sec. <ref>, we present a test case, the deformation energy 
due to the vertical displacement of a carbon atom outside graphene's plane, as a
 first application of the proposed scheme and compare the prediction of Models 1 and 2 with DFT results. 
Finally, a summary and conclusions are given in Sect. <ref>.



§ STRUCTURES AND DFT CALCULATIONS




 §.§ Torsion angles in graphene


In Fig. <ref>, we show a part of  the honeycomb 
structure of graphene and a few carbon-atom positions labeled as i, j, k, l, m.  The quadruple
(i, j, k, l), with 3 of these positions belonging to the same hexagonal ring and one to an 
adjacent is 
customary called “trans" while the quadruple  (i, j, k, m), with all belonging to the same ring,
“cis". In the case that the structure is distorted and the atoms in the quadruple are no longer 
co-planar, we can define a torsion angle, which we label as (i-j-k-l),
as the dihedral angle of the planes i-j-k, j-k-l. The torsion angles can then be classified 
as “trans" or “cis" accordingly. 
The dihedral angle between two planes, e.g. i-j-k and j-k-l can be defined as the angle between
the vectors perpendicular to the planes. We assumed that the perpendicular vectors are pointing 
inwards for clockwise triples (like i-j-k) or outwards for anti-clockwise triples (like j-k-l). 
Under this assumption, torsional angles are in the range [0,π] with “cis" angles smaller
than π/2 and “trans" angles larger than π/2.






 §.§ DFT results


DFT calculations were performed for two distorted graphene structures, one that graphene is folded around an armchair
axis and one around a zig-zag. These structures are shown in Fig. <ref> where
we label all atoms relevant to the present discussion. They are periodic along the folding axis while on
the vertical they are not.
The  folding angle around either the armchair or zig-zag axis is denoted by ϕ. 
With  symbol θ, we denote “usual” angles between carbon bonds (bending angles) and
with ω torsion angles as defined  in the previous subsection.

All calculations were performed with Quantum-Espresso periodic-DFT code <cit.>, with the same
pseudopotential <cit.> as in Ref. kalosakas. The wave-function and density plane-wave
cutoffs were chosen 40 Ry and 400 Ry respectively. The unit cell we chose is minimal in the periodic
direction (that of the folding axis) while, in the vertical, it is appropriately large to avoid edge-effects. 
Thus, the simulated structures are nanoribbons that are folded around their middle line direction. 
In the case of the armchair folding (Fig. <ref>(top)), the vertical unit cell direction is such that 
neighbors up to the 5th in the vertical direction were included. In the case of the zig-zag folding that size is
long enough to include up to 8th neighbors. Thus, the unit cells contain 22 and 18 atoms for the armchair and
the zig-zag folding, respectively. In the reciprocal space, we used a mesh of
1×24×1, i.e. 24 points were assumed along the bending direction that the structure is periodic.





In Fig. <ref>, we show with filled circles the total deformation energy per unit cell along the folding direction,
as calculated by DFT, as a function of the folding angle ϕ, for both armchair  and zig-zag folding actions,
E_d^ (a) and E_d^ (z), respectively. The total deformation energy per unit-cell is taken as the
energy difference  between the folded structure and the not folded one (ϕ=0). 
Apparently, the two structure distortions due to the considered foldings are complex and consist of several individual
angle-bending and  torsional deformations. Note that bond lengths are not altered so there is no bond-stretching contribution 
in the total deformation energy. As we see in Fig. <ref>, the contribution from angle bending 
is significant for ϕ larger than 0.2 rad. In order to perform a fitting for the torsional terms alone, we first need 
to exclude angle-bending contributions from the total deformation energy. In order to do so, (i) we identify all 
angle-bending terms and express analytically their corresponding bending angles θ in terms of ϕ and then 
(ii) we remove the angle-bending terms using the analytic terms in the Ref. <cit.>. 
The residual, torsional energy per unit cell, when the contribution from angle-bending is subtracted, as a function of the out-of-plane
folding angle ϕ is shown  in Fig. <ref> (diamonds), for the two folding directions.
In order to fit an analytic expression to 
the torsional terms, we also have to (i) identify all the individual torsional terms that contribute 
for each of the zig-zag and armchair cases, and subsequently (ii) express the corresponding torsion
angles as functions of the folding angle, ϕ. 
These steps are described below where we provide all necessary analytical expressions. 



 

 


§ REMOVING THE ANGLE-BENDING TERMS





 §.§ For the folding around armchair axis


The folding around the armchair direction (Fig. <ref>, top)  alters 
two bond angles per unit cell, (32̂4),(61̂5), which are equal. One can show that 
these angles, in terms of ϕ, are given by 


    θ^ (a)=2arcsin(√(3/8)√(cosϕ+1))


The angle-bending energy that one needs to remove from the total energy is  

    U_b^ (a) = 2 V_b ( θ^ (a)(ϕ) ) ,

where V_b(θ) is the analytical expression for the angle bending given in Ref. <cit.>, i.e.

    V_b(θ) = k/2( θ-2π/3)^2 - k^'/3( θ-2π/3)^3 ,

with k=7.0   eV/ rad^2 and k^'=4   eV/ rad^3.

Removing these terms from the total deformation energies, E_d^ (a) we find the 
total torsional energy

    E_t^  (a) = E_d^ (a) - U_b^ (a) ,


 shown in Fig. <ref>. 




 §.§ For the folding around zig-zag axis


Similarly, the folding around the zig-zag direction (Fig. <ref>, (bottom)) 
affects two angles per unit cell, (21̂6),(21̂5), that are also equal. 
In terms of ϕ, these angles are given by 


    θ^ (z)=2arcsin(1/2√(2+cosϕ))


Again, the angle-bending energy that one needs to remove from the total energy is  

    U_b^ (z) = 2 V_b ( θ^ (z)(ϕ) ) ,

we remove these terms from the total deformation energy and we find the total torsional energy,

    E_t^  (z) = E_d^ (z) - U_b^ (z) ,

shown in Fig. <ref>.




§ ANALYTICAL EXPRESSIONS FOR THE TORSIONAL TERMS



Here we provide analytical expressions, U_t^ (a)(ϕ), U_t^ (z)(ϕ), for the total torsional 
energies, as functions of ϕ that will contain parameters to be fitted so that these expressions 
reproduce as close as possible the E_t^ (a), E_t^ (z) points shown in Fig. <ref>.
To arrive to such analytical expressions we need first to identify all altered torsion angles (per unit cell) 
and express them in terms of the folding angle ϕ. Then, U_t^ (a)(ϕ), U_t^ (z)(ϕ) will
be just the sum of all individual torsional terms that correspond to these altered torsional angles.

Regarding the individual torsional energy term, V_t(ω), two different functional forms would be considered.
The most frequently used formula, referred as Model 1 here, is 

    V_t(ω)    =   1/2 V_1 [ 1 + cosω] 
       +   1/2 V_2 [ 1 - cos (2 ω ) ]  .
 
An alternative model that we considered, which we call Model 2,
assumes a different fitting formula for cis or trans dihedral angles ω

    [     V_t^ (cis)(ω) = K_ cis[ 1 - cos (2 ω ) ] ,; V_t^ (trans)(ω) = K_ trans[ 1 - cos (2 ω ) ] , ]
 
where either the first or the second expression is used for cis or trans torsion angles, respectively.
Below we use both Models 1 and 2 to fit their parameters to the obtained DFT results.




 §.§ For the folding around armchair axis



Inspecting the Fig. <ref> (top) we identify the following
torsion (dihedral) angles per unit cell that are altered by folding along the armchair axis

  
  ∙  2 trans dihedral angles, (5-1-2-3), (4-2-1-6), with
  
    ω_1^ (a)(ϕ)=arccos(-cosϕ)

  
  ∙ 4 cis dihedral angles, (11-4-2-3), (14-5-1-6), (7-3-2-4), (10-6-1-5), with
  
    ω_2^ (a)(ϕ)=arccos(√(3)1+cosϕ/√(9sin^2ϕ+6(1+cosϕ)))

  
  ∙ 4 trans dihedral angles, (12-4-2-3), (13-5-1-6), (8-3-2-4), (9-6-1-5), with
  
    ω_3^(a)(ϕ)=arccos(-√(3)1+cosϕ/√(9sin^2ϕ+6(1+cosϕ)))



Through the angle expressions given above, the total torsional energy, U_t^ (a) within the Model 1,
 becomes an analytic function of ϕ:

    U_t^ (a)(ϕ)     =     2 V_t (ω_1^ (a)(ϕ)) + 4 V_t ( ω_2^(a)(ϕ) ) 
       +    4 V_t ( ω_3^(a)(ϕ) )

where V_t(ω) is the individual torsional term Eq. (<ref>).

For the Model 2, the corresponding expression of the total torsional energy U_t^ (a) is

    U_t^ (a)(ϕ)    =    2 V_t^ (trans) (ω_1^(a)(ϕ)) + 4 V_t^ (cis) ( ω_2^(a)(ϕ) ) 
       +    4 V_t^ (trans) ( ω_3^(a)(ϕ) )

with V_t^ (cis) and V_t^ (trans) given by Eq. (<ref>).




 §.§ For the folding around zig-zag axis



Inspecting the Fig. <ref> (bottom) we identify the following
dihedral angles per unit cell that are affected 


  
  ∙  2 cis dihedral angles, (3-2-1-6), (4-2-1-5), with
  
    ω_1^ (z)(ϕ)=arccos(√(3/sin^2ϕ+3))


  
  ∙ 2 trans dihedral angles, (4-2-1-6), (3-2-1-5), with
  
    ω_2^ (z)(ϕ)=arccos(-√(3/sin^2ϕ+3))


  
  ∙ 2 cis dihedral angles, (2-1-5-13), (2-1-6-9), with
  
    ω_3^ (z)(ϕ)=arccos(√(3/sin^2ϕ+3)cosϕ)


  
  ∙ 2 trans dihedral angles, (2-1-5-14), (2-1-6-10), with
  
    ω_4^ (z)(ϕ)=arccos(-√(3/sin^2ϕ+3)cosϕ)


   

And the total torsional energy is given by

    U_t^ (z)(ϕ)     =     2 V_t (ω_1^ (z)(ϕ)) + 2 V_t ( ω_2^ (z)(ϕ) ) 
       +    2 V_t ( ω_3^ (z)(ϕ)) + 2 V_t (ω_4^ (z)(ϕ)) ,

where V_t is given by Eq. (<ref>) for Model 1.
In case of Model 2, the above formula becomes 

    U_t^ (z)(ϕ)    =   2 V_t^ (cis) (ω_1^ (z)(ϕ)) + 2 V_t^ (trans) ( ω_2^ (z)(ϕ) ) 
       +    2 V_t^ (cis) ( ω_3^ (z)(ϕ)) + 
    2 V_t^ (trans) (ω_4^ (z)(ϕ)).






§ FITTING PROCEDURE



The total torsional energy data, (ϕ_i, E_t,i^ (a)) and (ϕ_i, E_t,i^ (z)),
shown in red and black diamonds in the Fig. <ref>, and the analytical (to be fitted) 
expressions U_ t^ (a), U_ t^ (z) given in the Eqs. (<ref>) 
and (<ref>) for the Model 1 (or the Eqs. (<ref>) and (<ref>) for the Model 2)
can be used to obtain the optimal parameters V_1 and V_2 (or K_ cis and K_ trans) of the 
individual torsional terms so that U_ t^ (a), U_ t^ (z) reproduce the 
dependence of E_t^ (a)(ϕ) and E_t^ (z)(ϕ) as close as possible. 
For this purpose, adopting a standard procedure, we minimize an objective function 
O(V_1,V_2) which is the equal-weighted sum of the square differences,

    O(V_1,V_2)    =   ∑_i=1^ϕ_i < ϕ_ max[ E_t,i^ (a) - U_ t^ (a)(ϕ_i)]^2 
       +   ∑_i=1^ϕ_i < ϕ_ max[E_t,i^ (z) - U_ t^ (z)(ϕ_i)]^2  .

The sums in the above expression runs over all i for which ϕ_i is smaller than an 
upper-limit angle ϕ_ max.

For the Model 1, the fitted total torsional energies 
U_ t^ (a)(ϕ) and U_ t^ (z)(ϕ) given in the Eqs. (<ref>) and (<ref>)
depend on V_1 and V_2 through the dependence of the individual
terms V_t of the Eq. (<ref>). 

In the case of Model 2, we are optimizing K_ cis and K_ trans parameters, and the expressions 
(<ref>) and (<ref>) are used instead and the individual terms 
V_t are given by the Eq. (<ref>).

The choice of ϕ_ max is expected to affect the quality of fitting for small and large ϕ.
We are interested in seeing whether the fitting parameters depend on ϕ_ max and, if so,
at what extend.




 §.§ Model 1: fitting Results for V_1, V_2






We performed fitting of V_1, V_2 of Eq. (<ref>), for three different 
ϕ_ max values: 10^∘, 20^∘ and 30^∘.
The optimal parameters V_1 and V_2 given in table <ref>.



In Fig. <ref>, we show, the total torsional energies, E_t^ (a) and 
E_t^ (z) and the fitted lines for these three values of ϕ_ max. 
As we see, the fitting is better for small values of ϕ and deteriorates as ϕ increases.
In all cases, it reproduces the armchair data in closer agreement than
the zig-zag, i.e. for a larger range of ϕ.
For ϕ_ max=10^∘, there is a satisfactory agreement for values of 
ϕ up to 0.55 rad for the armchair case and 0.4 rad for the zig-zag case. For  
ϕ_ max=20^∘, the range of satisfactory agreement increases roughly up to 0.65 rad 
and 0.45 rad for the armchair and zig-zag cases respectively. Finally for 
ϕ_ max=30^∘, 
the agreement range increases further up to 0.7 rad and 0.5 rad, respectively. Although by 
increasing ϕ_ max, the range of satisfactory agreement also increases, this is at the
cost of the agreement for smaller angles. As the code is trying to fit better at larger values
of ϕ the quality for smaller angles deteriorates. This deterioration, however, is rather 
small as we observe in the insets of the Figs. <ref>(a), (b), (c) where we zoom in that 
region. 



However, in general, for angles ϕ up to 0.4 rad (22^∘) corresponding to energies  0.2 to  0.3 eV,
all fittings are satisfactory. On the other hand, as we see in Table <ref>,
the fitted values of V_1 and V_2 are not so sensitive to the
value of ϕ_ max: V_1 remains close to 0 while V_2 is in the range 0.22-0.23 eV. 
In addition, the large value of V_2, i.e. 0.23 eV, obtained for ϕ_ max=20^∘,30^∘ 
performs better for larger angles, up to 0.5 rad (∼30^∘), corresponding to energies of 
0.4-0.5 eV, while on the other hand the fitted results in the region of small ϕ remain 
satisfactory. These considerations suggest that it is quite reasonable to adopt as optimal V_1=0 and V_2=0.23 eV
and our proposed torsional potential has the simple form

    V_t(ω)=1/2V_2(1-cos(2ω)),   V_2=0.23  eV .

This potential is shown in Fig. <ref> together with E_t^ (a) and E_t^ (z).






 §.§ Model 2: fitting results for K_ cis, K_ trans






As in the previous section, we performed the fitting of K_ cis and K_ trans
of Eq. (<ref>) for the same values of ϕ_ max, i.e. 10^∘, 20^∘ and 30^∘.  
The optimal parameters K_ cis and K_ trans for each of these cases are given in 
Table <ref>.
In Fig. <ref>, we show the fitting lines for the armchair and zig-zag folding cases 
for all three values of ϕ_ max
compared to the data points for the total torsional energy per unit cell, E_t^ (a) and 
E_t^ (z). The fitting quality is quite similar to that of the previous section. Again,
 for ϕ_ max=20^∘ and 30^∘, the quality improves for larger 
values of ϕ and at the same time the fit for smaller ϕ does not deteriorate substantially.
Thus, we propose a rounded optimal set 
K_ cis=0.14 eV and K_ trans=0.10 eV which is close to the values obtained for 
ϕ_ max=20^∘ and 30^∘. In Fig. <ref>, we show the torsional energy 
obtained with Model 2 and these values for K_ cis and K_ trans compared with the
data E_t^ (a) and E_t^ (z).





To summarize, for the Model 2, we propose

    [    V_t^ (cis)(ω) = K_ cis[ 1 - cos (2 ω ) ],    K_ cis=0.14  ,; V_t^ (trans)(ω) = K_ trans[ 1 - cos (2 ω ) ],  K_ trans=0.1  , ]
 
where either the first or the second expression is used depending on whether the
torsional angle ω is cis or trans.





 §.§ Comparison of the two models



In Fig. <ref> we show the torsional energies per unit cell, for both fitting forms of Models 1 and 2,
for the  case of the optimal parameters we arrived at.  We notice that the two models
are of the same quality. They almost coincide for the armchair case, while for the
zig-zag, the Model 2 is slightly better for larger ϕ's and the Model 1 marginally better for smaller.
The differences however are not significant for ϕ's up to 0.5 rad. 





Note that the obtained optimal parameters of Model 2, K_ cis=0.14 eV and K_ trans=0.10 eV,
are close to each other, indicating that a single parameter with value the average of them 
would offer a reasonable description. Moreover, this average value is almost equal to V_2 /2.
Thus, it is rather unnecessary to assume different parameters for cis and trans dihedral angles 
and, to keep things as simple as possible, the simple form of the Eq. (<ref>), is quite sufficient to describe
all torsional distortions. Therefore, the Model 1 of Eq. (<ref>) is our proposed one, containing a
single parameter V_2. We should mention that our modeling describes accurately the energy of torsional
angles ϕ up to 0.5 rad which is already a sufficiently large value, corresponding to rather unphysical structural 
deformations.





§ APPLICATION TO THE VERTICAL DISPLACEMENT OF A CARBON ATOM IN GRAPHENE




In order to test the accuracy of the proposed parameters for the torsional terms, 
we consider the deformation energy of graphene due to a vertical, out-of-plane displacement 
of a single carbon atom. We consider that apart from the vertically displaced atom, all other atoms remain
fixed at their equilibrium positions within graphene's plane. The task is to compare 
the deformation energy obtained by the present potential, along with the in-plane force field of Ref. <cit.>,
with that obtained by DFT calculations (using the same method that was used to produce the data discussed above). 

The process of moving a carbon atom vertically outside graphene's plane is described by a deformation 
energy consisting of all kinds of individual terms, i.e. bond-stretching, 
angle-bending and of course torsional terms. 



Concerning the bond-stretching terms, the vertical movement of a carbon atom at a displacement z 
over the plane, alters only the three bonds of that atom (see Fig. <ref>). If the length 
of these bonds at z=0 is d their altered length d' becomes 

    d'(z)=√(d^2+z^2)


At the same time, two different kinds of angle-bending terms appear corresponding to: (i) the three angles, θ_1,
between the atom's bonds (marked in red in Fig. <ref>) and (ii) the six angles, θ_2 between these bonds
and the bonds marked in blue in Fig. <ref>. These angles can be expressed in terms of the displacement z as

    θ_1(z)=arccos(2z^2-d^2/2(d^2+z^2))

and 

    θ_2(z)=arccos(-d/2√(d^2+z^2))


Finally, several torsional terms also contribute. There are rotations around the 3 bonds 
of the displaced atom and its first neighbors (marked in red in Fig. <ref>) as well as 
rotations around the 6 bonds of the first neighbors and the second neighbors. These
rotations correspond to the following torsional angles


  * 6 cis dihedral angles around the bonds marked in red in Fig. <ref> given by
      
    ω_ cis^(1)(z)=arccos(3/4/√((z/d)^2+
              3/4)√(3(z/d)^2+3/4))
 

  * 6 trans dihedral angles around the bonds marked in red in Fig. <ref> given by
      
    ω_ trans^(1)(z)=arccos(-3/2(z/d)^2-
              3/4/√((z/d)^2+3/4)√(3(z/d)^2+3/4))


  * 6 cis dihedral angles around the bonds marked in blue in Fig. <ref> given by
      
    ω^(2)_ cis(z)=arccos(1/√(4/3(z/d)^2+1))


  * 6 trans dihedral angles around the bonds marked in blue in Fig. <ref> given by
       
    ω^(2)_ trans(z)=arccos(-1/√(4/3(z/d)^2+1))






Then, the deformation energy is given by

    E_d(z)    =    3 V_s (d'(z)) + 3 V_b (θ_1(z)) + 6 V_b (θ_2(z))
       +    6 V_t(ω_ cis^(1)(z)) + 6 V_t(ω_ trans^(1)(z)) 
       +    6 V_t(ω_ cis^(2)(z)) + 6 V_t(ω_ trans^(2)(z)) ,

where V_s (d), V_b (θ) are respectively the bond stretching and angle bending terms 
given in Ref. <cit.>. V_t is the individual torsional term of Eq. (<ref>)
for Model 1, while for Model 2 it should be replaced by V_t^ (cis) or V_t^ (trans)
of Eq. (<ref>) for the two ω_ cis and the two ω_ trans
respectively.
In Eq. (<ref>), the deformation energy, E_d, becomes an analytic function of z through the
explicit dependence on z of the bonds d', the angles θ_1, θ_2 and the dihedral angles
ω_ cis^(1), ω_ trans^(1), ω_ cis^(2), ω_ trans^(2).
 
In Fig. <ref>,  we show E_d in comparison with the corresponding DFT results. 
As we see, both models perform equally well and the error does not exceed 0.05 eV for the 
deformation range shown.  There are small differences in their 
agreement with DFT, for example Model 2 seems slightly better for indermediate-size 
displacements while Model 1 for larger ones, however, these small differences are 
insignificant validating our preference for Model 1 on the basis of its simplicity. 



§ CONCLUSION


In summary, we present a simple torsional force field for graphene and other sp^2 carbon nanostructures.
To obtain this potential we performed two sets of DFT calculations by folding two different graphene nanoribbon
 structures around their middle line. The first set of calculations concern the folding of an armchair nanoribbon 
around its middle line, an armchair direction. The second concerns the folding of a zig-zag nanoribbon around its 
middle line which is a zig-zag direction. From the deformation energies we isolated the “pure” torsional
contribution by removing angle bending terms with the use of our previously proposed angle bending terms<cit.>.
The purified torsional deformation energy was then fit two different analytic forms with two parameters each that we
call Models 1 and 2. The first (Model 1) is that of Eq. (<ref>) and does not distinguish torsional angles, while
the second (Model 2) of Eq. (<ref>) treats differently “cis” and “trans” torsional angles. We found that 
the form of Model 1 reduces to one parameter form, see Eq. (<ref>), which was found to be an average of the “cis” 
and “trans” terms of Model 2 which differ very little from each other, see Eq. (<ref>). That suggests that the 
use of two different terms is redundant and the single term of Model 1 suffices at a reasonable level or accuracy. We found 
that both models reproduce accurately the torsional deformation energy of graphene nanoribbons due to the folding we 
considered up to θ_ max of the order of 30^o (≈ 0.5 rad).

As an additional validation test we considered the case of the deformation energy due to the vertical displacement of a 
single C atom outside of graphene's plane. For this task we used the torsional terms of either Model 1 or 2 combined with 
the bond stretching and angle bending terms of Ref. kalosakas. For all terms in this formula we provide 
analytic expressions in terms of the displacement z. We found that both models perform equally well in this case with errors
not exceeding 0.05 eV for a relatively large range of z, up to 0.4-0.5 Å. The good performance of both Models in this case
validates our choice for the simpler Model 1. 

The torsional force field presented here, in combination with the bond stretching and angle bending terms of Ref. kalosakas,
which were also fitted to DFT results using the same density functional approximation provide a complete, accurate, but simple in form
atomistic potential, which is computationally efficient due to its simplicity. Thus, we expect that 
it will be proven a very useful tool for large scale atomistic simulation of graphene and other sp^2 nanostructures.   

Acknowledgements:
We acknowledge helpful discussions with E. N. Koukaras.
The research leading to the present results has received funding from Thales project
“GRAPHENECOMP”, co-financed by the European Union (ESF) and the Greek Ministry of Education 
(through ΕΣΠΑ program). NNL  acknowledges support from the Hellenic Ministry of Education/GSRT (ESPA), through 
"Advanced Materials and Devices" program (MIS:5002409) and EU H2020 ETN project ‘Enabling Excellence’ Grant 
Agreement 642742.
 
10

Novoselov666  K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva,
A. A. Firsov, Science 306, 666 (2004).

castroneto A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).

synthesis  L. Chen, Y. Hernandez, X. Feng, K. Mullen, Angew. Chem. Int. Ed. 51, 7640 (2012).

roadmap  A. C. Ferrari et al, Nanoscale 7, 4598 (2015).

elasticRev  C. Daniels, A. Horning, A. Phillips, D. V. P. Massote, L. Liang, Z. Bullard, B. G. Sumpter, V. Meunier, J. Phys.: Condens. Matter 27, 373002 (2015).

avouris  Y.-M. Lin, et al, Science 332, 1294 (2011).

Kusmartsev2015  F. V. Kusmartsev, W. M. Wu, M. P. Pierpoint, K. C. Yung, in Applied Spectroscopy and the Science of Nanomaterials, Ed. P. Misra, Springer Singapore 2015.

Celebi289  K. Celebi, et al, Science 344, 289 (2014).

batteries  M. F. El-Kady, Y. Shao, R. B. Kaner, Nature Rev. Mater. 1, 16033 (2016).

Lee385   C. Lee, X. Wei, J.W. Kysar, J. Hone, Science  321, 385 (2008).

SMLL  G. Tsoukleri, J. Parthenios, K. Papagelis, R. Jalil, A.C. Ferrari, A.K. Geim, K.S. Novoselov, C. Galiotis, Small 21, 2397 (2009).

kalosakas G. Kalosakas, N. N. Lathiotakis, C. Galiotis, K. Papagelis, J. App. Phys. 113, 134307 (2013).

fasolinoNM  A. Fasolino, J. H. Los, M. I. Katsnelson, Nature Mater. 6, 858 (2007).

fasolino  K. V. Zakharchenko, M. I. Katsnelson, A. Fasolino, Phys. Rev. Lett. 102, 046808 (2009).

bending P. Liu, Y. W. Zhang, Appl. Phys. Lett. 94, 231912 (2009)

buechler  Z. Xu, M. J. Buechler, ACS Nano 4, 3869 (2010).

peeters1  M. Neek-Amal, F. M. Peeters, Phys. Rev. B 82, 085432 (2010).

peeters2  M. Neek-Amal, F. M. Peeters, Appl. Phys. Lett. 97, 153118 (2010).

indent  X. Tan, J. Wu, K. Zhang, X. Peng, L. Sun, J. Zhong, Appl. Phys. Lett. 102, 071908 (2013).

campbell  Z. Qi, D. A. Bahamon, V. M. Pereira, H. S. Park, D. K. Campbell, A. H. Castro Neto, Nano Lett. 13, 2692 (2013).

arisRSC  A. P. Sgouros, G. Kalosakas, M. M. Sigalas, K. Papagelis, RSC Adv. 5, 39930 (2015).

koukaras E. N. Koukaras, G. Kalosakas, C. Galiotis, K. Papagelis, Sci. Rep. 5, 12923 (2015).

aris2d A. P. Sgouros, G. Kalosakas, C. Galiotis, K. Papagelis, 2D Mater. 3, 025033 (2016).

Tersoff J. Tersoff, Phys. Rev. Lett. 61, 2879 (1988).

Tersoff1 J. Tersoff, Phys. Rev. B 37, 6991 (1988).

Brenner D. W. Brenner, Phys. Rev. B 42, 9458 (1990).

T2010   L. Lindsay, D. A. Broido, Phys. Rev. B 81, 205441 (2010).

LCBOP  J. H. Los, A. Fasolino, Phys. Rev. B 68, 024107 (2003).

LCBOPII  J. H. Los, L. M. Ghiringhelli, E. J. Meijer, A. Fasolino, Phys. Rev. B 72, 214102 (2005).

AIREBO   S. J. Stuart. A. B. Tutein, J. A. Harrison, J. Chem. Phys. 112, 6472 (2000).

Wei  D. Wei, Y. Song, and F. Wang, J. Chem. Phys. 134, 184704 (2011).

allo   Z. G. Fthenakis, G. Kalosakas, G. D. Chatzidakis, C. Galiotis, K. Papagelis, N. N. Lathiotakis, preprint (2017). 

QE-2009   P. Giannozzi et al, J. Phys.: Cond. Matt. 21, 395502 (2009).

pseudo   A. Dal Corso, http://www.quantum-espresso.org/wp-content/uploads/upf_files/C.pbe-rrkjus.UPF




  