Examples of an exact diagonalization study of a two-orbital model
Part of package SNEG,http://nrgljubljana.ijs.si/sneg
Rok Zitko,rok.zitko@ijs.si, Feb 2007, Mar 2007

<< sneg`sneg`;

<< Graphics`Graphics`

<< diagonalize.m

Off[General::spell1];
snegfermionoperators[c, d]
snegrealconstants[ϵs, δd, J, U]

H = ϵs (number[c[1]] - 
                1) + δd (number[d[1]] - 1) + 
          U/2 pow[number[d[1]] - 1, 2] + J spinspin[c[1], d[1]] // 
        FullSimplify

basisops = {c[1], d[1]};
diagSetupBasisQSZ[];

diagSetupHam[];

(* 
      Sample calculation *) 
diagonalize[{ϵs -> 1, δd -> 0, J -> 1/5, 
        U -> 1}]

deg[{q_, sz_}] := 
        1;

\(\(\( (*\ Sum\ of\ degeneracies\ *) \)\(
\)\(sum1[eig, 
      1 &]

 (*\ 
      Number\ of\ states\ *) 

    sum1[eig, Length@First[#] &]

 (*\ 
      Number\ of\ states\ *) 

    Z = sum1expr[eig2, 1]

    9

    16

    16\)\)\)

\(\(\( (*\ Partition\ function\ *) \)\(
\)\(Z = 
      tdtr1expr[eig2, 1]

    tdtr1expr[eig2, pow[number[d[1]], 2]]/Z

    % /. beta -> 10. \)\)\)

sum2[
      eig2, 1 &]

  

    sum2[eig2, 
      Length @ First[#1]  Length @ 
            First[#2] &]

  

    sum2expr[eig2, 1]

Z = 
      sum2expr[eig2, 1]

    sum2expr[eig2, number[d[1]]]/Z

    sum2expr[eig2, pow[number[d[1]], 2]]/Z

Z = 
      tdtr2expr[eig2, 1]

    tdtr2expr[eig2, pow[number[d[1]], 2]]/Z

    % /. beta -> 0.

tdtr2exprALT[eig2, pow[number[d[1]], 2]]/Z

specd = spectral2expr[eig2, d[CR, 1, UP]]/Z
Chop @ Simplify[specd /. beta -> 100. ]

specc = spectral2expr[eig2, c[CR, 1, UP]]/Z
Chop @ Simplify[specc /. beta -> 100. ]

evl = {number[d[1]], pow[number[d[1]], 2], number[c[1]], 
          pow[number[c[1]], 2], 
          spinspin[c[1], d[1]]};
spl = {d[CR, 1, UP], c[CR, 1, UP]};
calculate[evl, spl]

Simplify[expvs /. beta -> 100. ] // 
      Chop
Simplify[specs /. beta -> 100. ] // Chop

ClearAll[expvresults, specresults];
doit[value_, CALCCODE_:  calculate] := 
        Module[{}, 
diagonalize[{ϵs -> 
                value, δd -> 0, J -> 1/5, 
              U -> 1}]; 
expvs = specs = 
              Infinity; 
CALCCODE[evl, 
            spl]; 
  
expvs = 
            Simplify[expvs /. beta -> 100. ] // 
              Chop; 
specs = 
            Simplify[specs /. beta -> 100. ] // 
              Chop; 
AppendTo[values, 
            value]; 
expvresults[value] = 
            expvs; 
specresults[value] = 
            specs;
];

values = {};
Scan[(Print[#]; doit[#]) &, 
      Table[x, {x, -1, 1, 0.2}]]

MatrixForm @ Map[{#, expvresults[#]} &, values]

rule[x_] := 
        weight_ Delta[omega + freq_. ] :> 
          weight  Point[{x, freq}];
cfn[x_] := RGBColor[x, 0, 1 - x];
pointswithsize[lst_] := 
        Flatten @ ( lst /. Plus -> List) /. 
          weight_ p_Point :> {cfn[weight], 
              PointSize[scale weight], p};

scale = 0.1;
pts = 
        Map[specresults[#] /. rule[#] &, 
          values];
pts1 = 
        pointswithsize @  pts[[All, 1]];
pts2 = pointswithsize @  pts[[All, 2]];

Show [ 
      GraphicsArray[ {Graphics [pts1, Frame -> True], 
          Graphics[pts2, Frame -> True]}], 
      ImageSize -> 12 72]

Faster approach

Zmatrix = 
      makematricesbzvc[1, bvc]

    x1 = makematricesbzvc[number[d[1]], bvc]

tdtr1FAST[eig, Zmatrix]
tdtr1expr[eig2, 1]

tdtr1FAST[eig, x1]
tdtr1expr[eig2, number[d[1]]]

makeallmatricesbzvc[1, bvc]
% /. {{qn1_?VectorQ, qn2_?VectorQ}, mat_?ArrayQ} :> {qn1, 
          mat}
makematricesbzvc[1, bvc]

cmat = makeallmatricesbzvc[c[CR, 1, UP], bvc]

spectral2FAST[eig, cmat]
Expand @ spectral2expr[eig2, c[CR, 1, UP]]

calculateSetup[evl, spl]

values = {};
Scan[(doit[#, calculateFAST]) &, 
      Table[x, {x, -0.4, 0.4, 0.02}]]

MatrixForm @ Map[{#, expvresults[#]} &, values]

{ϵs -> value, δd -> 0, J -> 1/5, 
      U -> 1}

scale = 0.03;
pts = 
        Map[specresults[#] /. rule[#] &, 
          values];
pts1 = pointswithsize @  pts[[All, 1]];
pts2 = 
        pointswithsize @  pts[[All, 2]];
Show [ 
      GraphicsArray[ {Graphics [pts1, Frame -> True], 
          Graphics[pts2, Frame -> True]}], 
      ImageSize -> 12 72]

values = {};
Scan[(doit[#, calculateFAST]) &, 
      Table[x, {x, 0.1, 0.2, 0.005}]]

scale = 0.03;
pts = 
        Map[specresults[#] /. rule[#] &, 
          values];
pts1 = pointswithsize @  pts[[All, 1]];
pts2 = 
        pointswithsize @  pts[[All, 2]];
Show [ 
      GraphicsArray[ {Graphics [pts1, Frame -> True], 
          Graphics[pts2, Frame -> True]}], 
      ImageSize -> 12 72]