-- We define H0 as a 4 site tight binding model with periodic boundary conditions
--
--    ... (0) --- (3) --- (6) --- (9) --- (0) ...
-- 
H0 = NewOperator(12,0,{{  0, -3,1},{ 3, -0,1},
                       {  3, -6,1},{ 6, -3,1},
                       {  6, -9,1},{ 9, -6,1},
                       {  9, -0,1},{ 0, -9,1}})
print("The tight binding Hamiltonian for a 4 cite chain is given as")
print(H0)   
print("In matrix form")   
print(OperatorToMatrix(H0))        
          
-- We define H1 as coupling to the local sites
--     
--        (0)     (3)     (6)     (9)     
--         |       |       |       |       
--        (1)     (4)     (7)     (10)    
--         |       |       |       |       
--        (2)     (5)     (8)     (11)    
--
H1 = NewOperator(12,0,{{  0, -1,1},{ 1, -0,1},
                       {  1, -2,1},{ 2, -1,1},
                       
                       {  3, -4,1},{ 4, -3,1},
                       {  4, -5,1},{ 5, -4,1},
                       
                       {  6, -7,1},{ 7, -6,1},
                       {  7, -8,1},{ 8, -7,1},
                       
                       {  9,-10,1},{10, -9,1},
                       { 10,-11,1},{11,-10,1}})
print("The full Hamiltonian we take as the one dimensional chain with side chains")
print("The side chains are given by the Hamiltonian")
print(H1)
print("In matrix form")   
print(OperatorToMatrix(H1))  

-- We now can define the full Hamiltonian
H = H0 + H1
print("The full Hamiltonian is H0 + H1")
print(H)
print("In matrix form")   
print(OperatorToMatrix(H))  

print("We now calcualte the bare (G0) and full (G) vacuum Green's functions")
-- In order to define the Green's function G0 and G we need the vacuum state
psi0 = NewWavefunction(12,0,{{"000 000 000 000",1}})

-- And the transition operators
-- We want a^{dag}_k = sum_{x=1}^N \sqrt{1/N} e^{i k x} a^{dag}_x
adag={}
x = {0,1,2,3}
k = {0/4,1/4,2/4,3/4}
for ik=1,4 do
  adag[ik] = NewOperator(12,0,{{ 0, sqrt(1/4)*exp(I*x[1]*k[ik]*2*Pi)},
                               { 3, sqrt(1/4)*exp(I*x[2]*k[ik]*2*Pi)},
                               { 6, sqrt(1/4)*exp(I*x[3]*k[ik]*2*Pi)},
                               { 9, sqrt(1/4)*exp(I*x[4]*k[ik]*2*Pi)}})
end
-- Create Green's function
S0, G0 = CreateSpectra(H0, adag, psi0,{{"Tensor",true},{"DenseBorder",0}})
G0.Chop()
S,  G  = CreateSpectra(H , adag, psi0,{{"Tensor",true},{"DenseBorder",0}})
G.Chop()

print("G0 is")
print(G0)
print("in matrix form that becomes")
print(ResponseFunction.ToMatrix(G0))

print("G is")
print(G)
print("in matrix form that becomes")
print(ResponseFunction.ToMatrix(G))

print("The self energy can be calcualted as Sigma = G0^{-1} - G^{-1}")
Sigma = ResponseFunction.CalculateSelfEnergy(G0,G)
Sigma = ResponseFunction.ChangeType(Sigma,"Tri")
Sigma.Chop()

print(Sigma)

print("We now can compare G0(w), G(w) and G0(w-Sigma(w))")
omega = 0.1
Gamma = 0.3

print("G(omega)")
print(G(omega,Gamma/2))

print("G(omega-Sigma(omega))")
print(Matrix.Inverse( Matrix.Inverse(G0(omega,Gamma/2)) - Sigma(omega,Gamma/2) ))

print("G(omega)-Sigma(omega)")
print(G(omega,Gamma/2) - Matrix.Inverse( Matrix.Inverse(G0(omega,Gamma/2)) - Sigma(omega,Gamma/2) ))
print("The previous should have been zero but computer math is different. Have a look at")
print("0.1+0.2-0.3")
print(0.1+0.2-0.3)

print("If we \"Chop\" the previous result it becomes zero")
print(Chop( G(omega,Gamma/2) - Matrix.Inverse( Matrix.Inverse(G0(omega,Gamma/2)) - Sigma(omega,Gamma/2) )) )