From 4bd57b6849ed74f888a9eb76ba7b8eaafccbe6dd Mon Sep 17 00:00:00 2001 From: Kun Chen Date: Thu, 8 Jun 2023 13:10:42 -0400 Subject: [PATCH 1/7] T matrix still has a problem --- src/interaction.jl | 63 ++++++++++++++++++++++++++++++++++++++ src/legendreinteraction.jl | 55 ++++++++++++++++++++++++++++++--- src/parameter.jl | 5 ++- src/polarization.jl | 40 ++++++++++++++++++++---- src/selfenergy.jl | 36 ++++++++++++++++++++++ test/ladder.jl | 4 +-- 6 files changed, 190 insertions(+), 13 deletions(-) diff --git a/src/interaction.jl b/src/interaction.jl index a45d56c..d9ae3a1 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -525,4 +525,67 @@ function RPA_total(q, n, param; pifunc=Polarization0_ZeroTemp, Vinv_Bare=coulomb return Ws, Wa end +""" + function Tmatrix(q, n, param; pifunc = Polarization0_ZeroTemp, landaufunc = landauParameterTakada, Vinv_Bare = coulombinv, regular = false, kwargs...) + +Dynamic part of the T-matrix. Returns the spin symmetric part and asymmetric part separately. + +#Arguments: + - q: momentum + - n: matsubara frequency given in integer s.t. ωn=2πTn (either one frequency or array of frequencies) + - param: other system parameters + - pifunc: caller to the polarization function + - landaufunc: caller to the Landau parameter (exchange-correlation kernel) + - Vinv_Bare: caller to the bare Coulomb interaction + - regular: regularized RPA or not + +# Return: +```math + \\frac{1} {\\frac{m}{4\\pi a_s} +\\Gamma(q, n)} - \\frac{4\\pi a_s}{m}. +``` +""" +function Tmatrix(q, n, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) + as = 4π * param.as / param.me + if param.as < 1e-6 + return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) + else + return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...)) + end +end + +""" + function Tmatrix_wrapped(Euv, rtol, sgrid::SGT, param; int_type=:ko, + pifunc=Polarization0_ZeroTemp, landaufunc=landauParameterTakada, Vinv_Bare=coulombinv, kwargs...) where {SGT} + +Return dynamic part and instant part of KO-interaction Green's function separately. Each part is a MeshArray with inner state +(1: spin symmetric part, 2: asymmetric part), and ImFreq and q-grid mesh. + +#Arguments: + - Euv: Euv of DLRGrid + - rtol: rtol of DLRGrid + - sgrid: momentum grid + - param: other system parameters + - pifunc: caller to the polarization function + - landaufunc: caller to the Landau parameter (exchange-correlation kernel) + - Vinv_Bare: caller to the bare Coulomb interaction +""" +function Tmatrix_wrapped(Euv, rtol, sgrid::SGT, param; + ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT} + # TODO: innerstate should be in the outermost layer of the loop. Hence, the functions such as KO and Vinv_Bare should be fixed with inner state as argument. + @unpack β = param + as = 4π * param.as / param.me + wn_mesh = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:none) + green_dyn = GreenFunc.MeshArray(wn_mesh, sgrid; dtype=ComplexF64) + green_tail = GreenFunc.MeshArray(wn_mesh, sgrid; dtype=ComplexF64) + + for (ki, k) in enumerate(sgrid) + for (ni, n) in enumerate(wn_mesh.grid) + green_tail[ni, ki] = Tmatrix(k, n, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...) + green_dyn[ni, ki] = Tmatrix(k, n, param; ladderfunc=ladderfunc, kwargs...) - green_tail[ni, ki] + end + end + + return green_dyn, green_tail +end + end diff --git a/src/legendreinteraction.jl b/src/legendreinteraction.jl index d06cc1a..13512be 100644 --- a/src/legendreinteraction.jl +++ b/src/legendreinteraction.jl @@ -169,7 +169,7 @@ function helper_function_grid(ygrid, intgrid, n::Int, W, param) end -struct DCKernel{C} +struct DCKernel{C,T} int_type::Symbol spin_state::Symbol channel::Int @@ -180,13 +180,16 @@ struct DCKernel{C} qgrids::Vector{CompositeGrid.Composite} dlrGrid::DLRGrid - kernel_bare::Array{Float64,2} - kernel::Array{Float64,3} + kernel_bare::Array{T,2} + kernel::Array{T,3} function DCKernel(int_type, spin_state, channel, param, kgrid::C, qgrids, dlrGrid, kernel_bare, kernel) where {C} - return new{C}(int_type, spin_state, channel, param, kgrid, qgrids, dlrGrid, kernel_bare, kernel) + return new{C,Float64}(int_type, spin_state, channel, param, kgrid, qgrids, dlrGrid, kernel_bare, kernel) end + function DCKernel{T}(int_type, spin_state, channel, param, kgrid::C, qgrids, dlrGrid, kernel_bare, kernel) where {C,T} + return new{C,T}(int_type, spin_state, channel, param, kgrid, qgrids, dlrGrid, kernel_bare, kernel) + end end function DCKernel_2d(param, Euv, rtol, Nk, maxK, minK, order, int_type, channel, spin_state=:auto; @@ -382,4 +385,48 @@ function DCKernel0(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state return DCKernel(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel) end + +function DCKernel_Ladder(param; Euv=param.EF * 100, rtol=1e-10, Nk=8, maxK=param.kF * 10, minK=param.kF * 1e-7, order=4, int_type=:rpa, spin_state=:auto, kwargs...) + return DCKernel_Ladder(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state; kwargs...) +end + +function DCKernel_Ladder(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state=:sigma; + kgrid=CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order), + kwargs...) + # use helper function + @unpack kF, β = param + channel = 0 + + if spin_state == :sigma + # for self-energy, always use ℓ=0 + channel = 0 + elseif spin_state == :auto + error("not implemented!") + end + + bdlr = DLRGrid(Euv, β, rtol, false, :ph) + qgrids = [CompositeGrid.LogDensedGrid(:gauss, [0.0, maxK], [k, kF], Nk, minK, order) for k in kgrid.grid] + qgridmax = maximum([qg.size for qg in qgrids]) + #println(qgridmax) + + kernel_bare = zeros(ComplexF64, (length(kgrid.grid), (qgridmax))) + kernel = zeros(ComplexF64, (length(kgrid.grid), (qgridmax), length(bdlr.n))) + + helper_grid = CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.1 * maxK], [0.0, kF], 2Nk, 0.01minK, 2order) + intgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, helper_grid[end]], [0.0, kF], 2Nk, 0.01minK, 2order) + + # dynamic + for (ni, n) in enumerate(bdlr.n) + helper = helper_function_grid(helper_grid, intgrid, 1, u -> Interaction.Tmatrix(u, n, param; kwargs...), param) + for (ki, k) in enumerate(kgrid.grid) + for (pi, p) in enumerate(qgrids[ki].grid) + Hp, Hm = Interp.interp1D(helper, helper_grid, k + p), Interp.interp1D(helper, helper_grid, abs(k - p)) + kernel[ki, pi, ni] = (Hp - Hm) + end + end + end + + return DCKernel{ComplexF64}(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel) +end + end diff --git a/src/parameter.jl b/src/parameter.jl index b431d3f..580e6ce 100644 --- a/src/parameter.jl +++ b/src/parameter.jl @@ -24,6 +24,7 @@ using ..Parameters EF::Float64 = 1.0 #kF^2 / (2me) β::Float64 = 200 # bare inverse temperature μ::Float64 = 1.0 + as::Float64 = 0.0 # s-wave scattering length (local contact interaction) # artificial parameters Λs::Float64 = 0.0 # Yukawa-type spin-symmetric interaction ~1/(q^2+Λs) @@ -130,7 +131,7 @@ end # end """ - function fullUnit(ϵ0, e0, me, EF, β) + function fullUnit(ϵ0, e0, me, EF, β, dim=3, spin=2; kwargs...) generate Para with a complete set of parameters, no value presumed. @@ -140,6 +141,8 @@ generate Para with a complete set of parameters, no value presumed. - me: electron mass - EF: Fermi energy - β: inverse temperature +- dim: dimension +- spin: spin """ @inline function fullUnit(ϵ0, e0, me, EF, β, dim=3, spin=2; kwargs...) # μ = try diff --git a/src/polarization.jl b/src/polarization.jl index 81e6c87..5115356 100644 --- a/src/polarization.jl +++ b/src/polarization.jl @@ -5,7 +5,7 @@ using ..Parameter using ..Parameters using ..Convention -export Polarization0_ZeroTemp, Polarization0_FiniteTemp +export Polarization0_ZeroTemp, Polarization0_FiniteTemp, Ladder0_FiniteTemp # Analytical calculated integrand of Π0 in 2D. @inline function _ΠT2d_integrand(k, q, ω, param) @@ -50,14 +50,21 @@ end @unpack me, β, μ = param # ω only appears as ω^2 so no need to check sign of ω k = x / (1 - x) - nk = 1.0 / (exp(β * (k^2 / 2 / me - μ)) + 1) + ek = β * (k^2 / 2 / me - μ) + nk = ek < 0.0 ? 1.0 / (exp(ek) + 1) : exp(-ek) / (1 + exp(-ek)) if q > 1e-6 * param.kF a = ω * 1im - k^2 / me - q^2 / 2 / me + 2μ b = q * k / me return me / (2π)^2 * (k / q * (log(a + b) - log(a - b)) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian else - return me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + if ω != 0.0 || abs(ek) > 1e-4 + return me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + else # ω==0.0 && abs(ek)<1e-4 + #expansion of 1/ek*(2/(exp(ek)+1)-1) + f = -0.5 + ek^2 / 24 - ek^4 / 240 + 17 * ek^6 / 40320 + return me / (2π)^2 * (2k^2 / me * (f / 2 * β) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + end end # if q is too small, use safe form @@ -90,14 +97,15 @@ function finitetemp_kgrid_ladder(μ::Float64, m::Float64, q::Float64, kF::Float6 mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF) mixX = mink / (1 + mink) if q >= sqrt(8 * μ * m) - kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0], [0.0, kF], scaleN, mink, gaussN) + x = kF / (1 + kF) + kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0 - 1e-8], [0.0, x], scaleN, mink, gaussN) else k1 = abs(q - sqrt(8 * μ * m - q^2)) / 2 k2 = (q + sqrt(8 * μ * m - q^2)) / 2 x1 = k1 / (1 + k1) x2 = k2 / (1 + k2) x3 = kF / (1 + kF) - kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0], sort([x1, x2, x3]), scaleN, mink, gaussN) + kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, 1.0 - 1e-8], sort([x1, x2, x3]), scaleN, mink, gaussN) end return kgrid end @@ -142,13 +150,21 @@ function Ladder0_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6 if dim == 3 for (ki, k) in enumerate(kgrid.grid) integrand[ki] = _LadderT3d_integrand(k, q, 2π * n / β, param) - @assert !isnan(integrand[ki]) "nan at k=$k, q=$q" + @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n" end end return Interp.integrate1D(integrand, kgrid) end +function Ladder0_FreeElectron(q::Float64, n::Int, param) + @unpack dim, kF, β, me, μ = param + if dim != 3 + error("No support for finite-temperature polarization in $dim dimension!") + end + return -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * n / β) +end + """ function Polarization0_FiniteTemp(q::Float64, n::Int, param, maxk=20, scaleN=20, minterval=1e-6, gaussN=10) @@ -252,6 +268,18 @@ function Ladder0_FiniteTemp(q::Float64, n::AbstractVector, param; scaleN=20, min return Interp.integrate1D(integrand, kgrid; axis=1) end +function Ladder0_FreeElectron(q::Float64, n::AbstractVector, param) + @unpack dim, kF, β, me, μ = param + if dim != 3 + error("No support for finite-temperature polarization in $dim dimension!") + end + integrand = zeros(ComplexF64, length(n)) + for (mi, m) in enumerate(n) + integrand[mi] = Ladder0_FreeElectron(q, m, param) + end + return integrand +end + """ function Polarization0_FiniteTemp(q::Float64, n::AbstractVector, param, maxk=20, scaleN=20, minterval=1e-6, gaussN=10) diff --git a/src/selfenergy.jl b/src/selfenergy.jl index 64b7511..1790fa3 100644 --- a/src/selfenergy.jl +++ b/src/selfenergy.jl @@ -255,6 +255,42 @@ function G0W0(param, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing; return G0W0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid; kwargs...) end +function G0Γ0(param, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing; Euv=100 * param.EF, rtol=1e-14, Nk=12, maxK=6 * param.kF, minK=1e-8 * param.kF, order=8, int_type=:none, + kwargs...) + return G0Γ0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid; kwargs...) +end + +function G0Γ0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing; kwargs...) + @unpack dim = param + # kernel = SelfEnergy.LegendreInteraction.DCKernel_old(param; + # Euv = Euv, rtol = rtol, Nk = Nk, maxK = maxK, minK = minK, order = order, int_type = int_type, spin_state = :sigma) + # kernel = SelfEnergy.LegendreInteraction.DCKernel0(param; + # Euv = Euv, rtol = rtol, Nk = Nk, maxK = maxK, minK = minK, order = order, int_type = int_type, spin_state = :sigma) + # G0 = G0wrapped(Euv, rtol, kernel.kgrid, param) + + kGgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order) + if dim == 2 + error("No support for G0Γ0 in 2d dimension!") + elseif dim == 3 + if isnothing(kgrid) + kernel = SelfEnergy.LegendreInteraction.DCKernel_Ladder(param; + Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kwargs...) + else + if (kgrid isa AbstractVector) + kgrid = SimpleG.Arbitrary{eltype(kgrid)}(kgrid) + end + kernel = SelfEnergy.LegendreInteraction.DCKernel_Ladder(param; + Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kgrid=kgrid, kwargs...) + end + G0 = G0wrapped(Euv, rtol, kGgrid, param) + Σ, Σ_ins = calcΣ_3d(G0, kernel) + else + error("No support for G0W0 in $dim dimension!") + end + + return Σ +end + """ function zfactor(param, Σ::GreenFunc.MeshArray; kamp=param.kF, ngrid=[0, 1]) diff --git a/test/ladder.jl b/test/ladder.jl index cc985cc..70aa7f1 100644 --- a/test/ladder.jl +++ b/test/ladder.jl @@ -1,9 +1,9 @@ using MCIntegration using ElectronGas -const para = Parameter.rydbergUnit(1.0 / 10, 4.0, 3) +const para = Parameter.rydbergUnit(1.0 / 100, 1.0, 3) const n = 0 -const k = 0.1 * para.kF +const k = 0.0 * para.kF function integrand(vars, config) n, k, para = config.userdata From e5067740e9dcc9a29003ed1bb322457b345f232d Mon Sep 17 00:00:00 2001 From: Kun Chen Date: Thu, 8 Jun 2023 14:58:49 -0400 Subject: [PATCH 2/7] G0Gamma0 may still have a bug --- src/interaction.jl | 22 +++++++++++----- src/legendreinteraction.jl | 53 ++++++++++++++++++++++++++++++++++---- src/selfenergy.jl | 13 +++++++--- 3 files changed, 72 insertions(+), 16 deletions(-) diff --git a/src/interaction.jl b/src/interaction.jl index d9ae3a1..3db0c56 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -574,18 +574,26 @@ function Tmatrix_wrapped(Euv, rtol, sgrid::SGT, param; # TODO: innerstate should be in the outermost layer of the loop. Hence, the functions such as KO and Vinv_Bare should be fixed with inner state as argument. @unpack β = param as = 4π * param.as / param.me - wn_mesh = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:none) - green_dyn = GreenFunc.MeshArray(wn_mesh, sgrid; dtype=ComplexF64) - green_tail = GreenFunc.MeshArray(wn_mesh, sgrid; dtype=ComplexF64) + wn_mesh_ph = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:ph) + green_dyn_r = GreenFunc.MeshArray(wn_mesh_ph, sgrid; dtype=ComplexF64) + green_tail_r = GreenFunc.MeshArray(wn_mesh_ph, sgrid; dtype=ComplexF64) + + wn_mesh_pha = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:pha) + green_dyn_i = GreenFunc.MeshArray(wn_mesh_pha, sgrid; dtype=ComplexF64) + green_tail_i = GreenFunc.MeshArray(wn_mesh_pha, sgrid; dtype=ComplexF64) for (ki, k) in enumerate(sgrid) - for (ni, n) in enumerate(wn_mesh.grid) - green_tail[ni, ki] = Tmatrix(k, n, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...) - green_dyn[ni, ki] = Tmatrix(k, n, param; ladderfunc=ladderfunc, kwargs...) - green_tail[ni, ki] + for (ni, n) in enumerate(wn_mesh_ph.grid) + green_tail_r[ni, ki] = real.(Tmatrix(k, n, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...)) + green_dyn_r[ni, ki] = real.(Tmatrix(k, n, param; ladderfunc=ladderfunc, kwargs...)) - green_tail_r[ni, ki] + end + for (ni, n) in enumerate(wn_mesh_pha.grid) + green_tail_i[ni, ki] = 1im .* imag.(Tmatrix(k, n, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...)) + green_dyn_i[ni, ki] = 1im .* imag.(Tmatrix(k, n, param; ladderfunc=ladderfunc, kwargs...)) - green_tail_i[ni, ki] end end - return green_dyn, green_tail + return green_dyn_r, green_dyn_i, green_tail_r, green_tail_i end end diff --git a/src/legendreinteraction.jl b/src/legendreinteraction.jl index 13512be..e5c11cc 100644 --- a/src/legendreinteraction.jl +++ b/src/legendreinteraction.jl @@ -386,11 +386,11 @@ function DCKernel0(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state end -function DCKernel_Ladder(param; Euv=param.EF * 100, rtol=1e-10, Nk=8, maxK=param.kF * 10, minK=param.kF * 1e-7, order=4, int_type=:rpa, spin_state=:auto, kwargs...) - return DCKernel_Ladder(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state; kwargs...) +function DCKernel_Ladder_r(param; Euv=param.EF * 100, rtol=1e-10, Nk=8, maxK=param.kF * 10, minK=param.kF * 1e-7, order=4, int_type=:rpa, spin_state=:auto, kwargs...) + return DCKernel_Ladder_r(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state; kwargs...) end -function DCKernel_Ladder(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state=:sigma; +function DCKernel_Ladder_r(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state=:sigma; kgrid=CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order), kwargs...) # use helper function @@ -417,7 +417,7 @@ function DCKernel_Ladder(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin # dynamic for (ni, n) in enumerate(bdlr.n) - helper = helper_function_grid(helper_grid, intgrid, 1, u -> Interaction.Tmatrix(u, n, param; kwargs...), param) + helper = helper_function_grid(helper_grid, intgrid, 1, u -> real(Interaction.Tmatrix(u, n, param; kwargs...)), param) for (ki, k) in enumerate(kgrid.grid) for (pi, p) in enumerate(qgrids[ki].grid) Hp, Hm = Interp.interp1D(helper, helper_grid, k + p), Interp.interp1D(helper, helper_grid, abs(k - p)) @@ -426,7 +426,50 @@ function DCKernel_Ladder(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin end end - return DCKernel{ComplexF64}(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel) + return DCKernel(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel) +end + +function DCKernel_Ladder_i(param; Euv=param.EF * 100, rtol=1e-10, Nk=8, maxK=param.kF * 10, minK=param.kF * 1e-7, order=4, int_type=:rpa, spin_state=:auto, kwargs...) + return DCKernel_Ladder_i(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state; kwargs...) +end + +function DCKernel_Ladder_i(param, Euv, rtol, Nk, maxK, minK, order, int_type, spin_state=:sigma; + kgrid=CompositeGrid.LogDensedGrid(:cheb, [0.0, maxK], [0.0, param.kF], Nk, minK, order), + kwargs...) + # use helper function + @unpack kF, β = param + channel = 0 + + if spin_state == :sigma + # for self-energy, always use ℓ=0 + channel = 0 + elseif spin_state == :auto + error("not implemented!") + end + + bdlr = DLRGrid(Euv, β, rtol, false, :ph) + qgrids = [CompositeGrid.LogDensedGrid(:gauss, [0.0, maxK], [k, kF], Nk, minK, order) for k in kgrid.grid] + qgridmax = maximum([qg.size for qg in qgrids]) + #println(qgridmax) + + kernel_bare = zeros(ComplexF64, (length(kgrid.grid), (qgridmax))) + kernel = zeros(ComplexF64, (length(kgrid.grid), (qgridmax), length(bdlr.n))) + + helper_grid = CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.1 * maxK], [0.0, kF], 2Nk, 0.01minK, 2order) + intgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, helper_grid[end]], [0.0, kF], 2Nk, 0.01minK, 2order) + + # dynamic + for (ni, n) in enumerate(bdlr.n) + helper = helper_function_grid(helper_grid, intgrid, 1, u -> imag(Interaction.Tmatrix(u, n, param; kwargs...)), param) + for (ki, k) in enumerate(kgrid.grid) + for (pi, p) in enumerate(qgrids[ki].grid) + Hp, Hm = Interp.interp1D(helper, helper_grid, k + p), Interp.interp1D(helper, helper_grid, abs(k - p)) + kernel[ki, pi, ni] = (Hp - Hm) + end + end + end + + return DCKernel(int_type, spin_state, channel, param, kgrid, qgrids, bdlr, kernel_bare, kernel) end end diff --git a/src/selfenergy.jl b/src/selfenergy.jl index 1790fa3..01c1768 100644 --- a/src/selfenergy.jl +++ b/src/selfenergy.jl @@ -273,22 +273,27 @@ function G0Γ0(param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid::Union{A error("No support for G0Γ0 in 2d dimension!") elseif dim == 3 if isnothing(kgrid) - kernel = SelfEnergy.LegendreInteraction.DCKernel_Ladder(param; + kernel_r = SelfEnergy.LegendreInteraction.DCKernel_Ladder_r(param; + Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kwargs...) + kernel_i = SelfEnergy.LegendreInteraction.DCKernel_Ladder_i(param; Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kwargs...) else if (kgrid isa AbstractVector) kgrid = SimpleG.Arbitrary{eltype(kgrid)}(kgrid) end - kernel = SelfEnergy.LegendreInteraction.DCKernel_Ladder(param; + kernel_r = SelfEnergy.LegendreInteraction.DCKernel_Ladder_r(param; + Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kgrid=kgrid, kwargs...) + kernel_i = SelfEnergy.LegendreInteraction.DCKernel_Ladder_i(param; Euv=Euv, rtol=rtol, Nk=Nk, maxK=maxK, minK=minK, order=order, int_type=int_type, spin_state=:sigma, kgrid=kgrid, kwargs...) end G0 = G0wrapped(Euv, rtol, kGgrid, param) - Σ, Σ_ins = calcΣ_3d(G0, kernel) + Σr, Σ_ins = calcΣ_3d(G0, kernel_r) + Σi, Σ_ins = calcΣ_3d(G0, kernel_i) else error("No support for G0W0 in $dim dimension!") end - return Σ + return Σr, Σi end """ From db57ae05e7b2ad9d778e7d319f03b0ade3712535 Mon Sep 17 00:00:00 2001 From: njqmforum Date: Sat, 1 Jul 2023 23:35:47 -0400 Subject: [PATCH 3/7] Delta Tmatrix in imaginary time not smooth --- src/interaction.jl | 132 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 128 insertions(+), 4 deletions(-) diff --git a/src/interaction.jl b/src/interaction.jl index 3db0c56..eddb117 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -526,7 +526,7 @@ function RPA_total(q, n, param; pifunc=Polarization0_ZeroTemp, Vinv_Bare=coulomb end """ - function Tmatrix(q, n, param; pifunc = Polarization0_ZeroTemp, landaufunc = landauParameterTakada, Vinv_Bare = coulombinv, regular = false, kwargs...) + function oldTmatrix(q, n, param; pifunc = Polarization0_ZeroTemp, landaufunc = landauParameterTakada, Vinv_Bare = coulombinv, regular = false, kwargs...) Dynamic part of the T-matrix. Returns the spin symmetric part and asymmetric part separately. @@ -544,7 +544,7 @@ Dynamic part of the T-matrix. Returns the spin symmetric part and asymmetric par \\frac{1} {\\frac{m}{4\\pi a_s} +\\Gamma(q, n)} - \\frac{4\\pi a_s}{m}. ``` """ -function Tmatrix(q, n, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) +function oldTmatrix(q, n, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) as = 4π * param.as / param.me if param.as < 1e-6 return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) @@ -554,7 +554,7 @@ function Tmatrix(q, n, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs end """ - function Tmatrix_wrapped(Euv, rtol, sgrid::SGT, param; int_type=:ko, + function oldTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; int_type=:ko, pifunc=Polarization0_ZeroTemp, landaufunc=landauParameterTakada, Vinv_Bare=coulombinv, kwargs...) where {SGT} Return dynamic part and instant part of KO-interaction Green's function separately. Each part is a MeshArray with inner state @@ -569,7 +569,7 @@ Return dynamic part and instant part of KO-interaction Green's function separate - landaufunc: caller to the Landau parameter (exchange-correlation kernel) - Vinv_Bare: caller to the bare Coulomb interaction """ -function Tmatrix_wrapped(Euv, rtol, sgrid::SGT, param; +function oldTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT} # TODO: innerstate should be in the outermost layer of the loop. Hence, the functions such as KO and Vinv_Bare should be fixed with inner state as argument. @unpack β = param @@ -596,4 +596,128 @@ function Tmatrix_wrapped(Euv, rtol, sgrid::SGT, param; return green_dyn_r, green_dyn_i, green_tail_r, green_tail_i end + +""" + function T0(q::Float64, n::Int, param) + +Vacuum two-particle T-matrix in matsubara frequency. +#Arguments: + - q: total incoming momentum + - n: matsubara frequency given in integer s.t. Ωn=2πTn (either one frequency or array of frequencies) + - param: other system parameters +""" + +function T0matrix(q::Float64, n::Int, param) + β, me = param.β, param.me + as = 4π * param.as / me + B = me^(3/2)/(4π) + if as < 1e-6 + return as/(1. - as * B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β))) + else + return 1/((1/as) - B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β))) + end +end + +""" + function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) + +Many-body T-matrix in matsubara frequency. +#Arguments: + - q: total incoming momentum + - n: matsubara frequency given in integer s.t. Ωn=2πTn (either one frequency or array of frequencies) + - param: other system parameters +""" + +function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) + as = 4π * param.as / param.me + if param.as < 1e-6 + return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) + else + return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...)) + end +end + +""" + function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT} + +Difference between Many-body T-matrix and vacuum two-particle T-matrix in imaginary time. +#Arguments: +- Euv: Euv of DLRGrid +- rtol: rtol of DLRGrid +- sgrid: momentum grid +- param: other system parameters +- ladderfunc: Particle-Particle bubble (ladder-diagram of G0G0) +""" + +function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT} + β= param.β + Ωn_mesh_ph = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:ph) + δTmatrix_n = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64) + + for (ki, k) in enumerate(sgrid) + for (ni, n) in enumerate(Ωn_mesh_ph.grid) + δTmatrix_n[ni, ki] = (Tmatrix(k, n, param; ladderfunc = ladderfunc) - T0matrix(k, n, param)) + end + end + δTmatrix_tau = δTmatrix_n |> to_dlr |> to_imtime + return δTmatrix_tau +end + +""" + function T0maxtrix_wrapped(tau::Float64, q::Float64, param) + +Difference between Many-body T-matrix and vacuum two-particle T-matrix in imaginary time. +#Arguments: +- tau: imaginary time +- q: total incoming momentum +- param: other system parameters +""" + +function T0maxtrix_wrapped(tau::Float64, q::Float64, param) + """ Two particle scattering in a vacuum T-matrix in imaginary time """ + β, me, as = param.β, param.me, param.as + B = me^(3/2)/(4π) + m = me * as^2 + + # Pole Term contribution + poletermnumerator = exp(tau/ m) + poletermdenominator = 1- exp(β * (1. / (me*as^2) - q^2/ (4 * me))) + poleterm = - poletermnumerator/ poletermdenominator + # Integral Terms + + function integrand1(vars, q, tau, param) + me, as = param.me, param.as + m = me * as^2 + u = vars[1][1] + y = u/(1-u) + fraction = y^2/(y^2 + tau/ m) + return 2 * exp(-y^2) * fraction + end + + function integrand2(vars, q, tau, param) + β, me = param.β, param.me + u = vars[1][1] + y = u/(1-u) + factor = 1/(exp(β * (y^2/ tau + q^2/ (4 * me)))-1) + return factor * integrand1(vars, q, tau, param) + end + + tgrid = CompositeGrid.LogDensedGrid( + :gauss,# The top layer grid is :gauss, optimized for integration. For interpolation use :cheb + [0.0, 1.0],# The grid is defined on [0.0, 1.0] s.t. the integral in the original variable is from [0.0 , infinity] + [0.0],# and is densed at 0.0 and β, as given by 2nd and 3rd parameter. + 5,# N of log grid + 0.005, # minimum interval length of log grid + 5 # N of bottom layer + ) + data = [(integrand1(t, q, tau, param)+integrand2(t, q, tau, param)) for (ti, t) in enumerate(tgrid.grid)] + int_result = Interp.integrate1D(data, tgrid) + integralfactor = 1/(B * pi * sqrt(tau)) + + # Propagator factor + propfactor = exp(-q^2/ (4 * me) * tau) + + return propfactor*(poleterm + int_result * integralfactor) +end + end From 6c33dbb8583a485309fb156a5b97254e696021c7 Mon Sep 17 00:00:00 2001 From: njqmforum Date: Sun, 2 Jul 2023 10:34:13 -0400 Subject: [PATCH 4/7] Subtracted the vacuum ladder diagrams, had to change as in parameter to see nonzero plot --- src/interaction.jl | 22 ++++++------ src/parameter.jl | 2 +- src/polarization.jl | 86 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 99 insertions(+), 11 deletions(-) diff --git a/src/interaction.jl b/src/interaction.jl index eddb117..ffdce85 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -612,28 +612,30 @@ function T0matrix(q::Float64, n::Int, param) as = 4π * param.as / me B = me^(3/2)/(4π) if as < 1e-6 - return as/(1. - as * B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β))) + return as/(1. - as * B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β * im))) else - return 1/((1/as) - B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β))) + return 1/((1/as) - B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β * im))) end end """ - function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) + function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, vacuumladderfunc=Polarization.LadderVacuum0_FiniteTemp, kwargs...) Many-body T-matrix in matsubara frequency. #Arguments: - q: total incoming momentum - n: matsubara frequency given in integer s.t. Ωn=2πTn (either one frequency or array of frequencies) - param: other system parameters + - ladderfunc: Particle-particle bubble (ladder-diagram of G0G0) + - vacuumladderfunc: Vacuum particle-particle bubble """ -function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) +function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, vacuumladderfunc=Polarization.Ladder0Vacuum_FiniteTemp, kwargs...) as = 4π * param.as / param.me if param.as < 1e-6 - return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) + return as ./ (1 .+ as * (ladderfunc(q, n, param; kwargs...) - vacuumladderfunc(q, n, param; kwargs...))) else - return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...)) + return 1 ./ (1 / as .+ (ladderfunc(q, n, param; kwargs...) - vacuumladderfunc(q, n, param; kwargs...))) end end @@ -649,14 +651,14 @@ Difference between Many-body T-matrix and vacuum two-particle T-matrix in imagin - ladderfunc: Particle-Particle bubble (ladder-diagram of G0G0) """ -function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT} +function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, vacuumladderfunc=Polarization.Ladder0Vacuum_FiniteTemp, kwargs...) where {SGT} β= param.β Ωn_mesh_ph = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:ph) δTmatrix_n = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64) for (ki, k) in enumerate(sgrid) for (ni, n) in enumerate(Ωn_mesh_ph.grid) - δTmatrix_n[ni, ki] = (Tmatrix(k, n, param; ladderfunc = ladderfunc) - T0matrix(k, n, param)) + δTmatrix_n[ni, ki] = (Tmatrix(k, n, param; ladderfunc = ladderfunc, vacuumladderfunc = vacuumladderfunc) - T0matrix(k, n, param)) end end δTmatrix_tau = δTmatrix_n |> to_dlr |> to_imtime @@ -664,7 +666,7 @@ function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization end """ - function T0maxtrix_wrapped(tau::Float64, q::Float64, param) + function T0maxtrix_imtime(tau::Float64, q::Float64, param) Difference between Many-body T-matrix and vacuum two-particle T-matrix in imaginary time. #Arguments: @@ -673,7 +675,7 @@ Difference between Many-body T-matrix and vacuum two-particle T-matrix in imagin - param: other system parameters """ -function T0maxtrix_wrapped(tau::Float64, q::Float64, param) +function T0maxtrix_imtime(tau::Float64, q::Float64, param) """ Two particle scattering in a vacuum T-matrix in imaginary time """ β, me, as = param.β, param.me, param.as B = me^(3/2)/(4π) diff --git a/src/parameter.jl b/src/parameter.jl index 580e6ce..aa67222 100644 --- a/src/parameter.jl +++ b/src/parameter.jl @@ -24,7 +24,7 @@ using ..Parameters EF::Float64 = 1.0 #kF^2 / (2me) β::Float64 = 200 # bare inverse temperature μ::Float64 = 1.0 - as::Float64 = 0.0 # s-wave scattering length (local contact interaction) + as::Float64 = 1.0 # s-wave scattering length (local contact interaction) # artificial parameters Λs::Float64 = 0.0 # Yukawa-type spin-symmetric interaction ~1/(q^2+Λs) diff --git a/src/polarization.jl b/src/polarization.jl index 5115356..85ff666 100644 --- a/src/polarization.jl +++ b/src/polarization.jl @@ -81,6 +81,43 @@ end # end end +# Analytical calculated integrand of the vacuum ladder in 3D. +@inline function _LadderTVacuum3d_integrand(x, q, ω, param) + me, β= param.me, param.β + # μ = 0.0 for vacuum + μ = 0.0 + # ω only appears as ω^2 so no need to check sign of ω + k = x / (1 - x) + ek = β * (k^2 / 2 / me - μ) + nk = ek < 0.0 ? 1.0 / (exp(ek) + 1) : exp(-ek) / (1 + exp(-ek)) + + if q > 1e-6 * param.kF + a = ω * 1im - k^2 / me - q^2 / 2 / me + 2μ + b = q * k / me + return me / (2π)^2 * (k / q * (log(a + b) - log(a - b)) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + else + if ω != 0.0 || abs(ek) > 1e-4 + return me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + else # ω==0.0 && abs(ek)<1e-4 + #expansion of 1/ek*(2/(exp(ek)+1)-1) + f = -0.5 + ek^2 / 24 - ek^4 / 240 + 17 * ek^6 / 40320 + return me / (2π)^2 * (2k^2 / me * (f / 2 * β) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + end + end + + # if q is too small, use safe form + # if q < 1e-16 && ω ≈ 0 + # if abs(q - 2 * k)^2 < 1e-16 + # return 0.0 + # else + # return -k * me / (4 * π^2) * nk * ((8 * k) / ((q - 2 * k)^2)) + # end + # elseif q < 1e-16 && !(ω ≈ 0) + # return -k * me / (4 * π^2) * nk * ((8 * k * q^2) / (4 * me^2 * ω^2 + (q^2 - 2 * k * q)^2)) + # else + # return -k * me / (4 * π^2 * q) * nk * log1p((8 * k * q^3) / (4 * me^2 * ω^2 + (q^2 - 2 * k * q)^2)) + # end +end function finitetemp_kgrid(q::Float64, kF::Float64, maxk=20, scaleN=20, minterval=1e-6, gaussN=10) mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF) @@ -165,6 +202,55 @@ function Ladder0_FreeElectron(q::Float64, n::Int, param) return -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * n / β) end +""" + function Ladder0Vacuum_FiniteTemp(q::Float64, n::Int, param, scaleN=20, minterval=1e-6, gaussN=10) + +Finite temperature vacuum ladder function for matsubara frequency and momentum. Analytically sum over total incoming frequency and angular +dependence of momentum, and numerically calculate integration of magnitude of momentum. +Assume G_vacuum^{-1} = iω_n - (k^2/(2m)) + +The ladder function is defined as +```math +\\int \\frac{d^3 \\vec{p}}{\\left(2π^3\\right)} T \\sum_{i ω_n} \\frac{1}{iω_n+iΩ_n-\\frac{(\\vec{k}+\\vec{p})^2}{2 m}} \\frac{1}{-iω_n-\\frac{p^2}{2 m}} - \\frac{m}{2π^2}Λ +``` +where we subtract the UV divergent term that is a constant proportional to the UV cutoff Λ. + +#Arguments: + - q: momentum + - n: matsubara frequency given in integer s.t. Ωn=2πTn + - param: other system parameters + - maxk: optional, upper limit of integral -> maxk*kF + - scaleN: optional, N of Log grid in LogDensedGrid, check CompositeGrids for more detail + - minterval: optional, actual minterval of grid is this value times min(q,kF) + - gaussN: optional, N of GaussLegendre grid in LogDensedGrid. +""" +function Ladder0Vacuum_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6, gaussN=10) + dim, kF, β, me = param.dim, param.kF, param.β, param.me + # μ = 0.0 for vacuum + μ = 0.0 + if dim != 3 + error("No support for finite-temperature polarization in $dim dimension!") + end + + # check sign of q, use -q if negative + if q < 0 + q = -q + end + + # mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF) + # kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, maxk * kF], [0.5 * q, kF], scaleN, mink, gaussN) + kgrid = finitetemp_kgrid_ladder(μ, me, q, kF, scaleN, minterval, gaussN) + integrand = zeros(ComplexF64, kgrid.size) + if dim == 3 + for (ki, k) in enumerate(kgrid.grid) + integrand[ki] = _LadderTVacuum3d_integrand(k, q, 2π * n / β, param) + @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n" + end + end + + return Interp.integrate1D(integrand, kgrid) +end + """ function Polarization0_FiniteTemp(q::Float64, n::Int, param, maxk=20, scaleN=20, minterval=1e-6, gaussN=10) From 10135d33aac77052198bb59b9945149d19ad1c37 Mon Sep 17 00:00:00 2001 From: njqmforum Date: Tue, 11 Jul 2023 13:01:44 +0100 Subject: [PATCH 5/7] replaced high frequency tail of the finite temperature ladder with the two-particle ladder contribution --- src/interaction.jl | 12 +++++------- src/parameter.jl | 2 +- src/polarization.jl | 38 +++++++++++++++++++++++++++++--------- 3 files changed, 35 insertions(+), 17 deletions(-) diff --git a/src/interaction.jl b/src/interaction.jl index ffdce85..bf13156 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -627,15 +627,13 @@ Many-body T-matrix in matsubara frequency. - n: matsubara frequency given in integer s.t. Ωn=2πTn (either one frequency or array of frequencies) - param: other system parameters - ladderfunc: Particle-particle bubble (ladder-diagram of G0G0) - - vacuumladderfunc: Vacuum particle-particle bubble """ -function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, vacuumladderfunc=Polarization.Ladder0Vacuum_FiniteTemp, kwargs...) +function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) as = 4π * param.as / param.me if param.as < 1e-6 - return as ./ (1 .+ as * (ladderfunc(q, n, param; kwargs...) - vacuumladderfunc(q, n, param; kwargs...))) - else - return 1 ./ (1 / as .+ (ladderfunc(q, n, param; kwargs...) - vacuumladderfunc(q, n, param; kwargs...))) + return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) + return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...) ) end end @@ -651,14 +649,14 @@ Difference between Many-body T-matrix and vacuum two-particle T-matrix in imagin - ladderfunc: Particle-Particle bubble (ladder-diagram of G0G0) """ -function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, vacuumladderfunc=Polarization.Ladder0Vacuum_FiniteTemp, kwargs...) where {SGT} +function δTmatrix_wrapped(Euv, rtol, sgrid::SGT, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) where {SGT} β= param.β Ωn_mesh_ph = GreenFunc.ImFreq(β, BOSON; Euv=Euv, rtol=rtol, symmetry=:ph) δTmatrix_n = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64) for (ki, k) in enumerate(sgrid) for (ni, n) in enumerate(Ωn_mesh_ph.grid) - δTmatrix_n[ni, ki] = (Tmatrix(k, n, param; ladderfunc = ladderfunc, vacuumladderfunc = vacuumladderfunc) - T0matrix(k, n, param)) + δTmatrix_n[ni, ki] = (Tmatrix(k, n, param; ladderfunc = ladderfunc) - T0matrix(k, n, param)) end end δTmatrix_tau = δTmatrix_n |> to_dlr |> to_imtime diff --git a/src/parameter.jl b/src/parameter.jl index aa67222..c5634c3 100644 --- a/src/parameter.jl +++ b/src/parameter.jl @@ -24,7 +24,7 @@ using ..Parameters EF::Float64 = 1.0 #kF^2 / (2me) β::Float64 = 200 # bare inverse temperature μ::Float64 = 1.0 - as::Float64 = 1.0 # s-wave scattering length (local contact interaction) + as::Float64 = -2.5 # s-wave scattering length (local contact interaction) # artificial parameters Λs::Float64 = 0.0 # Yukawa-type spin-symmetric interaction ~1/(q^2+Λs) diff --git a/src/polarization.jl b/src/polarization.jl index 85ff666..f99b345 100644 --- a/src/polarization.jl +++ b/src/polarization.jl @@ -59,7 +59,7 @@ end return me / (2π)^2 * (k / q * (log(a + b) - log(a - b)) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian else if ω != 0.0 || abs(ek) > 1e-4 - return me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian + me / (2π)^2 * (2k^2 / me / (ω * 1im - k^2 / me + 2μ) * (2 * nk - 1) - 2) / (1 - x)^2 # additional factor 1/(1-x)^2 is from the Jacobian else # ω==0.0 && abs(ek)<1e-4 #expansion of 1/ek*(2/(exp(ek)+1)-1) f = -0.5 + ek^2 / 24 - ek^4 / 240 + 17 * ek^6 / 40320 @@ -179,24 +179,32 @@ function Ladder0_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6 if q < 0 q = -q end - + Ωn = 2π * n / β # mink = (q < 1e-16 / minterval) ? minterval * kF : minterval * min(q, kF) # kgrid = CompositeGrid.LogDensedGrid(:gauss, [0.0, maxk * kF], [0.5 * q, kF], scaleN, mink, gaussN) kgrid = finitetemp_kgrid_ladder(μ, me, q, kF, scaleN, minterval, gaussN) integrand = zeros(ComplexF64, kgrid.size) if dim == 3 for (ki, k) in enumerate(kgrid.grid) - integrand[ki] = _LadderT3d_integrand(k, q, 2π * n / β, param) - @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n" + if Ωn <1000 + integrand[ki] = _LadderT3d_integrand(k, q, Ωn, param) + @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n" + else + integrand[ki] = 0 # Do not calculate integrand for high frequency tail + end end end - - return Interp.integrate1D(integrand, kgrid) + result = Interp.integrate1D(integrand, kgrid) + if Ωn> 1000 + result = -me^(3 / 2) / (4π) * sqrt(Complex(q^2 / (4me) - 1im * Ωn - 2μ)) # Replace high frequency tail with two-particle particle-particle term + end + return result end function Ladder0_FreeElectron(q::Float64, n::Int, param) @unpack dim, kF, β, me, μ = param if dim != 3 + error("No support for finite-temperature polarization in $dim dimension!") end return -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * n / β) @@ -345,13 +353,25 @@ function Ladder0_FiniteTemp(q::Float64, n::AbstractVector, param; scaleN=20, min if dim == 3 for (ki, k) in enumerate(kgrid.grid) for (mi, m) in enumerate(n) - integrand[ki, mi] = _LadderT3d_integrand(k, q, 2π * m / β, param) - @assert !isnan(integrand[ki, mi]) "nan at k=$k, q=$q" + Ωm = 2π * m / β + if Ωm < 1000 + integrand[ki, mi] = _LadderT3d_integrand(k, q, Ωm, param) + @assert !isnan(integrand[ki, mi]) "nan at k=$k, q=$q" + else + integrand[ki, mi] = 0 # Don't integrate for high frequencies, we use abrupt two-particle expression first + end end end end + result = Interp.integrate1D(integrand, kgrid; axis=1) + for (mi, m) in enumerate(n) + Ωm = 2π * m / β + if Ωm > 1000 + result[mi] = -me^(3 / 2) / (4π) * sqrt(Complex(q^2 / (4me) - 1im * Ωm - 2μ)) # Replace with two-particle particle-particle high frequency tail + end + end - return Interp.integrate1D(integrand, kgrid; axis=1) + return result end function Ladder0_FreeElectron(q::Float64, n::AbstractVector, param) From d51b8c1ff62c41b4928e0733d1858e05484648a9 Mon Sep 17 00:00:00 2001 From: LiamLau1 Date: Tue, 11 Jul 2023 13:22:45 +0100 Subject: [PATCH 6/7] slight cleanup --- src/interaction.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/interaction.jl b/src/interaction.jl index bf13156..db6c47d 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -685,7 +685,7 @@ function T0maxtrix_imtime(tau::Float64, q::Float64, param) poleterm = - poletermnumerator/ poletermdenominator # Integral Terms - function integrand1(vars, q, tau, param) + function integrand1(vars, tau, param) me, as = param.me, param.as m = me * as^2 u = vars[1][1] From d188e2adc268b01ebc7a5be119cfb72cbe7d0f7a Mon Sep 17 00:00:00 2001 From: LiamLau1 Date: Wed, 4 Oct 2023 15:22:09 -0400 Subject: [PATCH 7/7] Added galilean invariance test file in test which calculates the self energy from T-matrix to first order --- src/interaction.jl | 30 +++++----- src/parameter.jl | 2 +- src/polarization.jl | 16 +++--- test/galileaninvariance.jl | 111 +++++++++++++++++++++++++++++++++++++ 4 files changed, 137 insertions(+), 22 deletions(-) create mode 100644 test/galileaninvariance.jl diff --git a/src/interaction.jl b/src/interaction.jl index db6c47d..fa8f406 100644 --- a/src/interaction.jl +++ b/src/interaction.jl @@ -598,7 +598,7 @@ end """ - function T0(q::Float64, n::Int, param) + function T0(q::Float64, n::Int, param; ; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...) Vacuum two-particle T-matrix in matsubara frequency. #Arguments: @@ -607,14 +607,13 @@ Vacuum two-particle T-matrix in matsubara frequency. - param: other system parameters """ -function T0matrix(q::Float64, n::Int, param) - β, me = param.β, param.me - as = 4π * param.as / me - B = me^(3/2)/(4π) - if as < 1e-6 - return as/(1. - as * B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β * im))) + +function T0matrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...) + as = (4π * param.as) / param.me + if param.as < 1e-6 + return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) else - return 1/((1/as) - B * sqrt(Complex(q^2/ (4 * me) - 2π * n/β * im))) + return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...) ) end end @@ -630,9 +629,10 @@ Many-body T-matrix in matsubara frequency. """ function Tmatrix(q::Float64, n::Int, param; ladderfunc=Polarization.Ladder0_FiniteTemp, kwargs...) - as = 4π * param.as / param.me + as = (4π * param.as) / param.me if param.as < 1e-6 return as ./ (1 .+ as * ladderfunc(q, n, param; kwargs...)) + else return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...) ) end end @@ -673,7 +673,7 @@ Difference between Many-body T-matrix and vacuum two-particle T-matrix in imagin - param: other system parameters """ -function T0maxtrix_imtime(tau::Float64, q::Float64, param) +function T0matrix_imtime(tau::Float64, q::Float64, param) """ Two particle scattering in a vacuum T-matrix in imaginary time """ β, me, as = param.β, param.me, param.as B = me^(3/2)/(4π) @@ -688,7 +688,7 @@ function T0maxtrix_imtime(tau::Float64, q::Float64, param) function integrand1(vars, tau, param) me, as = param.me, param.as m = me * as^2 - u = vars[1][1] + u = vars y = u/(1-u) fraction = y^2/(y^2 + tau/ m) return 2 * exp(-y^2) * fraction @@ -696,10 +696,10 @@ function T0maxtrix_imtime(tau::Float64, q::Float64, param) function integrand2(vars, q, tau, param) β, me = param.β, param.me - u = vars[1][1] + u = vars y = u/(1-u) factor = 1/(exp(β * (y^2/ tau + q^2/ (4 * me)))-1) - return factor * integrand1(vars, q, tau, param) + return factor * integrand1(vars, tau, param) end tgrid = CompositeGrid.LogDensedGrid( @@ -710,13 +710,15 @@ function T0maxtrix_imtime(tau::Float64, q::Float64, param) 0.005, # minimum interval length of log grid 5 # N of bottom layer ) - data = [(integrand1(t, q, tau, param)+integrand2(t, q, tau, param)) for (ti, t) in enumerate(tgrid.grid)] + data = [(integrand1(t, tau, param)) for (ti, t) in enumerate(tgrid.grid)] + # data = [0.0 for (ti,t) in enumerate(tgrid.grid)] int_result = Interp.integrate1D(data, tgrid) integralfactor = 1/(B * pi * sqrt(tau)) # Propagator factor propfactor = exp(-q^2/ (4 * me) * tau) + # return propfactor*(poleterm + int_result * integralfactor) return propfactor*(poleterm + int_result * integralfactor) end diff --git a/src/parameter.jl b/src/parameter.jl index c5634c3..37f5416 100644 --- a/src/parameter.jl +++ b/src/parameter.jl @@ -24,7 +24,7 @@ using ..Parameters EF::Float64 = 1.0 #kF^2 / (2me) β::Float64 = 200 # bare inverse temperature μ::Float64 = 1.0 - as::Float64 = -2.5 # s-wave scattering length (local contact interaction) + as::Float64 = 2.5 # s-wave scattering length (local contact interaction) # artificial parameters Λs::Float64 = 0.0 # Yukawa-type spin-symmetric interaction ~1/(q^2+Λs) diff --git a/src/polarization.jl b/src/polarization.jl index f99b345..bc6dca5 100644 --- a/src/polarization.jl +++ b/src/polarization.jl @@ -186,7 +186,8 @@ function Ladder0_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6 integrand = zeros(ComplexF64, kgrid.size) if dim == 3 for (ki, k) in enumerate(kgrid.grid) - if Ωn <1000 + #integrand[ki] = _LadderT3d_integrand(k, q, Ωn, param) + if Ωn <100 integrand[ki] = _LadderT3d_integrand(k, q, Ωn, param) @assert !isnan(integrand[ki]) "nan at k=$k, q=$q, n=$n" else @@ -195,7 +196,7 @@ function Ladder0_FiniteTemp(q::Float64, n::Int, param; scaleN=20, minterval=1e-6 end end result = Interp.integrate1D(integrand, kgrid) - if Ωn> 1000 + if Ωn> 100 result = -me^(3 / 2) / (4π) * sqrt(Complex(q^2 / (4me) - 1im * Ωn - 2μ)) # Replace high frequency tail with two-particle particle-particle term end return result @@ -207,7 +208,7 @@ function Ladder0_FreeElectron(q::Float64, n::Int, param) error("No support for finite-temperature polarization in $dim dimension!") end - return -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * n / β) + return -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * n / β - 2μ) end """ @@ -354,6 +355,7 @@ function Ladder0_FiniteTemp(q::Float64, n::AbstractVector, param; scaleN=20, min for (ki, k) in enumerate(kgrid.grid) for (mi, m) in enumerate(n) Ωm = 2π * m / β + #integrand[ki, mi] = _LadderT3d_integrand(k, q, Ωm, param) if Ωm < 1000 integrand[ki, mi] = _LadderT3d_integrand(k, q, Ωm, param) @assert !isnan(integrand[ki, mi]) "nan at k=$k, q=$q" @@ -366,7 +368,7 @@ function Ladder0_FiniteTemp(q::Float64, n::AbstractVector, param; scaleN=20, min result = Interp.integrate1D(integrand, kgrid; axis=1) for (mi, m) in enumerate(n) Ωm = 2π * m / β - if Ωm > 1000 + if Ωm > 100 result[mi] = -me^(3 / 2) / (4π) * sqrt(Complex(q^2 / (4me) - 1im * Ωm - 2μ)) # Replace with two-particle particle-particle high frequency tail end end @@ -379,11 +381,11 @@ function Ladder0_FreeElectron(q::Float64, n::AbstractVector, param) if dim != 3 error("No support for finite-temperature polarization in $dim dimension!") end - integrand = zeros(ComplexF64, length(n)) + result = zeros(ComplexF64, length(n)) for (mi, m) in enumerate(n) - integrand[mi] = Ladder0_FreeElectron(q, m, param) + result[mi] = -me^(3 / 2) / (4π) * sqrt(q^2 / (4me) - 1im * 2π * m / β - 2μ) end - return integrand + return result end diff --git a/test/galileaninvariance.jl b/test/galileaninvariance.jl new file mode 100644 index 0000000..126f777 --- /dev/null +++ b/test/galileaninvariance.jl @@ -0,0 +1,111 @@ +# using ElectronGas +using Parameters +using GreenFunc +using CompositeGrids +# print(dirname(pwd())) +include("./ElectronGas.jl/src/convention.jl") +include("./ElectronGas.jl/src/parameter.jl") +include("./ElectronGas.jl/src/polarization.jl") +include("./ElectronGas.jl/src/interaction.jl") +include("./ElectronGas.jl/src/twopoint.jl") +using Plots + +para = Parameter.rydbergUnit(0.01, 1, 3) +Ωn_mesh_ph = GreenFunc.ImFreq(para.β, BOSON; Euv=100.0, symmetry=:ph) +sgrid = 0.0 + +δladderarray_n = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64) +δTarray_n2 = GreenFunc.MeshArray(Ωn_mesh_ph, sgrid; dtype=ComplexF64) + +for (ki, k) in enumerate(sgrid) + for (ni, n) in enumerate(Ωn_mesh_ph.grid) + δTarray_n2[ni, ki] = Interaction.Tmatrix(k, n, para) - Interaction.T0matrix(k, n, para) + end +end + + + + +# Many Body Particle-Particle +mbladderarray = Polarization.Ladder0_FiniteTemp(0.0, Ωn_mesh_ph.grid, para) + +# Vacuum Particle-Particle +vacladderarray = Polarization.Ladder0_FreeElectron(0.0, Ωn_mesh_ph.grid, para) + +for (ki, k) in enumerate(sgrid) + for (ni, n) in enumerate(Ωn_mesh_ph.grid) + δladderarray_n[ni, ki] = mbladderarray[ni,ki] - vacladderarray[ni,ki] + end +end + +# for (ki, k) in enumerate(sgrid) +# for (ni, n) in enumerate(Ωn_mesh_ph.grid) +# δTarray_n[ni, ki] = Interaction.Tmatrix(k, n, para) - Interaction.T0matrix(k, n, para) +# end +# end + +δladderarray_τ = δladderarray_n |> to_dlr |> to_imtime +# δTarray_τ = δTarray_n |> to_dlr |> to_imtime + +# Analytic expression for T0(k, τ) # + +function T0_τ(q::Float64, τ::Float64, param) + + me, as, μ, β = param.me, param.as, param.μ, param.β + B = me^(3/2) / (4 * π) + # Generating a log densed composite grid with LogDensedGrid() + ygrid = CompositeGrid.LogDensedGrid(:gauss,[0.0,50.],[0.0],5,0.005,5) + + # Functions to be integrated + f(y) = ( exp(- y^2 ) * y^2 ) / (τ / (me * as) + y^2) + g(y) = ( exp(- y^2 ) * y^2 ) / (τ / (me * as) + y^2) * 1. / ( exp( β * ( y^2/ τ + (q^2 / (4me) - 2μ)) ) - 1) + + # Generate data for integration + fdata = [f(y) for (yi, y) in enumerate(ygrid.grid)] + gdata = [g(y) for (yi, y) in enumerate(ygrid.grid)] + + fint_result = Interp.integrate1D(fdata, ygrid) + gint_result = Interp.integrate1D(gdata, ygrid) + + return exp(-(q^2 / (4 * me) - 2μ) * τ) * (2/( B* π * sqrt(τ) ) * ( fint_result + gint_result ) - exp( τ / (me*as^2) ) / ( 1 - exp( β* ( 1 / (me*as^2) - (q^2 / (4me) - 2μ)) ) )) + # return - exp( τ / (me*as^2) ) / ( 1 - exp( β* ( 1 / (me*as^2) - (q^2 / (4me) - 2μ)) ) ) + +end + +function Σ_τ(q::Float64, param) + # Output self energy for the mesh imaginary values + + kgrid = Polarization.finitetemp_kgrid_ladder(param.μ, param.me, q, param.kF, 20, 1e-6, 10) + δTarray_n = GreenFunc.MeshArray(Ωn_mesh_ph, kgrid; dtype=ComplexF64) + + for (ki, k) in enumerate(kgrid) + for (ni, n) in enumerate(Ωn_mesh_ph.grid) + δTarray_n[ni, ki] = Interaction.Tmatrix(k+q, n, param) - Interaction.T0matrix(k+q, n, param) + end + end + + δTarray_τ = δTarray_n |> to_dlr |> to_imtime + integrand = GreenFunc.MeshArray(δTarray_τ, kgrid; dtype=ComplexF64) + result = GreenFunc.MeshArray(δTarray_τ.mesh[1]; dtype=ComplexF64) + + for (ti, t) in enumerate(δTarray_τ.mesh[1][:]) + for (ki, k) in enumerate(kgrid.grid) + integrand[ti, ki] = (T0_τ(k + q, t, param) + δTarray_τ[ti]) * exp(t * ( k^2 / (2*param.me) - param.μ)) * 1 / (exp( param.β * (k^2 / (2* param.me) - param.μ)) + 1) + end + result[ti] = Interp.integrate1D(integrand[ti, :], kgrid) + end + return result +end + +qarray = LinRange(0,2*para.kF, 5) +#test_τ = #[Σ_τ(qi, para) for qi in qarray] +#test_ω = test_τ |> to_dlr |> to_imfreq +#plot(test_ω.mesh[1][:], imag(test_ω)) + +#G0_τ = GreenFunc.MeshArray(test_τ.mesh[1]; dtype=ComplexF64) + +#for (ti, t) in enumerate(G0_τ.mesh[1][:]) +# G0_τ[ti] = T0_τ(0.0, t, para) #exp(t * ( 0.0^2 / (2*para.me) - para.μ)) * 1 / (exp( para.β * (0.0^2 / (2* para.me) - para.μ)) + 1) + +#end +#G0_ω = G0_τ |> to_dlr |> to_imfreq \ No newline at end of file