0
$\begingroup$

I am trying to implement the continuum model Hamiltonian for transition metal dichalcogenide systems as described in the reference: PRL 122, 086402 (2019)

   (* Define the constants*)
   (* twist angle*)
   \[Theta] = 1.2 Degree;
   (* For MoTe2 V, psi, w*)
   V = 8 (* meV*);
   psi = -89.6 ;
   w = -8.5 (*meV*);
   Vz = -2 (* displacement field meV*);
   meff = 0.62 (* valence band effective mass*);
   a0 = 3.472 (*angstrom*);
   NN = 6(*truncate range*);

   \[Theta] = 
   1.2 Degree; (*\[Theta]\[Pi]/180conversion from degress to radian*)

    V = V/1000;
    w = w/1000;
    Vz = Vz/1000;
    psi = psi Degree;
    hbar = 4.135667*10^-15/2/\[Pi];
    me = 0.51099895*10^6;
    meff = meff*me;
    kin = -hbar^2/2/meff*9*10^16;
    a0 = 3.472 ;
    a0 = a0*10^-10;
    a1 = a0 {1, 0};
    a2 = a0 {Cos[\[Pi]/3], Sin[\[Pi]/3]};
    (* Reciprocal lattice vectors*)
    G1 = (4 \[Pi])/((Sqrt[3]) a0) {0, 1};
    G2 = (4 \[Pi])/((Sqrt[3]) a0) {-Sin[\[Pi]/3], Cos[\[Pi]/3]};
    (*G2=(4\[Pi])/((Sqrt[3])a0){Cos[\[Pi]/3],Sin[\[Pi]/3]};*)

    g1 = \[Theta] {G1[[2]], -G1[[1]]};
    g2 = \[Theta] {G2[[2]], -G2[[1]]};
    kD = g1[[1]]/Sqrt[3];

    Kb = {-kD*Sqrt[3]/2, -kD/2}; (* K botom layer *)
    Kt = {-kD*Sqrt[3]/2, 
     kD/2} ;   (* K top layer *)

     (* Define the lattice*)
      n = 6;
      nlayer = (2 n + 1)^2;
      L = Table[0, {i, 1, 2 (2 n + 1)^2}, {j, 1, 2}];
     invL = Table[0, {i, 1, 2 n + 1}, {j, 1, 2 n + 1}];
     k = 1;
     count = 0;
     Do[Do[
    (*bottom layer*)
     L[[k, 1]] = i; L[[k, 2]] = j;
     invL[[i + n + 1, j + n + 1]] = count;
     count = count + 1;

     (* top layer*)
     L[[k + nlayer, 1]] = i; L[[k + nlayer, 2]] = j;
     k = k + 1;
     , {j, -n, n}], {i, -n, n}]

      (* check the form of inVL*)
       MatrixForm[invL];

       (* Hamiltonian implementation *)
        Clear[kx, ky]

         HHHH[kx_, ky_] :=
          Module[{ h, NN},

      NN = 6;
      siteN = (2 NN + 1)^2;
      h = ConstantArray[0, {2 siteN, 2 siteN}];
     Do[ 
      ix = L[[i, 1]];
      iy = L[[i, 2]];

      (* diagnal elements in the Hamiltonian*)
  
       ax = kx - Kb[[1]] + ix g1[[1]] + iy g2[[1]];
         ay = ky - Kb[[2]] + ix g1[[2]] + iy g2[[2]];
         h[[i, i]] += kin (ax^2 + ay^2) + Vz/2;

       ax = kx - Kt[[1]] + ix g1[[1]] + iy g2[[1]];
         ay = ky - Kt[[2]] + ix g1[[2]] + iy g2[[2]];
         h[[i + siteN, i + siteN]] += kin (ax^2 + ay^2) - Vz/2;

        h[[i + siteN, i]] += w;
         h[[i, i + siteN]] += w;

        If[ix != NN,
         j = invL[[ix + 1 + NN + 1, iy + NN + 1]];
          h[[j + 1, i]] += V Exp[I*psi];
          h[[j + siteN + 1, i + siteN]] += V Exp[-I*psi];
         ];

          If [iy != -NN,
          j = invL[[ix + NN + 1, iy - 1 + NN + 1]];
                 h[[j + 1, i]] += V Exp[I psi];
                 h[[j + siteN + 1, i + siteN]] += V*Exp[-I psi];
                 h[[j + 1, i + siteN]] += w;
           ];

             If [ix != -NN && iy != NN,

               j = invL[[ix - 1 + NN + 1, iy + 1 + NN + 1]];
                 h[[j + 1, i]] += V*Exp[I*psi];
                 h[[j + siteN + 1, i + siteN]] += V*Exp[-I*psi];
                 h[[j + siteN + 1, i]] += w;
                 ];

            If [ix != -NN,
             j = invL[[ix - 1 + NN + 1, iy + NN + 1]];
                 h[[j + 1, i]] += V*Exp[-I*psi];
                h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
              ];
          If [iy != NN,
          j = invL[[ix + NN + 1, iy + 1 + NN + 1]];
                 h[[j + 1, i]] += V*Exp[-I*psi];
                h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
                h[[j + siteN + 1, i]] += w;
              ];

         If [ix != NN && iy != -NN,

                 j = invL[[ix + 1 + NN + 1, iy - 1 + NN + 1]];
                 h[[j + 1, i]] += V*Exp[-I*psi];
                 h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
                 h[[j + 1, i + siteN]] += w;];

               , {i, 1, siteN}];

             Return[h]

              ]

                Num = 50;
                KXX1 = N@Subdivide[g1[[1]]/2, 0, Num - 1];
                 KXX2 = N@Subdivide[0, -g1[[1]]/2, Num - 1];
                KXX3 = N@Subdivide[-g1[[1]]/2, -g1[[1]]/2, Num - 1];
                KXX4 = N@Subdivide[-g1[[1]]/2, g1[[1]]/2, (2 Num - 1)];
                KX = Join[KXX1, KXX2, KXX3, KXX4];

                KYY1 = N@Subdivide[-g1[[1]]/(Sqrt [3]/2), 0, Num - 1];
               KYY2 = N@Subdivide[0, g1[[1]]/(Sqrt[3]/2), Num - 1];
                KYY3 = N@
                Subdivide[g1[[1]]/(Sqrt[3]/2), -g1[[1]]/(Sqrt [3]/2), Num - 1];
              KYY4 = N@
              Subdivide[-g1[[1]]/(Sqrt [3]/2), -g1[[1]]/(Sqrt [3]/2), (2 Num - 
                    1)];
              KY = Join[KYY1, KYY2, KYY3, KYY4];

           Eigen = Table[0, {i1, Length[KX]}, {j1, 2 siteN}];

          Do[Eigen[[k]] = Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1, 
          Length[KX]}];

        out = ParallelTable[
         Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1, Length[KX]}];
           ListLinePlot[Transpose[out], PlotRange -> {{0, 250}, {0, 0.025}}]

The resulting band structure is look like this

enter image description here

I recently came across the github code which reproduce the band structure from the reference. (https://github.com/zihaophys/twisted_bilayer_graphene/blob/master/TMD/homobilayer/moireTMD.py)

I have studied and trying to reproduce in exact same way. Can someone please point out how to fix the band structure.

Thanks

New contributor
Shenaya is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
3
  • 1
    $\begingroup$ Welcome to the Mathematica Stack Exchange. The code you have uploaded does not evaluate or plot anything (yet). This stack site is about the technical software called Mathematica and the associated Wolfram Language. Stack sites are not project assistance sites. Please present a minimal example and ask an on-topic question that can be answered. Where exactly are you having difficulty? $\endgroup$
    – Syed
    Mar 17 at 4:49
  • 1
    $\begingroup$ There several posts about twisted bilayer graphene - see mathematica.stackexchange.com/search?q=graphene $\endgroup$ 2 days ago
  • $\begingroup$ Dear All, I was trying to get the continuum band structure for TMDs. So I did implement the Hamiltonian. My main problem is getting the band structure along the given high symmetry line. The band structure getting from this code is different than the one in the reference. I'm thinking may be the way I define the high symmetry line is not correct. Will look into you. Yes I did refer the examples about twisted bilayer system. Meantime I would really appreciate if someone can give me a comment who is familiar with continuum models of TMDs. Thanks $\endgroup$
    – Shenaya
    11 hours ago

0

Your Answer

Shenaya is a new contributor. Be nice, and check out our Code of Conduct.

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.