fixed greens function at finite temperature??

Author: Anjishnubose

src/BdGModel.jl

@@ -6,7 +6,7 @@ module BdG
     using ..TightBindingToolkit.DesignUCell: ModifyIsotropicFields!
     using ..TightBindingToolkit.BZone:BZ, MomentumPhaseFFT
     using ..TightBindingToolkit.Hams:Hamiltonian, ModifyHamiltonianField!
-    
+
     using LinearAlgebra, Tullio, TensorCast, Logging, Statistics
 
     import ..TightBindingToolkit.TBModel:FindFilling, GetMu!, GetFilling!, GetGk!, GetGr!, SolveModel!, GetGap!, FreeEnergy
@@ -27,7 +27,7 @@ module BdG
     - `Gk      ::  Array{Matrix{ComplexF64}}` : An Array (corresponding to the grid of k-points in `BZ`) of Greens functions.
     - `Fk      ::  Array{Matrix{ComplexF64}}` : An Array (corresponding to the grid of k-points in `BZ`) of anomalous Greens functions.
 
-    Initialize this structure using 
+    Initialize this structure using
     ```julia
     BdGModel(uc_hop::UnitCell, uc_pair::UnitCell, bz::BZ, Ham::Hamiltonian ; T::Float64=1e-3, filling::Float64=-1.0, mu::Float64=0.0, stat::Int64=-1)
     ```
@@ -53,7 +53,7 @@ module BdG
         Fk      ::  Array{Matrix{ComplexF64}}
         Gr      ::  Array{Matrix{ComplexF64}}
         Fr      ::  Array{Matrix{ComplexF64}}
-        
+
         BdGModel(uc_hop::UnitCell, uc_pair::UnitCell, bz::BZ, Ham::Hamiltonian ; T::Float64=1e-3, filling::Float64=-1.0, mu::Float64=0.0, stat::Int64=-1) = new{}(uc_hop, uc_pair, bz, Ham, T, filling, mu, -999.0, stat ,Array{Matrix{ComplexF64}}(undef, zeros(Int64, length(uc_hop.primitives))...), Array{Matrix{ComplexF64}}(undef, zeros(Int64, length(uc_hop.primitives))...), Array{Matrix{ComplexF64}}(undef, zeros(Int64, length(uc_hop.primitives))...),Array{Matrix{ComplexF64}}(undef, zeros(Int64, length(uc_hop.primitives))...))
         ##### Chosen the default value of filling to be -1 which is unphysical so that the code knows when filling has not been provided and has to be calculated from mu instead!
     end
@@ -75,10 +75,12 @@ Because of this, if you want to calculate the filling at a different chemical po
         Eks     =   reshape(getindex.(M.Ham.bands, Ref(1:N÷2)), prod(M.bz.gridSize))
         U11s    =   reshape(getindex.(M.Ham.states, Ref(1:N÷2), Ref(1:N÷2)), prod(M.bz.gridSize))
         U21s    =   reshape(getindex.(M.Ham.states, Ref(N÷2 + 1: N), Ref(1:N÷2)), prod(M.bz.gridSize))
+        U12s    =   reshape(getindex.(M.Ham.states, Ref(1:N÷2), Ref(N÷2 + 1: N)) , prod(M.bz.gridSize))
+        U22s    =   reshape(getindex.(M.Ham.states, Ref(N÷2 + 1: N), Ref(N÷2 + 1: N)) , prod(M.bz.gridSize))
 
         nFs     =   DistFunction.(Eks; T=M.T, mu=0.0, stat=M.stat)
 
-        @tullio filling := ((conj(U11s[k1][i, j]) * U11s[k1][i, j] * nFs[k1][j]) 
+        @tullio filling := ((conj(U11s[k1][i, j]) * U11s[k1][i, j] * nFs[k1][j])
                           + (conj(U21s[k1][i, j]) * U21s[k1][i, j] * (1 - nFs[k1][j])))
 
         return real(filling) / (prod(M.bz.gridSize) * M.uc_hop.localDim * length(M.uc_hop.basis))
@@ -109,9 +111,9 @@ GetMu!(M::BdGModel ;  tol::Float64=0.001)
 Function to get the correct chemical potential given a filling.
 """
     function GetMu!(M::BdGModel ;  guess::Float64 = 0.0, mu_tol::Float64 = 0.001, filling_tol::Float64 = 1e-6)
-        M.mu    =   BinarySearch(M.filling, M.Ham.bandwidth, FindFilling, (M,) ; initial = guess, x_tol = mu_tol, target_tol = filling_tol)
+        M.mu    =   BinarySearch(M.filling, (2*M.Ham.bandwidth[1], 2*M.Ham.bandwidth[2]), FindFilling, (M,) ; initial = guess, x_tol = mu_tol, target_tol = filling_tol)
         @info "Found chemical potential μ = $(M.mu) for given filling = $(M.filling)."
-    
+
     end
 
 
@@ -127,16 +129,16 @@ Finding the Greens functions, and anomalous greens functions in momentum space a
         Eks     =   reshape(getindex.(M.Ham.bands, Ref(1:N÷2)) , prod(M.bz.gridSize))   ##### Only the negative energies from the bdG spectrum
 
         U11s    =   reshape(getindex.(M.Ham.states, Ref(1:N÷2), Ref(1:N÷2)) , prod(M.bz.gridSize))          ##### The 4 different quadrants in the Unitary relating the BdG quasiparticles and the nambu basis
-        U21s    =   reshape(getindex.(M.Ham.states, Ref(N÷2 + 1: N), Ref(1:N÷2)) , prod(M.bz.gridSize)) 
-        U12s    =   reshape(getindex.(M.Ham.states, Ref(1:N÷2), Ref(N÷2 + 1: N)) , prod(M.bz.gridSize)) 
+        U21s    =   reshape(getindex.(M.Ham.states, Ref(N÷2 + 1: N), Ref(1:N÷2)) , prod(M.bz.gridSize))
+        U12s    =   reshape(getindex.(M.Ham.states, Ref(1:N÷2), Ref(N÷2 + 1: N)) , prod(M.bz.gridSize))
         U22s    =   reshape(getindex.(M.Ham.states, Ref(N÷2 + 1: N), Ref(N÷2 + 1: N)) , prod(M.bz.gridSize))
 
         nFs     =   DistFunction.(Eks; T=M.T, mu=0.0, stat=M.stat)
 
-        @reduce Gk[k1][i, j] |= sum(l) ((conj(U11s[k1][i, l]) * U11s[k1][j, l] * nFs[k1][l]) 
-                                      + (conj(U12s[k1][i, l]) * U12s[k1][j, l] * (1 - nFs[k1][l])))
+        @reduce Gk[k1][i, j] |= sum(l) ((conj(U11s[k1][i, l]) * U11s[k1][j, l] * nFs[k1][l])
+                                      + (conj(U21s[k1][i, l]) * U21s[k1][j, l] * (1 - nFs[k1][l])))
 
-        @reduce Fk[k1][i, j] |= sum(l) ((conj(U11s[k1][i, l]) * U21s[k1][j, l] * nFs[k1][l]) 
+        @reduce Fk[k1][i, j] |= sum(l) ((conj(U11s[k1][i, l]) * U21s[k1][j, l] * nFs[k1][l])
                                       + (conj(U12s[k1][i, l]) * U22s[k1][j, l] * (1 - nFs[k1][l])))
         M.Gk    =   reshape(Gk, M.bz.gridSize...)
         M.Fk    =   reshape(Fk, M.bz.gridSize...)
@@ -152,14 +154,14 @@ Finding the equal-time Greens functions and anomalous Greens function in real sp
 
 """
     function GetGr!(M::BdGModel)
-        
+
         phaseShift   =  MomentumPhaseFFT(M.bz, M.uc_hop)
 
-        Gr           =  FFTArrayofMatrix(M.Gk)  
-        M.Gr         =  Gr .* phaseShift   
+        Gr           =  FFTArrayofMatrix(M.Gk)
+        M.Gr         =  Gr .* phaseShift
 
-        Fr           =  FFTArrayofMatrix(M.Fk)  
-        M.Fr         =  Fr .* phaseShift   
+        Fr           =  FFTArrayofMatrix(M.Fk)
+        M.Fr         =  Fr .* phaseShift
     end
 
 
@@ -169,9 +171,10 @@ SolveModel!(M::BdGModel)
 ```
 one-step function to find all the attributes in  BdGModel after it has been initialized.
 """
-    function SolveModel!(M::BdGModel ; mu_guess::Float64 = M.Ham.bandwidth[1] + M.filling * (M.Ham.bandwidth[2] - M.Ham.bandwidth[1]), get_correlations::Bool = true, get_gap::Bool = false, verbose::Bool = true, mu_tol::Float64 = 1e-3, filling_tol::Float64 = 1e-6)
+    function SolveModel!(M::BdGModel ; mu_guess::Float64 = 2*M.Ham.bandwidth[1] + 2*M.filling * (M.Ham.bandwidth[2] - M.Ham.bandwidth[1]), get_correlations::Bool = true, get_gap::Bool = false, verbose::Bool = true, mu_tol::Float64 = 1e-3, filling_tol::Float64 = 1e-6)
         @assert M.Ham.is_BdG==true "Use other format for pure hopping Hamiltonian"
-
+        # println(mu_guess)
+        # println(M.Ham.bandwidth)
         if M.filling<0    ##### Must imply that filling was not provided by user and hence needs to be calculated from given mu
             GetFilling!(M)
         else
@@ -188,7 +191,7 @@ one-step function to find all the attributes in  BdGModel after it has been init
             GetGk!(M)
             GetGr!(M)
         end
-        
+
         if verbose
             @info "System Filled!"
         end
@@ -225,4 +228,14 @@ Calculate the free energy of the given `BdGModel`.
     end
 
 
-end
+# function entropy(M::BdGModel)
+#     Es      =   2*reduce(vcat, M.Ham.bands)
+#     ps = DistFunction(Es; T=M.T, mu=M.mu, stat=M.stat)
+
+
+#     return S
+
+
+end
+
+# end