====== CalculateG ======

###
Warning! Page is under development.

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

===== Input =====

  * //HTB// : Tight-binding object, which can be created using the function //[[documentation:language_reference:functions:NewTightBinding|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 //[[documentation:language_reference:objects:responsefunction:functions:changetype|ResponseFunction.ChangeType()]]//). Default value is "ListOfPoles"

===== Output =====

  * G : Response Function in the Block List of Poles (If not chosen otherwise) 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 real numbers, 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//.
===== Example =====

###
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.
###

==== Input ====
<code Quanty 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()

</code>

==== Result ====
Diagonal Elements of the plots below show Partial DOS of Cu,O1,O2. 
{{ :documentation:language_reference:functions:pes_ij.png?1250 |}}

==== Used functions ====
<code Quanty function.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
</code>
===== Table of contents =====
{{indexmenu>.#1}}