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

$\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, 2023 at 4:49
  • 1
    $\begingroup$ There several posts about twisted bilayer graphene - see mathematica.stackexchange.com/search?q=graphene $\endgroup$ Mar 18, 2023 at 6:50
  • $\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
    Mar 20, 2023 at 0:32

1 Answer 1

0
$\begingroup$

The ky is wrong, it should be like KYY1 = N@Subdivide[-g1[[1]]/Sqrt [3]/2, 0, Num - 1] to create the correct one.

$\endgroup$
1
  • 3
    $\begingroup$ Welcome to Mathematica. I recommend that you show the answer that results from your correction. $\endgroup$
    – bbgodfrey
    Nov 25, 2023 at 17:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.