Atomic and Electronic Structure of 
Carbon Nanotubes

written by 
Mads Brandbyge, mbr@mic.dtu.dk
Mikroelectronik Centret (MIC)
Danish Technical University,  Build. 345 east
DK-2800 Lyngby, DENMARK

Misc.

Turn off unneeded error messages first

Off[General::spell1, Graphics::realu]

Lines connecting the nearest neighbour atoms:

\(Clear[makelines]; \ 

    makelines[l_, range_] := 
      Module[{lines, n, i, j, dist}, 
		lines = {}; 
		n = Length[l]; 

        		For[i = 1, i <= n, \(i++\), 
			
          \(For[j = 1, j < i, \(j++\), 
				
            dist = Sqrt[
                \((l[\([i]\)] - l[\([j]\)])\) . 
                  \((l[\([i]\)] - l[\([j]\)])\)]; 
		
            If[dist\  < \ range, \ \ 	
              \(lines = 
                Append[lines, \ \ \ 
                  Line[{\ \ l[\([i]\)], l[\([j]\)]\ }]\ \ \ ]; \)

              					]; 					
			]; \)
		]; 
		
        Return[lines]]\)

Introduction

Carbon nanotubes are huge tube-like molecules consisting only of Carbon atoms with a diameter in the nanometer range and lengths up to microns. 
The atoms are sitting in a hexagonal arrangement with 3-fold coordination (3 nearest neighbours). The nanotube can be viewed as a rolled up sheet of
Graphene i.e. a 2 dimensional hexagonal lattice of Carbon atoms.

The interatomic distance in Å:

d = 1.42;

Unit vectors of the hexagonal lattice (like in a sheet of graphite = graphene) lying in the xy-plane:

e1=(2*d*Cos[Pi/6])*{Cos[-Pi/6],Sin[-Pi/6],0};

e2=(2*d*Cos[Pi/6])*{Cos[Pi/6],Sin[Pi/6],0};

The atoms in the basis of the hexagonal lattice are often called A-atoms and B-atoms. A-atoms are equivalent with the
atom at {0,0,0} and B-atoms are equivalent with the atom at {d,0,0}.

\(\(Acoor = 
      Flatten[Table[n*e1\  + \ m*e2, {n, 0, 5}, 
{m, 0, 5}], 1] // N; \)\)

Bcoor = {d, 0, 0} + #&/@Acoor;

grapheneplot=
Show[Graphics3D[{
{RGBColor[0,0,1],AbsolutePointSize[20],
Point[#]}&/@Acoor,
{RGBColor[1,0,0],AbsolutePointSize[20],
Point[#]}&/@Bcoor,
{AbsoluteThickness[6],#}&/@(makelines[Flatten[{Acoor,Bcoor},1],1.1*d])
}]]

The wrap up of the nanotube is done by identifying the A-atom at {0,0,0} with the A-atom at nn1*e1+nn2*e2, 
where nn1 and nn2 are integer numbers: These are indices which describe the nanotube geometry completely.
The vector  nn1*e1+nn2*e2 is the circumferential vector and the tube will be directed perpendicular to this. 
The structure formed will in general have chiral (translation along the tube axis + rotation around the tube axis) 
and rotational symmetry. These symmetries are used to classify the quantum energy eigenstates of the tubes below.

The Carbon nanotubes are e.g. described in  the paper by White & Mintmire, Carbon,Vol. 33, No. 7, pp. 893-902, 1995. In the notebook the notation will primarily be adapted from this paper. A nice introduction by C. Dekker can be found in Physics Today, May 1999 p. 22.
The homepage of C. Dekker (http://vortex.tn.tudelft.nl/~dekker/) contains nice pictures from this latter reference.

Atomic Structure

nn1 and nn2 is the index defining the carbon nanotube (defining its chiral vector), 
u is a measure of the length of the piece we want to plot.

\(Clear[makenanotube]; 

    makenanotube[nn1_, nn2_, u_] := 
      Module[
		{coors, i, j, theta, rotxy, E1, E2, NC, H, B, BE, h, phi, 
          Radius, 
			TA, Btmp, BEtmp, S, dbasis, hbasis, phibasis}, 

        		
		n1 = Min[nn1, nn2]; 
		n2 = Max[nn1, nn2]; 

	
        Btmp = \((n1*e1 + n2*e2)\); 
		BEtmp = Btmp/Sqrt[Btmp . Btmp]; 

        		 (*\ redefine\ e1, \ e2\ so\ B, BE\ is\ along\ the\ x - axis, \ 
          \(--\ the\)\ Tube\ is\ along\ y\ *) 
		
        theta = ArcCos[BEtmp . {1, 0, 0}]; 
		
        rotxy = {{Cos[theta], Sin[theta], 0}, {\(-Sin[theta]\), Cos[theta], 
              0}, {0, 0, 1}}; 
		E1 = rotxy . e1; 
		E2 = rotxy . e2; 

        		B = {Sqrt[Btmp . Btmp], 0, 0}; 
		BE = {1, 0, 0}; 
		

        			Radius = Sqrt[B . B]/\((2  Pi)\); 
	\ \ 
        Print["Radius = ", Radius]; 
		

        		 (*\ Determine\ m1, m2\ *) 
				NC = GCD[n1, n2]; 

        \ \ \ \ \ \ \ \ \ 
        Print["Rotational symmetry with ", NC, "-order axis"]; 

        				m1 = 0.5; 
		m2 = \(-1\); \ 
	
        While[IntegerPart[m1] < m1, 
			m2 = m2 + 1; 
			
          m1 = \((1 + m2*\((n1/NC)\))\)/\((n2/NC)\); 
		]; 
		
        Print["m1,m2 = ", m1, ",", m2]; 
		H = E1*m1 + E2*m2; 

        		h = Sqrt[\((Cross[H, B] . Cross[H, B])\)/\((B . B)\)]; 
		
        phi = 2  Pi \((H . B)\)/\((B . B)\); 

					
        Print["(h,phi) = (", h, ",", phi, ")"]; 


        		 (*\ Chiral\ 
            \(operator : \ 
              translate\ h\ in\ y\ dir . \ and\ rotate\ phi\ around\ y\ axis
                \)\ *) 
	
        S[r_] := 
          \ \ {0, h, 0}\  + 
            \((\ {{Cos[phi], 0, Sin[phi]}, {0, 1, 0}, {\(-Sin[phi]\), 0, 
                    Cos[phi]}} . r)\); \ 
		

        		 (*\ Rotation\ operator\ Cn\ *) 
		
        Cn[r_] = 
          \ {{Cos[2  Pi/NC], 0, Sin[2  Pi/NC]}, {0, 1, 0}, {
                \(-Sin[2  Pi/NC]\), 0, Cos[2  Pi/NC]}} . r; 
		

        		 (*\ unit - cell\ coordinates\ uc\ *) 
	
		
        aatom = {Radius, 0, 0}; \ 
		dbasis = \((E1 + E2)\)/3\ ; 
		
        hbasis = 
          Sqrt[\((Cross[dbasis, {1, 0, 0}] . Cross[dbasis, {1, 0, 0}])\)]; 

        		phibasis = 2  Pi*\((dbasis . B)\)/\((B . B)\); 
		
        batom = 	\ {0, hbasis, 0}\  + 
					
            \((\ {{Cos[phibasis], 0, Sin[phibasis]}, {0, 1, 0}, {
                    \(-Sin[phibasis]\), 0, Cos[phibasis]}} . {Radius, 0, 0})
              \); \ 
		

        	 (*\ Generate\ coordinates\ using\ symmetry\ operations\ *) 
		

        		coors = 
          Flatten[
				
            Table[{Nest[Cn, Nest[S, aatom, i], j], 
                Nest[Cn, Nest[S, batom, i], j]}, {i, 0, u}, {j, 1, NC}], 1]; 
        
		Return[coors]

	]\)

Plot it and make a small movie

Note that the distances are not quite equal to the d value
-- especially for small, highly curved tubes

ntcoor = makenanotube[5, 5, 10];

nanotubeplot=
Show[Graphics3D[{
{RGBColor[0,0,1],AbsolutePointSize[20],
Point[#]}&/@Transpose[ntcoor][[1]],
{RGBColor[1,0,0],AbsolutePointSize[20],
Point[#]}&/@Transpose[ntcoor][[2]],

{AbsoluteThickness[6],#}&/@(makelines[Flatten[ntcoor,1],1.1*d])
}],
ViewPoint-> {-2.136, -0.058, 0.112} ]

cameraat = {ViewPoint -> {-4.999,  0.048,  0.062}, 
        ViewPoint -> {-4.786,  -0.146,  1.441}, 
        ViewPoint -> {-3.971,  -0.121,  3.035}, 
        ViewPoint -> {-2.733,  0.026,  4.187}, 
        ViewPoint -> {-0.986,  -0.010,  4.902}, 
        ViewPoint -> {-0.392,  -0.004,  4.985}, 
        ViewPoint -> {-0.353,  -0.004,  4.988}, 
        ViewPoint -> {-0.329,  0.127,  4.988}, 
        ViewPoint -> {-0.286,  0.208,  4.988}, 
        ViewPoint -> {-0.214,  0.281,  4.988}, 
        ViewPoint -> {-0.135,  0.326,  4.988}, 
        ViewPoint -> {-0.022,  0.352,  4.988}, 
        ViewPoint -> {0.093,  0.341,  4.988}, 
        ViewPoint -> {0.203,  0.289,  4.988}, 
        ViewPoint -> {0.292,  0.198,  4.988}, 
        ViewPoint -> {0.346,  0.069,  4.988}, 
        ViewPoint -> {0.353,  0.014,  4.988}, 
        ViewPoint -> {1.529,  -0.001,  4.761}, 
        ViewPoint -> {3.108,  -0.002,  3.917}, 
        ViewPoint -> {4.597,  -0.003,  1.968}, 
        ViewPoint -> {4.871,  -0.004,  1.129}, 
        ViewPoint -> {4.998,  -0.004,  -0.137}, 
        ViewPoint -> {4.859,  -0.692,  -0.135}, 
        ViewPoint -> {4.338,  -2.297,  -0.135}, 
        ViewPoint -> {3.732,  -3.188,  -0.135}, 
        ViewPoint -> {3.039,  -3.854,  -0.135}, 
        ViewPoint -> {2.331,  -4.319,  -0.135}, 
        ViewPoint -> {1.295,  -4.734,  -0.135}, 
        ViewPoint -> {-0.116,  -4.907,  -0.135}, 
        ViewPoint -> {-0.100,  -4.257,  -0.117}, 
        ViewPoint -> {-0.093,  -3.957,  -0.109}, 
        ViewPoint -> {-0.071,  -3.018,  -0.083}, 
        ViewPoint -> {-0.059,  -2.518,  -0.069}, 
        ViewPoint -> {-0.055,  -2.318,  -0.064}, 
        ViewPoint -> {-0.046,  -1.949,  -0.054}, 
        ViewPoint -> {-0.039,  -1.669,  -0.046}, 
        ViewPoint -> {-0.036,  -1.519,  -0.042}, 
        ViewPoint -> {-0.028,  -1.189,  -0.033}, 
        ViewPoint -> {-0.024,  -1.009,  -0.028}, 
        ViewPoint -> {-0.026,  -1.089,  -0.030}, 
        ViewPoint -> {-0.027,  -1.149,  -0.032}, 
        ViewPoint -> {-0.029,  -1.249,  -0.034}, 
        ViewPoint -> {-0.033,  -1.389,  -0.038}, 
        ViewPoint -> {-0.037,  -1.569,  -0.043}, 
        ViewPoint -> {-0.041,  -1.759,  -0.048}, 
        ViewPoint -> {-0.053,  -2.229,  -0.061}, 
        ViewPoint -> {-0.068,  -2.878,  -0.079}, 
        ViewPoint -> {-0.088,  -3.738,  -0.103}, 
        ViewPoint -> {-0.122,  -5.187,  -0.143}, 
        ViewPoint -> {-1.052,  -5.080,  -0.143}, 
        ViewPoint -> {-2.848,  -4.336,  -0.143}, 
        ViewPoint -> {-4.336,  -2.848,  -0.143}, 
        ViewPoint -> {-4.971,  -1.487,  -0.143}, 
        ViewPoint -> {-5.130,  -0.771,  -0.143}, 
        ViewPoint -> {-5.188,  0.000,  -0.143}};

ntmovie = 
      Table[Show[nanotubeplot, cameraat[[i]]], {i, Length[cameraat]}]

Electronic Structure (only π-bands)

The electronic quantum energy eigenstates of the nanotubes can be classified/indexed by a rotational quantum number 
corresponding to the rotational symmetry of the nanotube and a chiral quantum number corresponding to the 
chiral symmetry. A nanotube will of cause also have translational symmetry but this symmetry is already contained in the (more complete) chiral symmetry. Below we will plot the eigenenergy versus chiral-symmetry index (κ). We will also make the corresponding plot where we use the translational symmetry index - or k-point which is the 1-dimensional band-structure known from basic solid state physics textbooks.

Hamiltonian (in eV):

Vpppi = -2.77;

ham[r1_, r2_] := 
      Vpppi*If[(r1 - r2) . (r1 - r2)  <  (1.1*d)^2, 1.0, 0.0];

\(\(Clear[nanotubebands]; 

    nanotubebands[nn1_, nn2_] := 
      Module[
		{othercoors, i, j, theta, rotxy, E1, E2, NC, H, B, BE, h, 
          phi, Radius, aatom, batom, m1, m2, n1, n2, 
			Btmp, BEtmp, S, 
          invS, Cn, dbasis, hbasis, phibasis, phaselist, 
		bandplotkappa, 
          bandplotky, allbandplots, HAB, M, kappa, NS, eig1, k1, kybands, T, 
          ky}, 
		
		n1 = Min[nn1, nn2]; 
		n2 = Max[nn1, nn2]; 


        	\ \ Btmp = \((n1*e1 + n2*e2)\); 
		
        BEtmp = Btmp/Sqrt[Btmp . Btmp]; 
		

        		 (*\ redefine\ e1, \ e2\ so\ B, BE\ is\ along\ the\ x - axis, \ 
          \(--\ the\)\ Tube\ is\ along\ y\ *) 
		
        theta = ArcCos[BEtmp . {1, 0, 0}]; 
		
        rotxy = {{Cos[theta], Sin[theta], 0}, {\(-Sin[theta]\), Cos[theta], 
              0}, {0, 0, 1}}; 
		E1 = rotxy . e1; 
		E2 = rotxy . e2; 

        		B = {Sqrt[Btmp . Btmp], 0, 0}; 
		BE = {1, 0, 0}; 
		

        			Radius = Sqrt[B . B]/\((2  Pi)\); 
	\ \ 
        Print["Radius = ", Radius]; 
		

        		 (*\ Determine\ m1, m2\ *) 
				NC = GCD[n1, n2]; 

        \ \ \ \ \ \ \ \ \ 
        Print["Rotational symmetry with ", NC, "-order axis"]; 

        				m1 = 0.5; 
		m2 = \(-1\); \ 
	
        While[IntegerPart[m1] < m1, 
			m2 = m2 + 1; 
			
          m1 = \((1 + m2*\((n1/NC)\))\)/\((n2/NC)\); 
		]; 
		
        Print["m1,m2 = ", m1, ",", m2]; 
		H = E1*m1 + E2*m2; 

        		h = Sqrt[\((Cross[H, B] . Cross[H, B])\)/\((B . B)\)]; 
		
        phi = 2  Pi \((H . B)\)/\((B . B)\); 
		
        Print["(h,phi) = (", h, ",", phi, ")"]; 


        		 (*\ Chiral\ 
            \(operator : \ 
              translate\ h\ in\ y\ dir . \ and\ rotate\ phi\ around\ y\ axis
                \)\ *) 
				Clear[S]; 
	
        S[r_] := 
          \ \ {0, h, 0}\  + 
            \((\ {{Cos[phi], 0, Sin[phi]}, {0, 1, 0}, {\(-Sin[phi]\), 0, 
                    Cos[phi]}} . r)\); \ 
			
				Clear[invS]; 
		
        invS[r_] := 
          \ \ {0, \(-h\), 0}\  + 
            \((\ {{Cos[\(-phi\)], 0, Sin[\(-phi\)]}, {0, 1, 0}, {
                    \(-Sin[\(-phi\)]\), 0, Cos[\(-phi\)]}} . r)\); 
		

        		 (*\ Rotation\ operator\ Cn\ *) 
		Clear[Cn]; 
		
        Cn[r_] := 
          \ {{Cos[2  Pi/NC], 0, Sin[2  Pi/NC]}, {0, 1, 0}, {
                \(-Sin[2  Pi/NC]\), 0, Cos[2  Pi/NC]}} . r; 
		

        		 (*\ unit - cell\ coordinates\ uc\ *) 
		
        aatom = {Radius, 0, 0}; \ 
		dbasis = \((E1 + E2)\)/3\ ; 
		
        hbasis = 
          Sqrt[\((Cross[dbasis, {1, 0, 0}] . Cross[dbasis, {1, 0, 0}])\)]; 

        		phibasis = 2  Pi*\((dbasis . B)\)/\((B . B)\); 
		
        batom = 	\ {0, hbasis, 0}\  + 
					
            \((\ {{Cos[phibasis], 0, Sin[phibasis]}, {0, 1, 0}, {
                    \(-Sin[phibasis]\), 0, Cos[phibasis]}} . {Radius, 0, 0})
              \); \ 
	Print["phibasis = ", phibasis]; 
	

        		 (*\ nanotube\ y\ unit - cell\ *) 
		
        dR = GCD[2*\((n1)\) + n2, 2*\((n2)\) + n1]; 
		
        T = Sqrt[3]*2  Pi*Radius/dR; 
		
        Print["Translational symmetry along y: T = ", T]; 	
		
        NT = Ceiling[T/h]; 
			
        Print["Corresponding to ", NT, \ " chiral operations"]; 	

        	 (*\ Generate\ bands\ from\ coordinates\ using\ symmetry\ operations
            \ *) 
		

        		 (*\ Note\ we\ use\ that\ we\ have\ chiral\ "inversion"\ 
            symmetry\ in\ the\ sense\ *) 

        		 (*\ that\ \((1/S)\)\ generates\ an\ equivalent\ atom
            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) \ 
		
        
		 (*\ 
          We\ must\ find\ matrix\ elements\ which\ give\ non - 
            zero\ value\ on\ the\ B\ basis\ \(state : \)\ *) 

        		 (*\ We\ search\ among\ NS\ A - 
            equivalent\ coordinates\ \((excl . \ A\ it\ self)\)\ *) 
		
        NS = 2*\((n1 + n2)\); 
		
        abhamlist = 
          Flatten[

            				{
              Table[ham[aatom, Nest[Cn, Evaluate[Nest[invS, batom, i]], j]], {
                  j, 0, NC - 1}, {i, 1, NS}], 
\ \ \ \ \ \ \ \ \ \ 
              Table[ham[aatom, Nest[Cn, Evaluate[Nest[\ S, batom, i]], j]], {
                  j, 0, NC - 1}, {i, 0, NS}]}
			\ \ \ \ \ \ \ \ \ ]; 
	
        
		 (*\ 
          Table\ of\ bandstructure\ plots\ for\ different\ rotation\ 
            \(symmetry : \ M\)\ *) 
		Clear[kappa]; 
		
        allbandplots = 
          Table[	
            Print[" ======= Bands for rotational quantum number M = ", M, 
              " ======= "]; \ 
				
            phaselist[kappa_] = 
					
              Flatten[{
                  Table[Exp[\(+I\)*kappa*i]*Exp[\(-2\) Pi*I*M*j/NC], {j, 0, 
                      NC - 1}, {i, 1, NS}], 
                  Table[Exp[\(-I\)*kappa*i]*Exp[\(-2\) Pi*I*M*j/NC], {j, 0, 
                      NC - 1}, {i, 0, NS}]}]; 
				
            HAB[kappa_] = 
              \((Flatten[phaselist[kappa]])\) . \((Flatten[abhamlist])\); 

            				bandplotkappa = 
              Plot[{\(-Abs[HAB[kappa]]\), Abs[HAB[kappa]]}, {kappa, \(-Pi\), 
                  Pi}, PlotStyle -> {AbsoluteThickness[0.15]}, Frame -> True, 
                FrameLabel -> {"κ", "ε"}, 

                				\ \ \ 
                FrameTicks -> {{\(-Pi\), \(-Pi\)/2, 0, Pi/2, Pi}, 
                    Automatic, {}, {}}]; 
				

            ky[kappa_] = 
              Mod[\((kappa)\)*T/h\ \  - \ 
                  \((2  Pi/NC)\)*M*\((\((phi*T/h)\)/\((2  Pi/NC)\))\), 
                2  Pi]; 	
				
						
            kybands = 
              Table[					
\ \ \ \ \ \ \ \ \ \ eig1 = Abs[HAB[kappa]]; \ 
                k1 = ky[kappa]; 

                	\ \ \ \ \ \ \ 		{{\(-k1\), \(-eig1\)}, {k1, \(-eig1\)}, {
                    \(-k1\), eig1}, {k1, eig1}}, 

                	\ \ \ \ \ \ \ 		{kappa, \(-Pi\), Pi, NC*2  Pi/400}]; \ 

            				bandplotky = 
              ListPlot[Flatten[kybands, 1], 
                PlotStyle -> {AbsolutePointSize[3]}, \ 
						
                PlotLabel -> M, Frame -> True, 
                FrameLabel -> {\*""\!\(k\_y\)"", 
                    "ε"}, 
                PlotRange -> {{\(-Pi\), Pi}, Automatic}, \ 
					
                FrameTicks -> {{\(-Pi\), \(-Pi\)/2, 0, Pi/2, Pi}, 
                    Automatic, {}, {}}]; 

            					{bandplotkappa, bandplotky}, {M, 0, NC - 1}]; 
		
        Print["All bands together:"]; 
	
        Show[\(Transpose[allbandplots]\)[\([1]\)]]; 
        Show[\(Transpose[allbandplots]\)[\([2]\)]]; 
					]
	\)\)

Example: (5,5) Tube

nanotubebands[5, 5]

Example: (9,2) Tube

nanotubebands[9, 2]

More detailed electronic structure calculation with sp-orbital 
basis using Slater-Koster matrix elements

Hamiltonian

White & Mintmire's (see e.g. Carbon,Vol. 33, No. 7, pp. 893-902, 1995 for a review)  parameters (in eV):

es = -6. ; ep = 0. ; Vpppi = -2.77; Vpps = 4.37; 
    Vsss = -4.76;  Vsps = 4.33;

The order of orbitals is: s, px, py, pz

SKham calculates 
          HAB  :   which is part of the atom A,  atom B 	(({{HAA, HAB}, {HBA, HBB}}))  hamiltonian,
       where A and B refers to the two basis 
            atoms .  Each		 block is (number of orbitals) x 
              (number of orbitals) = \(4\ x\ 4. \ If\ the\ distance\ does\ not\ fit\ a\ typical\ A - 
                B\ distance\ it\ returns\ zero . 
				
HAB\)=(({{Vssσ, l Vspσ, m Vspσ, n Vspσ}, {-l Vspσ, (1 - (l^2)) Vppπ + (l^2) Vppσ, l m (-Vppπ + Vppσ), l n (-Vppπ + Vppσ)}, {-m Vspσ, l m (-Vppπ + Vppσ), (1 - (m^2)) Vppπ + (m^2) Vppσ, m n (-Vppπ + Vppσ)}, {-n Vspσ, l n (-Vppπ + Vppσ), m n (-Vppπ + Vppσ), (1 - (n^2)) Vppπ + (n^2) Vppσ}}))  using px,py,pz
\(However\ we\ intend\ to\ use\ the\ {s, pr, pl = py, pt}\ basis\ which\ 
            is\ invariant\ with\ respect\ to\ the\ symmetry\ 
            \(operations . \)
	 | \(pr > \)\  = \ 
        Cos \((φ)\)\  | px > \ \(+\ \ Sin\) \((φ)\)\  | 
          \(pz > \), \ \ 	
      \(|\( pt > \)\)\  = \ 
        \(-Sin\) \((φ)\)\  | 
          px > \ \(+\ \ Cos\) \((φ)\)\  | \(px > \)\)

\(\(
Clear[SKHam, l, m, n]; 

    SKHam[l_, m_, n_] = 
      	{
				{Vsss, l*Vsps, m*Vsps, n*Vsps}, 

        		\ {\(-l\)*Vsps, l*l*Vpps + \((1 - l*l)\)*Vpppi, 
          l*m*\((Vpps - Vpppi)\), l*n*\((Vpps - Vpppi)\)}, 

        			{\(-m\)*Vsps, l*m*\((Vpps - Vpppi)\), 
          m*m*Vpps + \((1 - m*m)\) Vpppi, m*n*\((Vpps - Vpppi)\)}, 

        			{\(-n\)*Vsps, l*n*\((Vpps - Vpppi)\), m*n*\((Vpps - Vpppi)\), 
          n*n*Vpps + \((1 - n*n)\)*Vpppi}}; \)\)

SKhamOnsite = DiagonalMatrix[{es, ep, ep, ep, es, ep, ep, ep}];

\(\(
Clear[SKhamHOP12]; 

    SKhamHOP12[r1_, r2_] := 
      Module[{R, l, m, n, p, r1xz, r2xz, mtrans1, mtrans2}, 
		
        R = r1 - r2; 
		LR = Sqrt[R . R]; 
		
        If[LR\  > \ \((1.1*d)\), 
          Return[DiagonalMatrix[{0, 0, 0, 0, 0, 0, 0, 0}]]]; 
			
        If[LR\  < \ \((0.9*d)\), 
          Return[DiagonalMatrix[{0, 0, 0, 0, 0, 0, 0, 0}]]]; 
	UR = R/LR; 

        		{l, m, n} = UR; 	
		r1xz = {r1[\([1]\)], r1[\([3]\)]}; 
		
        r1xz = r1xz/Sqrt[r1xz . r1xz]; 
		
        r2xz = {r2[\([1]\)], r2[\([3]\)]}; 
	\ \ 
        r2xz = r2xz/Sqrt[r2xz . r2xz]; 
		
        mtrans1 = {{1, 0, 0, 0}, {0, r1xz[\([1]\)], 0, \(-r1xz[\([2]\)]\)}, {
              0, 0, 1, 0}, {0, r1xz[\([2]\)], 0, r1xz[\([1]\)]}}; 
		
        mtrans2 = {{1, 0, 0, 0}, {0, r2xz[\([1]\)], 0, \(-r2xz[\([2]\)]\)}, {
              0, 0, 1, 0}, {0, r2xz[\([2]\)], 0, r2xz[\([1]\)]}}; 

        			 (*\ the\ order\ in\ the\ matrix\ is\ s, pr, pl = py, pt\ *) 

        		m12 = Transpose[mtrans2] . SKHam[l, m, n] . mtrans1; 
		
        H12 = Transpose[{{1, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0, 0, 0}, {
                  0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0, 0}}] . 
            m12 . {{0, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 
                0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 1}}; 
	

        Return[H12]
		]

\)\)

nkpts is the number of kappa-points (chiral quantum number) we want to sample

\(\(Clear[SKnanotubebands]; 

    SKnanotubebands[nn1_, nn2_, nkpts_] := 
      Module[
		{othercoors, i, j, theta, rotxy, E1, E2, NC, H, B, BE, h, 
          phi, Radius, aatom, batom, m1, m2, n1, n2, 
			Btmp, BEtmp, S, 
          invS, Cn, dbasis, hbasis, phibasis, 
		bands, bandplots, 
          eiglist, HAM, M, kappa, NS, dottmp, abhamlist}, 
		
		
        n1 = Min[nn1, nn2]; 
		n2 = Max[nn1, nn2]; 

	\ \ 
        Btmp = \((n1*e1 + n2*e2)\); 
		BEtmp = Btmp/Sqrt[Btmp . Btmp]; 

        		
		 (*\ redefine\ e1, \ e2\ so\ B, 
          BE\ is\ along\ the\ x - axis, \ \(--\ the\)\ Tube\ is\ along\ y\ *) 
          
		theta = ArcCos[BEtmp . {1, 0, 0}]; 
		
        rotxy = {{Cos[theta], Sin[theta], 0}, {\(-Sin[theta]\), Cos[theta], 
              0}, {0, 0, 1}}; 
		E1 = rotxy . e1; 
		E2 = rotxy . e2; 

        		B = {Sqrt[Btmp . Btmp], 0, 0}; 
		BE = {1, 0, 0}; 
		

        			Radius = Sqrt[B . B]/\((2  Pi)\); 
	\ \ 
        Print["Radius = ", Radius]; 
		

        		 (*\ Determine\ m1, m2\ *) 
				NC = GCD[n1, n2]; 

        \ \ \ \ \ \ \ \ \ 
        Print["Rotational symmetry with ", NC, "-order axis"]; 

        				m1 = 0.5; 
		m2 = \(-1\); \ 
	
        While[IntegerPart[m1] < m1, 
			m2 = m2 + 1; 
			
          m1 = \((1 + m2*\((n1/NC)\))\)/\((n2/NC)\); 
		]; 
		
        Print["m1,m2 = ", m1, ",", m2]; 
		H = E1*m1 + E2*m2; 

        		h = Sqrt[\((Cross[H, B] . Cross[H, B])\)/\((B . B)\)]; 
		
        phi = 2  Pi \((H . B)\)/\((B . B)\); 
		
        Print["(h,phi) = (", h, ",", phi, ")"]; 


        		 (*\ Chiral\ 
            \(operator : \ 
              translate\ h\ in\ y\ dir . \ and\ rotate\ phi\ around\ y\ axis
                \)\ *) 
				Clear[S]; 
	
        S[r_] := 
          \ \ {0, h, 0}\  + 
            \((\ {{Cos[phi], 0, Sin[phi]}, {0, 1, 0}, {\(-Sin[phi]\), 0, 
                    Cos[phi]}} . r)\); \ 
			
				Clear[invS]; 
		
        invS[r_] := 
          \ \ {0, \(-h\), 0}\  + 
            \((\ {{Cos[\(-phi\)], 0, Sin[\(-phi\)]}, {0, 1, 0}, {
                    \(-Sin[\(-phi\)]\), 0, Cos[\(-phi\)]}} . r)\); 
		

        		 (*\ Rotation\ operator\ Cn\ *) 
		Clear[Cn]; 
		
        Cn[r_] := 
          \ {{Cos[2  Pi/NC], 0, Sin[2  Pi/NC]}, {0, 1, 0}, {
                \(-Sin[2  Pi/NC]\), 0, Cos[2  Pi/NC]}} . r; 
		

        		 (*\ basis\ 
            \(atoms : \ Note\ the\ tube\ axis\ is\ the\ y - axis\)\ *) 
	

        		aatom = {Radius, 0, 0}; \ 
		dbasis = \((E1 + E2)\)/3\ ; 
		
        hbasis = 
          Sqrt[\((Cross[dbasis, {1, 0, 0}] . Cross[dbasis, {1, 0, 0}])\)]; 

        		phibasis = 2  Pi*\((dbasis . B)\)/\((B . B)\); 
		
        batom = 	\ {0, hbasis, 0}\  + 
					
            \((\ {{Cos[phibasis], 0, Sin[phibasis]}, {0, 1, 0}, {
                    \(-Sin[phibasis]\), 0, Cos[phibasis]}} . {Radius, 0, 0})
              \); \ 
	Print["phibasis = ", phibasis]; 
	

        		 (*\ nanotube\ y\ unit - cell\ *) 
		
        dR = GCD[2*\((n1)\) + n2, 2*\((n2)\) + n1]; 
		
        T = Sqrt[3]*2  Pi*Radius/dR; 
		
        Print["Translational symmetry along y: T = ", T]; 	
		
        NT = T/h; 
			
        Print["Corresponding to ", NT, \ " chiral operations"]; 	

        	 (*\ Generate\ bands\ from\ coordinates\ using\ symmetry\ operations
            \ *) 
		

        		 (*\ Note\ we\ use\ that\ we\ have\ chiral\ "inversion"\ 
            symmetry\ in\ the\ sense\ *) 

        		 (*\ that\ \((1/S)\)\ generates\ an\ equivalent\ atom
            \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) \ 
		
        
		 (*\ 
          We\ must\ find\ matrix\ elements\ which\ give\ non - 
            zero\ value\ on\ the\ B\ basis\ \(state : \)\ *) 

        		 (*\ Lets\ search\ among\ NS\ A - 
            equivalent\ coordinates\ \((excl . \ A\ it\ self)\)\ *) 
		
        NS = Sqrt[n1^2\  + \ n2^2] + 5; 

		
        abhamlist = 
          Flatten[

            				{
              Table[SKhamHOP12[aatom, Nest[Cn, Nest[invS, batom, i], j]], {j, 
                  0, NC - 1}, {i, 1, NS}], 
\ \ \ \ \ \ \ \ \ \ 
              Table[\ SKhamHOP12[aatom, 
                  Nest[Cn, Nest[\ \ \ \ \ \ \ S, batom, i], j]], {j, 0, 
                  NC - 1}, {i, 0, NS}]}, 2]; 


        			 (*\ Bands\ using\ chiral\ kappa\ vector\ *) 	

        		 (*\ table\ of\ plots\ for\ different\ rotation\ symmetry\ *) 

        					
        bandplot = 
          Table[Print[
              " ======= Bands for rotational quantum number M = ", M, 
              " ======= "]; \ 
		
            bands = Table[
			
                phaselist = 
                  Flatten[

                    							{	
                      Table[Exp[\(+Pi\)*I*kappa*i\ ]*
                          Exp[\(-2\) Pi*I*M*j/NC]\ , {j, 0, NC - 1}, {i, 1, 
                          NS}], 
								\ \ \ 
                      Table[Exp[\(-Pi\)*I*kappa*i\ ]*
                          Exp[\(-2\) Pi*I*M*j/NC]\ , {j, 0, NC - 1}, {i, 0, 
                          NS}]}, 2]; 
						
						
                dottmp = phaselist . abhamlist; \ 
					

                										
                HAM = \ dottmp\  + \ Conjugate[Transpose[dottmp]]\  + \ 
                    SKhamOnsite\ ; 		
		
			
							
                If[M == 0, \ ky = Mod[Pi*\((kappa)\)*T/h, 2  Pi], 

                  							\ 
                  ky = Mod[
                      Pi*\((kappa)\)*T/h\ \  - \ 
                        \((2  Pi/NC)\)*M*\((\((phi*T/h)\)/\((2  Pi/NC)\))\), 
                      2  Pi]]; 	
						
						
                eigtmp = 	Eigenvalues[HAM]; 
		
							
                Transpose[{Pi*kappa*{1, 1, 1, 1, 1, 1, 1, 1}, Chop[eigtmp]}], 
                
						{kappa, \(-1\), 1, 1/nkpts}]; \ 


            						 (*
              		{\(Transpose[{\(-ky\)*{1, 1, 1, 1, 1, 1, 1, 1}, eigtmp}], 
                \)\ *) 

            						 (*
              	\ 	Transpose[{\ ky*{1, 1, 1, 1, 1, 1, 1, 1}, eigtmp}]}\ *) 

            			 (*	
              ListPlot[Flatten[bands, 2], 
                PlotStyle -> {AbsolutePointSize[3]}, \ *) 

            				 (*		PlotLabel -> M, Frame -> True, 
              FrameLabel -> {\*""\!\(k\_y\)"", "ε"}, 
              \ *) 
					 (*\ PlotRange -> {{\(-Pi\), Pi}, Automatic}, 
              \ *) 
				 (*
              	FrameTicks -> {{\(-Pi\), \(-Pi\)/2, 0, Pi/2, Pi}, 
                    Automatic, {}, {}}]	\ *) 
				
						
            ListPlot[Flatten[bands, 1], PlotStyle -> {AbsolutePointSize[3]}, 
              \ 
						PlotLabel -> M, Frame -> True, 
              FrameLabel -> {"κ", "ε"}, 
              PlotRange -> {{\(-Pi\), Pi}, Automatic}, \ 
					
              FrameTicks -> {{\(-Pi\), \(-Pi\)/2, 0, Pi/2, Pi}, 
                  Automatic, {}, {}}]	
							, \ \ {M, 0, NC - 1}]; 
        

	Return[bandplot]; 
]
				
	\)\)

Example: (5,5) Tube

SKnanotubebands[5, 5, 30]

Example: (9,2) Tube

SKnanotubebands[9, 2, 200]