Merge pull request #23 from numericalEFT/dev

Dev

Author: houpc

src/interaction.jl

@@ -61,7 +61,7 @@ end
 
 function bubbledyson(V::Float64, F::Float64, Π::Float64, n::Int)
     # V:bare interaction
-    # G:G^{+-} is local field factor,0 for RPA
+    # F:F^{+-} is local field factor,0 for RPA
     # Π:Polarization. 2*Polarization0 for spin 1/2
     # n:matfreq. special case for n=0
     # comparing to previous convention, an additional V is multiplied

src/parameter.jl

@@ -67,17 +67,17 @@ generate Para with a complete set of parameters, no value presumed.
  - EF: Fermi energy
  - β: inverse temperature
 """
-@inline function fullUnit(ϵ0, e0, me, EF, β, dim = 3, spin = 2)
-    return Para(
-        dim = dim,
+@inline function fullUnit(ϵ0, e0, me, EF, β, dim = 3, spin = 2; kwargs...)
+    para = Para(dim = dim,
         spin = spin,
         ϵ0 = ϵ0,
         e0 = e0,
         me = me,
         EF = EF,
         β = β,
-        μ = EF,
+        μ = EF
     )
+    return reconstruct(para, kwargs...)
 end
 
 """
@@ -89,13 +89,13 @@ assume 4πϵ0=1, me=0.5, EF=1
  - Θ: dimensionless temperature. Since EF=1 we have β=beta
  - rs: Wigner-Seitz radius over Bohr radius.
 """
-@inline function defaultUnit(Θ, rs, dim = 3, spin = 2)
+@inline function defaultUnit(Θ, rs, dim = 3, spin = 2; kwargs...)
     ϵ0 = 1 / (4π)
     e0 = (dim == 3) ? sqrt(2 * rs / (9π / (2spin))^(1 / 3)) : sqrt(sqrt(2) * rs)
     me = 0.5
     EF = 1
     β = 1 / Θ / EF
-    return fullUnit(ϵ0, e0, me, EF, β, dim, spin)
+    return fullUnit(ϵ0, e0, me, EF, β, dim, spin; kwargs...)
 end
 
 
@@ -108,14 +108,14 @@ assume 4πϵ0=1, me=0.5, e0=sqrt(2)
  - Θ: dimensionless temperature. beta could be different from β
  - rs: Wigner-Seitz radius over Bohr radius.
 """
-@inline function rydbergUnit(Θ, rs, dim = 3, spin = 2)
+@inline function rydbergUnit(Θ, rs, dim = 3, spin = 2; kwargs...)
     ϵ0 = 1 / (4π)
     e0 = sqrt(2)
     me = 0.5
     kF = (dim == 3) ? (9π / (2spin))^(1 / 3) / rs : sqrt(4 / spin) / rs
     EF = kF^2 / (2me)
     β = 1 / Θ / EF
-    return fullUnit(ϵ0, e0, me, EF, β, dim, spin)
+    return fullUnit(ϵ0, e0, me, EF, β, dim, spin; kwargs...)
 end
 
 export Para, Param

test/interaction.jl

@@ -3,17 +3,17 @@
     @testset "RPA and KO" begin
         beta = 1e4
         rs = 2.0
-        param = Interaction.Parameter.defaultUnit(1/beta,rs)
-        println(Interaction.RPA(0.01, 1, param))
-        println(Interaction.KO(0.01, 1, param))
-        println(Interaction.RPA(1.0, 1, param; pifunc = Interaction.Polarization.Polarization0_FiniteTemp))
-        println(Interaction.KO(1.0, 1, param; pifunc = Interaction.Polarization.Polarization0_FiniteTemp))
+        param = Interaction.Parameter.defaultUnit(1 / beta, rs)
+        # println(Interaction.RPA(0.01, 1, param))
+        # println(Interaction.KO(0.01, 1, param))
+        # println(Interaction.RPA(1.0, 1, param; pifunc = Interaction.Polarization.Polarization0_FiniteTemp))
+        # println(Interaction.KO(1.0, 1, param; pifunc = Interaction.Polarization.Polarization0_FiniteTemp))
 
         RPAs, RPAa = Interaction.RPAwrapped(100 * param.EF, 1e-8, [1e-8, 0.5, 1.0, 2.0, 10.0], param)
-        println(RPAs.dynamic)
-        println(RPAs.instant)
+        # println(RPAs.dynamic)
+        # println(RPAs.instant)
         KOs, KOa = Interaction.KOwrapped(100 * param.EF, 1e-8, [1e-8, 0.5, 1.0, 2.0, 10.0], param)
-        println(KOs.dynamic)
-        println(KOs.instant)
+        # println(KOs.dynamic)
+        # println(KOs.instant)
     end
 end

test/legendreinteraction.jl

@@ -1,51 +1,51 @@
 @testset "LegendreInteraction" begin
 
     # set parameters
-    param = LegendreInteraction.Parameter.defaultUnit(1/1000.0, 1.0)
+    param = LegendreInteraction.Parameter.defaultUnit(1 / 1000.0, 1.0)
     kF = param.kF
     Nk, minK, order, maxK = 16, 1e-8, 6, 10.0
-    print(param)
+    # print(param)
 
     @testset "Helper function" begin
         # test helper function
         # test with known result: bare coulomb interaction
-        H1=LegendreInteraction.helper_function(
-            1.0, 1, u->LegendreInteraction.interaction_instant(u,param,:sigma), param;
-            Nk=Nk,minK=minK,order=order
+        H1 = LegendreInteraction.helper_function(
+            1.0, 1, u -> LegendreInteraction.interaction_instant(u, param, :sigma), param;
+            Nk = Nk, minK = minK, order = order
         )
-        H2=LegendreInteraction.helper_function(
-            2.0, 1, u->LegendreInteraction.interaction_instant(u,param,:sigma), param;
-            Nk=Nk,minK=minK,order=order
+        H2 = LegendreInteraction.helper_function(
+            2.0, 1, u -> LegendreInteraction.interaction_instant(u, param, :sigma), param;
+            Nk = Nk, minK = minK, order = order
         )
         # println("$(H2-H1) == $(4*π*param.e0^2*log(2))")
-        @test isapprox(H2-H1, 4*π*param.e0^2*log(2), rtol = 1e-4)
+        @test isapprox(H2 - H1, 4 * π * param.e0^2 * log(2), rtol = 1e-4)
 
-        helper_grid = CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.1*maxK], [0.0, 2kF], 4, 0.001, 4)
-        intgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, helper_grid[end]], [0.0,2kF], 2Nk, 0.01minK, 2order)
-        helper = LegendreInteraction.helper_function_grid(helper_grid,intgrid, 1, u->LegendreInteraction.interaction_instant(u,param,:sigma),param)
-        helper_analytic = (4*π*param.e0^2) .* log.(helper_grid.grid)
+        helper_grid = CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.1 * maxK], [0.0, 2kF], 4, 0.001, 4)
+        intgrid = CompositeGrid.LogDensedGrid(:cheb, [0.0, helper_grid[end]], [0.0, 2kF], 2Nk, 0.01minK, 2order)
+        helper = LegendreInteraction.helper_function_grid(helper_grid, intgrid, 1, u -> LegendreInteraction.interaction_instant(u, param, :sigma), param)
+        helper_analytic = (4 * π * param.e0^2) .* log.(helper_grid.grid)
         helper_old = zeros(Float64, helper_grid.size)
         for (yi, y) in enumerate(helper_grid)
-            helper_old[yi] = LegendreInteraction.helper_function(y, 1, u->LegendreInteraction.interaction_instant(u,param,:sigma),param)
+            helper_old[yi] = LegendreInteraction.helper_function(y, 1, u -> LegendreInteraction.interaction_instant(u, param, :sigma), param)
         end
-        println(helper .- helper[1])
-        println(helper_old .- helper_old[1])
-        println(helper_analytic .- helper_analytic[1])
+        # println(helper .- helper[1])
+        # println(helper_old .- helper_old[1])
+        # println(helper_analytic .- helper_analytic[1])
         rtol = 1e-6
         for (yi, y) in enumerate(helper_grid)
-            @test isapprox(helper[yi]-helper[1], helper_analytic[yi]-helper_analytic[1], rtol=rtol)
-            @test isapprox(helper_old[yi]-helper_old[1], helper_analytic[yi]-helper_analytic[1], rtol=rtol)
+            @test isapprox(helper[yi] - helper[1], helper_analytic[yi] - helper_analytic[1], rtol = rtol)
+            @test isapprox(helper_old[yi] - helper_old[1], helper_analytic[yi] - helper_analytic[1], rtol = rtol)
         end
     end
 
     @testset "DCKernel" begin
         # test DCKernel
 
-        kernel = LegendreInteraction.DCKernel0(param; spin_state=:sigma, Nk=8, rtol=1e-10)
+        kernel = LegendreInteraction.DCKernel0(param; spin_state = :sigma, Nk = 8, rtol = 1e-10)
         # kernel = LegendreInteraction.DCKernel(param, 100*param.EF, 1e-8, 5, 10*param.kF, 1e-7*param.kF, 4, :rpa,0,:sigma)
         # kernel2 = LegendreInteraction.DCKernel(kernel, 500.0)
         kF_label = searchsortedfirst(kernel.kgrid.grid, kernel.param.kF)
         qF_label = searchsortedfirst(kernel.qgrids[kF_label].grid, kernel.param.kF)
-        println(kernel.kernel[kF_label,qF_label,:])
+        println(kernel.kernel[kF_label, qF_label, :])
     end
 end

test/polarization.jl

@@ -28,7 +28,7 @@
         println("n = $n,  ω_n = $(2π*n/param.β)")
         for (qi, q) in enumerate(qgrid)
             polar[qi] = Polarization.Polarization0_ZeroTemp(q, n, param)
-            println("q=$q, Π0=$(polar[qi])")
+            # println("q=$q, Π0=$(polar[qi])")
         end
     end
 
@@ -40,21 +40,21 @@
         println("q = $q")
         for (i, n) in enumerate(ngrid)
             polar[i] = Polarization.Polarization0_ZeroTemp(q, n, param)
-            println("n=$n, Π0=$(polar[i]) \n")
+            # println("n=$n, Π0=$(polar[i]) \n")
         end
     end
 
-    while true
-        print("Please enter a whole number between 1 and 5: ")
-        input = readline(stdin)
-        value = tryparse(Int, input)
-        if value !== nothing && 1 <= value <= 5
-            println("You entered $(input)")
-            break
-        else
-            @warn "Enter a whole number between 1 and 5"
-        end
-    end
+    # while true
+    #     print("Please enter a whole number between 1 and 5: ")
+    #     input = readline(stdin)
+    #     value = tryparse(Int, input)
+    #     if value !== nothing && 1 <= value <= 5
+    #         println("You entered $(input)")
+    #         break
+    #     else
+    #         @warn "Enter a whole number between 1 and 5"
+    #     end
+    # end
     # plt = plot(qgrid, polar, ytics = -0.08:0.01, ls = 1, Axes(grid = :on, key = "left"))
     # display(plt)
     # save(term = "png", output = "polar_2d.png",

test/selfenergy.jl

@@ -18,16 +18,16 @@
 
         ΣR = real(Σ.dynamic)
         ΣI = imag(Σ.dynamic)
-        println(Σ.instant[1, 1, :])
-        println(ΣR[1, 1, kF_label, :])
-        println(ΣI[1, 1, kF_label, :])
+        # println(Σ.instant[1, 1, :])
+        # println(ΣR[1, 1, kF_label, :])
+        # println(ΣI[1, 1, kF_label, :])
         Z0 = (SelfEnergy.zfactor(Σ))
         # @test isapprox(Z0, 0.862, rtol = 5e-3)
 
         println("θ = $θ,  rs= $rs")
         println("Z-factor = $Z0")
         G = SelfEnergy.Gwrapped(Σ, param)
-        println(G.dynamic[1, 1, kF_label, :])
+        # println(G.dynamic[1, 1, kF_label, :])
     end
     # @testset "default unit" begin
     #     # make sure everything works for different unit sets