From e547ed98c39b901ab20c22034d9cb20c5fb3d33a Mon Sep 17 00:00:00 2001 From: Anjishnu Bose <39477277+Anjishnubose@users.noreply.github.com> Date: Tue, 19 Mar 2024 14:11:15 -0400 Subject: [PATCH] fixed greens function at finite temperature?? --- src/BdGModel.jl | 53 ++++++++++++++++++++++++++++++------------------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/src/BdGModel.jl b/src/BdGModel.jl index cdb46bc..51fd033 100644 --- a/src/BdGModel.jl +++ b/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 \ No newline at end of file +# 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