-- 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()