Ver4 ad xc

example/ver4/test_PP.jl

@@ -8,7 +8,6 @@ using ElectronLiquid
 using ElectronLiquid.CompositeGrids
 using ElectronLiquid.UEG
 
-
 function compare(data, expect, ratio=5)
     # println(data, ", ", expect)
     @test isapprox(data.val, expect, atol=ratio * data.err)
@@ -27,28 +26,54 @@ function PP_interaction_dynamic(n, para::ParaMC, kamp=para.kF, kamp2=para.kF; ka
     return Interp.integrate1D(Wp, xgrid)
 end
 
+function yukawa_pp(para)
+    factor = para.e0^2 / para.ϵ0 * para.NF / 2
+    kF2 = para.kF^2
+    return factor / 2 / kF2 * log((para.mass2 + 4 * kF2) / para.mass2)
+end
+
 @testset "PP" begin
     seed = 1234
     # p = (1, 0, 0)
-    p = (2, 0, 0)
-    rs = 2.0
-    beta = 25
-    mass2 = 1e-8
-    neval = 1e7
-    para = ElectronLiquid.ParaMC(rs=rs, beta=beta, Fs=0.0, order=2, mass2=mass2, isDynamic=true)
-    UEG.MCinitialize!(para)
+    order = 4
+    # p = (order, 0, 0)
+    # p = (1, 0, 0)
+    rs = 1.0
+    beta = 100
+    mass2 = 3.5
+    neval = 1e5
+    # para = ElectronLiquid.ParaMC(rs=rs, beta=beta, Fs=0.0, order=order, mass2=mass2, isDynamic=true)
+    para = ElectronLiquid.ParaMC(rs=rs, beta=beta, Fs=0.0, order=order, mass2=mass2, isDynamic=false)
+    # partition = UEG.partition(para.order)
+    # partition = [(2, 0, 0), (1, 0, 1), (1, 0, 0)]
+    # partition = [(1, 0, 0), (2, 0, 0), (3, 0, 0),]
+    # partition = [(order, 0, 0),]
+    partition = [(2, 2, 0),]
+    # partition = [(2, 0, 0), (1, 0, 1)]
+    # partition = [(2, 0, 0), (3, 0, 0)]
+    # filter = [NoHartree, NoBubble]
+    filter = [NoHartree,]
+    # UEG.MCinitialize!(para)
+    println(partition)
     println(para)
-    # diagram = Ver4.diagram(para, [p,]; channel=[], filter=[])
-    # diagram = Ver4.diagram(para, [p, (2, 0, 0)]; channel=[PPr,], filter=[NoFock, NoBubble])
-    diagram = Ver4.diagram(para, [p,]; channel=[PPr,], filter=[NoFock, NoBubble])
 
     ############################ generic PH one-angle average ###########################
-    nlist = [0, 1, 2]
-    paras = [Ver4.OneAngleAveraged(para, [para.kF, para.kF], [[0, nlist[1], -1], [0, nlist[2], -1], [0, nlist[3], -1]], :PP, 0),]
+    # nlist = [0, 1, 2]
+    # paras = [Ver4.OneAngleAveraged(para, [para.kF, para.kF], [[0, nlist[1], -1], [0, nlist[2], -1], [0, nlist[3], -1]], :PP, 0),]
+    paras = [Ver4.OneAngleAveraged(para, [para.kF, para.kF], [[0, 0, -1],], :PP, 0),]
+
+    # diagram = Ver4.diagram(para, [p, (2, 0, 0)]; channel=[PPr,], filter=[NoFock, NoBubble])
+    # diagram = Ver4.diagramParquet(para, [p,]; channel=[PHr, PHEr, PPr,], filter=[NoFock,])
+    # data, result = Ver4.one_angle_averaged_ParquetAD(paras, diagram; neval=neval, print=-1, seed=seed)
+
+    # diagram = Ver4.diagram(para, [p,]; channel=[PHr, PHEr, PPr,], filter=[NoHartree, NoBubble])
+    diagram = Ver4.diagram(para, partition; channel=[PHr, PHEr, PPr,], filter=filter)
     data, result = Ver4.one_angle_averaged(paras, diagram; neval=neval, print=-1, seed=seed)
-    # obs = data[p]
-    obs2 = data[(2, 0, 0)]
-    println(obs2)
+    println([data[p] for p in partition])
+    # obs2 = data[p]
+    # println(obs2)
+
+    # println(yukawa_pp(para))
     # println("obs 1:", obs[:, 1, 1])
     # println("obs 2:", obs[:, 2, 1])
     # println("obs 3:", obs[:, 3, 1])
@@ -60,4 +85,14 @@ end
     #     println(real(obs[:, i, 1][2]), ", ", -PP_interaction_dynamic(nlist[i], para) / 2)
     #     # compare(real(obs[:, i, 1][2]), -PP_interaction_dynamic(nlist[i], para) / 2)
     # end
+
+    diagram = Ver4.diagramParquet(para, partition; channel=[PHr, PHEr, PPr,], filter=filter)
+    data, result = Ver4.one_angle_averaged_ParquetAD(paras, diagram; neval=neval, print=-1, seed=seed)
+    println([data[p] for p in partition])
+
+    diagram = Ver4.diagramParquet_load(para, partition; filter=filter)
+    data, result = Ver4.one_angle_averaged_ParquetAD_Clib(paras, diagram; neval=neval, print=-1, seed=seed)
+    println([data[p] for p in partition])
+    # obs2 = data[p]
+    # println(obs2)
 end

example/ver4/test_PP_short.jl

@@ -0,0 +1,10 @@
+using ElectronLiquid
+order = 2
+neval = 1e3
+para = UEG.ParaMC(rs=1.0, beta=100, mass2=3.5, order=order, isDynamic=false)
+partition = [(o, 0, 0) for o in 1:order]
+# partition = UEG.partition(para.order)
+ver4, result = Ver4.MC_PP_ParquetAD_vegas(para; partition=partition, neval=neval, seed=1234)
+# ver4, result = Ver4.MC_PP_ParquetAD_Clib(para; partition=partition, neval=neval)
+# ver4, result = Ver4.MC_PP_ParquetAD(para; partition=partition, neval=neval)
+println(ver4)

example/ver4/ver4compile.jl

@@ -0,0 +1,29 @@
+using ElectronLiquid, FeynmanDiagram
+
+const order = 3
+const isDynamic = true
+
+if !isDynamic
+    para = UEG.ParaMC(rs=1.0, beta=25, order=order, isDynamic=false)
+    partition = UEG.partition(para.order)
+
+    println("generating diagrams")
+    diagram = Ver4.diagramParquet(para, partition; channel=[PHr, PHEr, PPr,], filter=[NoHartree,])
+    println("diagram generated")
+    partition, diagpara, FeynGraphs, extT_labels, spin_conventions = diagram
+
+    MaxLoopNum = maximum([key[1] for key in partition]) + 2
+    Ver4.compileC_ParquetAD_toFiles(para.order, partition, FeynGraphs, MaxLoopNum)
+else
+    para = UEG.ParaMC(rs=1.0, beta=25, order=order, isDynamic=true)
+    partition = UEG.partition(para.order)
+
+    println("generating diagrams")
+    diagram = Ver4.diagramParquet(para, partition; channel=[PHr, PHEr, PPr,], filter=[NoHartree, NoBubble])
+    # diagram = Ver4.diagramParquet(para, partition; channel=channel = [PHr, PHEr, PPr,], filter=[NoHartree,])
+    println("diagram generated")
+    partition, diagpara, FeynGraphs, extT_labels, spin_conventions = diagram
+
+    MaxLoopNum = maximum([key[1] for key in partition]) + 2
+    Ver4.compileC_ParquetAD_toFiles_dynamic(para.order, partition, FeynGraphs, MaxLoopNum)
+end

example/ver4/ver4funcdict.jl

@@ -0,0 +1,32 @@
+using FeynmanDiagram, ElectronLiquid
+
+const isDynamic = true
+
+function func_dict_str(o)
+    str = ""
+    str *= "const evalfuncParquetAD_map = Dict(\n"
+    for order in 1:o
+        para = UEG.ParaMC(rs=1.0, beta=25, order=order, isDynamic=isDynamic)
+        partition = UEG.partition(para.order)
+        if isDynamic
+            diagram = Ver4.diagramParquetDynamic_load(para, partition)
+        else
+            diagram = Ver4.diagramParquet_load(para, partition)
+        end
+        _partition = diagram[1]
+        println(partition)
+        println(_partition)
+        for p in _partition
+            str *= "($order, $(p[1]), $(p[2]), $(p[3])) => eval_ver4O$(order)ParquetAD$(p[1])$(p[2])$(p[3])!, \n"
+        end
+    end
+    str *= ")\n"
+    return str
+end
+
+# fname = "./src/vertex4/source_codeParquetAD/func_dict_ParquetAD.jl"
+fname = "./src/vertex4/source_codeParquetAD/dynamic/func_dict_ParquetADDynamic.jl"
+
+open(fname, "w") do f
+    write(f, func_dict_str(3))
+end

src/common/counterterm.jl

@@ -113,7 +113,7 @@ By definition, the chemical potential renormalization is defined as
 function chemicalpotential_renormalization(order, data, δμ; offset::Int=0)
     # _partition = sort([k for k in keys(rdata)])
     # println(_partition)
-    @assert order <= 5 "Order $order hasn't been implemented!"
+    @assert order <= 6 "Order $order hasn't been implemented!"
     @assert length(δμ) + 1 >= order
     data = mergeInteraction(data)
     d = data
@@ -159,6 +159,26 @@ function chemicalpotential_renormalization(order, data, δμ; offset::Int=0)
             d[(2 + offset, 1)] .* δμ[3] +
             d[(1 + offset, 1)] .* δμ[4]
     end
+    if order >= 6
+        # Σ6 = Σ60 + Σ51*δμ1 + ...
+        z[6] =
+            d[(6 + offset, 0)] +
+            d[(5 + offset, 1)] .* δμ[1] +
+            d[(4 + offset, 2)] .* δμ[1] .^ 2 +
+            d[(3 + offset, 3)] .* δμ[1] .^ 3 +
+            d[(2 + offset, 4)] .* δμ[1] .^ 4 +
+            d[(1 + offset, 5)] .* δμ[1] .^ 5 +
+            d[(4 + offset, 1)] .* δμ[2] +
+            d[(3 + offset, 2)] .* 2 .* δμ[1] .* δμ[2] +
+            d[(2 + offset, 3)] .* 3 .* δμ[1] .^ 2 .* δμ[2] +
+            d[(1 + offset, 4)] .* 4 .* δμ[1] .^ 3 .* δμ[2] +
+            d[(3 + offset, 1)] .* δμ[3] +
+            d[(2 + offset, 2)] .* (δμ[2] .^ 2 + 2 * δμ[1] .* δμ[3]) +
+            d[(1 + offset, 3)] .* 5 .* (δμ[2] .^ 2 .* δμ[1] + δμ[1] .^ 2 .* δμ[3]) +
+            d[(2 + offset, 1)] .* δμ[4] +
+            d[(1 + offset, 2)] .* 2 .* (δμ[2] .* δμ[3] + δμ[1] .* δμ[4]) +
+            d[(1 + offset, 1)] .* δμ[5]
+    end
     return z
 end
 
@@ -308,6 +328,9 @@ function _inverse(z::AbstractVector{T}) where {T}
         zi[5] = -z[1] .^ 5 + 4z[1] .^ 3 .* z[2] - 3z[1] .* z[2] .^ 2 - 3z[1] .^ 2 .* z[3] + 2z[2] .* z[3] + 2z[1] .* z[4] - z[5]
     end
     if order >= 6
+        zi[6] = z[1] .^ 6 - 5z[1] .^ 4 .* z[2] + (4z[1] .^ 3 .* z[3] + 6z[1] .^ 2 .* z[2] .^ 2) - (3z[1] .^ 2 .* z[4] + 6z[1] .* z[2] .* z[3]) + 2z[1] .* z[5] - z[6]
+    end
+    if order >= 7
         error("order must be <= 5")
     end
     return zi

src/common/eval.jl

@@ -9,6 +9,96 @@ using ..ElectronGas
 
 const TAU_CUTOFF = 1e-10
 
+function green_derive_diagtree(τ, ϵ, β, order)
+    if order == 0
+        result = green(τ, ϵ, β)
+    elseif order == 1
+        result = -Spectral.kernelFermiT_dω(τ, ϵ, β)
+    elseif order == 2
+        result = Spectral.kernelFermiT_dω2(τ, ϵ, β)
+    elseif order == 3
+        result = -Spectral.kernelFermiT_dω3(τ, ϵ, β)
+    elseif order == 4
+        result = Spectral.kernelFermiT_dω4(τ, ϵ, β)
+    elseif order == 5
+        result = -Spectral.kernelFermiT_dω5(τ, ϵ, β)
+    else
+        error("not implemented!")
+        # result = Propagator.green(τ, ϵ, β) * 0.0
+    end
+    return result
+end
+
+function green_derive(τ, ϵ, β, order)
+    return green_derive_diagtree(τ, ϵ, β, order) / factorial(order)
+end
+
+function interaction_derive(τ1, τ2, K, p::ParaMC, idorder; idtype=Instant, tau_num=1, isLayered=false)
+    dim, e0, ϵ0, mass2 = p.dim, p.e0, p.ϵ0, p.mass2
+    qd = sqrt(dot(K, K))
+    if idorder[2] == 0
+        if idtype == Instant
+            if tau_num == 1
+                # return e0^2 / ϵ0 / (dot(K, K) + mass2)
+                if !(isLayered)
+                    return Coulombinstant(qd, p)
+                else
+                    q = sqrt(dot(K, K) + 1e-16)
+                    invK = 1.0 / q
+                    return e0^2 / 2ϵ0 * invK * tanh(mass2 * q)
+                end
+            elseif tau_num == 2
+                # println(id.extT)
+                if idorder[3] == 0
+                    return interactionStatic(p, qd, τ1, τ2)
+                else # return dR/df for the RG purpose. The static part is zero
+                    return 0.0
+                end
+            else
+                error("not implemented!")
+            end
+        elseif idtype == Dynamic
+            if idorder[3] == 0
+                return interactionDynamic(p, qd, τ1, τ2)
+            else # return dR/df for the RG purpose.
+                return UEG.interactionDynamic_df(p, qd, τ1, τ2)
+            end
+        else
+            error("not implemented!")
+        end
+    else # counterterm for the interaction
+        order = idorder[2]
+        if idtype == Instant
+            if tau_num == 1
+                if dim == 3
+                    invK = 1.0 / (qd^2 + mass2)
+                    return e0^2 / ϵ0 * invK * (mass2 * invK)^order
+                elseif dim == 2
+                    if !(isLayered)
+                        invK = 1.0 / (qd + mass2)
+                        return e0^2 / 2ϵ0 * invK * (mass2 * invK)^order
+                    else
+                        return 0.0
+                    end
+                else
+                    error("not implemented!")
+                end
+            else
+                # return counterR(qd, τ1, τ2, idorder[2])
+                return 0.0 #for dynamical interaction, the counter-interaction is always dynamic!
+            end
+        elseif idtype == Dynamic
+            if idorder[3] == 0
+                return counterR(p, qd, τ1, τ2, idorder[2])
+            else
+                return counterR_df(p, qd, τ1, τ2, idorder[2])
+            end
+        else
+            error("not implemented!")
+        end
+    end
+end
+
 function green(τ::T, ω::T, β::T) where {T}
     #generate green function of fermion
     if τ ≈ T(0.0)
@@ -80,19 +170,24 @@ function green3(Ek, τ, beta=β)
     return green
 end
 
-##################### propagator and interaction evaluation ##############
-function DiagTree.eval(id::BareGreenId, K, extT, varT, p::ParaMC)
-    kF, β, me, μ, massratio = p.kF, p.β, p.me, p.μ, p.massratio
-    τin, τout = varT[id.extT[1]], varT[id.extT[2]]
-    k = norm(K)
+function dispersion(k, p::ParaMC)
+    kF, me, μ, massratio = p.kF, p.me, p.μ, p.massratio
     if p.isFock
         fock = SelfEnergy.Fock0_ZeroTemp(k, p.basic) - SelfEnergy.Fock0_ZeroTemp(kF, p.basic)
         ϵ = k^2 / (2me * massratio) - μ + fock
     else
         ϵ = k^2 / (2me * massratio) - μ
         # ϵ = kF / me * (k - kF)
     end
+    return ϵ
+end
 
+##################### propagator and interaction evaluation ##############
+function DiagTree.eval(id::BareGreenId, K, extT, varT, p::ParaMC)
+    β = p.β
+    τin, τout = varT[id.extT[1]], varT[id.extT[2]]
+    k = norm(K)
+    ϵ = dispersion(k, p)
     # if k < 0.4 * kF || k > kF * 1.3
     #     return 0.0
     # end
@@ -101,28 +196,7 @@ function DiagTree.eval(id::BareGreenId, K, extT, varT, p::ParaMC)
     # ∂^n_μ g(ϵₖ - μ, τ) = (-1)^n ∂^n_ω g(ω, τ)
     τ = τout - τin
     order = id.order[1]
-    if order == 0  # g[μ]
-        if τ ≈ 0.0
-            return Spectral.kernelFermiT(-TAU_CUTOFF, ϵ, β)
-        else
-            return Spectral.kernelFermiT(τ, ϵ, β)
-        end
-    elseif order == 1  # ∂_μ g[μ]
-        return -Spectral.kernelFermiT_dω(τ, ϵ, β)
-    elseif order == 2  # ∂^2_μ g[μ]
-        # return Spectral.kernelFermiT_dω2(τ, ϵ, β) / 2.0
-        return Spectral.kernelFermiT_dω2(τ, ϵ, β)
-    elseif order == 3  # ∂^3_μ g[μ]
-        # return -Spectral.kernelFermiT_dω3(τ, ϵ, β) / 6.0
-        return -Spectral.kernelFermiT_dω3(τ, ϵ, β)
-        # return 0.0
-    elseif order == 4
-        return Spectral.kernelFermiT_dω4(τ, ϵ, β)
-    elseif order == 5
-        return -Spectral.kernelFermiT_dω5(τ, ϵ, β)
-    else
-        error("not implemented!")
-    end
+    return green_derive_diagtree(τ, ϵ, β, order)
 end
 
 # eval(id::InteractionId, K, varT) = e0^2 / ϵ0 / (dot(K, K) + mass2)
@@ -175,9 +249,9 @@ function DiagTree.eval(id::BareInteractionId, K, extT, varT, p::ParaMC)
             end
         elseif id.type == Dynamic
             if id.order[3] == 0
-                return counterR(p, qd, varT[id.extT[1]], varT[id.extT[2]], id.order[2])
+                return factorial(order) * counterR(p, qd, varT[id.extT[1]], varT[id.extT[2]], id.order[2])
             else
-                return counterR_df(p, qd, varT[id.extT[1]], varT[id.extT[2]], id.order[2])
+                return factorial(order) * counterR_df(p, qd, varT[id.extT[1]], varT[id.extT[2]], id.order[2])
             end
         else
             error("not implemented!")