CalculateG

Warning! Page is under development.

Function CalculateG(HTB, Options) calculates single-particle Green's function from a given Tight-Binding Object HTB.

• HTB : Tight-binding object, which can be created using the function NewTightBinding()
• Options : A table of options. Possible options are:

  • “Emin” : Real, minimum value of the energy. Default value is -10
  • “Emax” : Real, maximum value of the energy. Default value is 10
  • “NE” : Positive Integer defining number of grid points on the energy axis. Default value is 2000
  • “Nk” : Table of 3 integers {Nkx,Nky,Nkz}. Number of k-points along x,y,z directions. Default value {40,40,40}
  • “Type” : String, The type of output Green's function (see ResponseFunction.ChangeType()). Default value is “ListOfPoles”

• G : Response Function in the Block List of Poles representation $\{\{A_0, a_1, a_2,\dots,a_n\},\{ B_1,B_2, \dots, B_n \}, $ type=“ListOfPoles” $\dots \}$, which corresponds to

$$G(\omega) = A_0 + \sum_{k} \frac{B_k}{\omega-a_k+i\gamma/2}$$ where $a_1,\dots,a_n$ are reals, and $A_0,B_1,\dots, B_n$ are matrices of dimensions $N_O \times N_O $, where $N_O$ is the number of orbitals in HTB.Atoms.

This example creates Tight-Binding model for 2D layer of CuO2 also known as the Emery model. CalculateG() Is used to calculate and plot Green's functions. Code in the example relies on the user-written function ScalarResponseFunctionFromBlockListOfPolesResponseFunction() that is given at the end of the page.

Example.Quanty

-- define on-site energies and hopping parameters
ed = -1
ep1 = -3
ep2 = ep1 
t = 1 -- tdp
tdd = 0.1
tpp = 0.3
HTB = NewTightBinding()
HTB.Name = "Emery model"
 
HTB.Cell = {{1,0,0},{0,1,0},{0,0,1}}
 
HTB.Atoms = {{"Cu",{0.5,0.5,0},{{"d",{"x^2y^2"}}}},
             {"O1",{1,0.5,0},{{"p",{"x"}}}},
             {"O2",{0.5,1,0},{{"p",{"y"}}}}}
 
HTB.Hopping = { {"Cu.d","Cu.d",{0,0,0},{{ed}}},
                {"O1.p","O1.p",{0,0,0},{{ep1}}},
                {"O2.p","O2.p",{0,0,0},{{ep2}}},
                {"Cu.d","O1.p",{0.5,0,0},{{t}}},
                {"Cu.d","O1.p",{-0.5,0,0},{{t}}},
                {"O1.p","Cu.d",{0.5,0,0},{{t}}},
                {"O1.p","Cu.d",{-0.5,0,0},{{t}}},
                {"Cu.d","O2.p",{0,0.5,0},{{t}}},
                {"Cu.d","O2.p",{0,-0.5,0},{{t}}},
                {"O2.p","Cu.d",{0,0.5,0},{{t}}},
                {"O2.p","Cu.d",{0,-0.5,0},{{t}}},
                {"Cu.d","Cu.d",{1,0,0},{{tdd}}},
                {"Cu.d","Cu.d",{0,1,0},{{tdd}}},
                {"Cu.d","Cu.d",{-1,0,0},{{tdd}}},
                {"Cu.d","Cu.d",{0,-1,0},{{tdd}}},
                {"O1.p","O2.p",{0.5,0.5,0},{{tpp}}},
                {"O1.p","O2.p",{0.5,-0.5,0},{{tpp}}},
                {"O1.p","O2.p",{-0.5,0.5,0},{{tpp}}},
                {"O1.p","O2.p",{-0.5,-0.5,0},{{tpp}}},
                {"O2.p","O1.p",{0.5,0.5,0},{{tpp}}},
                {"O2.p","O1.p",{0.5,-0.5,0},{{tpp}}},
                {"O2.p","O1.p",{-0.5,0.5,0},{{tpp}}},
                {"O2.p","O1.p",{-0.5,-0.5,0},{{tpp}}}
            }
-- Calculate Block Greens Function
G0Block= CalculateG(HTB,{{"NE",1e4},{"Nk",{1000,1000,1}}})
-- Extract single orbital Green's functions for plotting
G0_Cu = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,1,1)
G0_O1 = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,2,2)
G0_O2 = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,3,3)
G0_Cu_O1 = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,1,2)
G0_Cu_O2 = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,1,3)
G0_O1_O2 = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,2,3)
G0_O2_O1 = ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,3,2)
G0_O1_Cu =  ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,2,1)
G0_O2_Cu =  ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0Block,3,1)
-- Plot Green's functions (Density of States)
Emin=-12
Emax=7
Ymin=-4
Ymax=3
 
Pl  = Graphics.Plot({G0_Cu,G0_Cu_O1,G0_Cu_O2,G0_O1_Cu,G0_O1,G0_O1_O2,G0_O2_Cu,G0_O2_O1,G0_O2,["GammaG"]=0.1}, {{"Nrow",3},{"NColumn",3},{"Frame",{{"Ymin",Ymin},{"Ymax",Ymax},{"Xmin",Emin},{"Xmax",Emax},{"dYTick",1},{"dXTick",2},{"FontSize",0.03},{"YFormat","%3.1f"},{"XLabel","E [t]"},{"YLabel","G"}}}})
PlSVG = Graphics.ToSVG(Pl)
file,err = io.open("DOS_ij.svg",'w')
file:write(PlSVG)
file:close()

Diagonal Elements of the plots below show Partial DOS of Cu,O1,O2.

Example.Quanty

-- this function extracts G_ij response function of Block List of Poles Response Function object
function ScalarResponseFunctionFromBlockListOfPolesResponseFunction(G0,i,j)
    local G0T = ResponseFunction.ToTable(G0)
    local k
    local A0_ij = G0T[1][1][i][j] 
    local ai = G0T[1] 
    table.remove(ai,1) --remove matrix A0
    table.insert(ai,1,A0_ij) -- add number A0_ij
    local bw_ij = {}
    for k=1,#G0T[2] do  -- from k=1 to NE
        bw_ij[#bw_ij+1] = G0T[2][k][i][j]
    end
    G_ij = ResponseFunction.New( {ai,bw_ij,mu=0,type="ListOfPoles", name="G0"} )
    return G_ij
end