Liam

src/interaction.jl

@@ -525,4 +525,201 @@ function RPA_total(q, n, param; pifunc=Polarization0_ZeroTemp, Vinv_Bare=coulomb
     return Ws, Wa
 end
 
+"""
+    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.
+
+#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 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...))
+    else
+        return 1 ./ (1 / as .+ ladderfunc(q, n, param; kwargs...))
+    end
+end
+
+"""
+    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 
+(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 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
+    as = 4π * param.as / param.me
+    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_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_r, green_dyn_i, green_tail_r, green_tail_i
+end
+
+
+"""
+    function T0(q::Float64, n::Int, param; ; ladderfunc=Polarization.Ladder0_FreeElectron, kwargs...)
+
+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; 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 .+ ladderfunc(q, n, param; kwargs...) )
+    end
+end
+
+"""
+    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)
+"""
+
+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_imtime(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 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π)
+    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, tau, param)
+        me, as = param.me, param.as
+        m = me * as^2
+        u = vars
+        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
+        y = u/(1-u)
+        factor = 1/(exp(β * (y^2/ tau + q^2/ (4 * me)))-1)
+        return factor * integrand1(vars, 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, 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
+
 end

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,91 @@ 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_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_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
+    @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 -> 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))
+                kernel[ki, pi, ni] = (Hp - Hm)
+            end
+        end
+    end
+
+    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

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 = 2.5   # 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