Merge pull request #29 from numericalEFT/xc-from-houpc

Xc from houpc

example/ver4/plot_test_PP.py

@@ -0,0 +1,12 @@
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+mass2 = np.array([1e-2,1e-4,1e-6,1e-8])
+val = np.array([1.179, 1.379, 1.422, 1.422])
+err = np.array([0.003,0.005,0.007,0.008])
+
+plt.figure()
+plt.xscale("log")
+plt.errorbar(mass2, val, yerr=err, fmt="o-", markersize=4)
+plt.savefig("testPP.pdf")

example/ver4/test_PP.jl

@@ -0,0 +1,63 @@
+using Test
+using ElectronLiquid
+using FeynmanDiagram
+using FiniteDifferences
+using Lehmann
+using Measurements
+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)
+end
+
+function PP_interaction_dynamic(n, para::ParaMC, kamp=para.kF, kamp2=para.kF; kawargs...)
+    kF = para.kF
+    xgrid = CompositeGrid.LogDensedGrid(:gauss, [-1.0, 1.0], [-1.0, 1.0], 16, 0.001, 16)
+    qs = [sqrt(kamp^2 + kamp2^2 - 2 * x * kamp * kamp2) for x in xgrid]
+
+    Wp = zeros(Float64, length(qs))
+    for (qi, q) in enumerate(qs)
+        Wp[qi] = UEG.KO_W(q, n, para)
+    end
+    Wp *= para.NFstar  # additional minus sign because the interaction is exchanged
+    return Interp.integrate1D(Wp, xgrid)
+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)
+    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),]
+    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("obs 1:", obs[:, 1, 1])
+    # println("obs 2:", obs[:, 2, 1])
+    # println("obs 3:", obs[:, 3, 1])
+
+    # println(PP_interaction_dynamic(nlist[1], para) / 2)
+    # println(PP_interaction_dynamic(nlist[2], para) / 2)
+    # println(PP_interaction_dynamic(nlist[3], para) / 2)
+    # for i in 1:length(nlist)
+    #     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
+end

src/common/eval.jl

@@ -7,10 +7,12 @@ using ..Lehmann
 using LinearAlgebra
 using ..ElectronGas
 
+const TAU_CUTOFF = 1e-10
+
 function green(τ::T, ω::T, β::T) where {T}
     #generate green function of fermion
     if τ ≈ T(0.0)
-        τ = -1e-10
+        τ = -TAU_CUTOFF
     end
     if τ > T(0.0)
         return ω > T(0.0) ?
@@ -25,7 +27,7 @@ end
 
 function green2(Ek, τ, beta)
     if τ ≈ 0.0
-        τ = -1.0e-10
+        τ = -TAU_CUTOFF
     end
 
     s = 1.0
@@ -54,7 +56,7 @@ end
 
 function green3(Ek, τ, beta=β)
     if τ ≈ 0.0
-        τ = -1.0e-10
+        τ = -TAU_CUTOFF
     end
 
     s = 1.0
@@ -101,7 +103,7 @@ function DiagTree.eval(id::BareGreenId, K, extT, varT, p::ParaMC)
     order = id.order[1]
     if order == 0  # g[μ]
         if τ ≈ 0.0
-            return Spectral.kernelFermiT(-1e-8, ϵ, β)
+            return Spectral.kernelFermiT(-TAU_CUTOFF, ϵ, β)
         else
             return Spectral.kernelFermiT(τ, ϵ, β)
         end

src/common/interaction.jl

@@ -5,6 +5,7 @@ using ..Lehmann
 
 export Coulombinstant, KOinstant, KOstatic, interactionDynamic, interactionStatic, counterR
 
+const Q_CUTOFF = 1e-10
 
 """
     function lindhard(x)
@@ -61,8 +62,8 @@ end
 
 #fast version
 @inline function KOinstant(q, e0, dim, mass2, fs, kF)
-    if abs(q) < 1e-6
-        q = 1e-6
+    if abs(q) < Q_CUTOFF
+        q = Q_CUTOFF
     end
     if dim == 3
         return 4π * e0^2 / (q^2 + mass2) + fs
@@ -140,17 +141,17 @@ Rq = r_q / (1 - r_q Π0) - f
 ```
 """
 function KO_W(q, n::Integer, para::ParaMC; Pi=polarKW(q, n, para))
-    if abs(q) < 1e-6
-        q = 1e-6
+    if abs(q) < Q_CUTOFF
+        q = Q_CUTOFF
     end
     invKOinstant = 1.0 / KOinstant(q, para)
     Rs = 1.0 ./ (invKOinstant - Pi) - para.fs
     return Rs
 end
 
 function KO_W_df(q, n::Integer, para::ParaMC; Pi=polarKW(q, n, para))
-    if abs(q) < 1e-6
-        q = 1e-6
+    if abs(q) < Q_CUTOFF
+        q = Q_CUTOFF
     end
     # invKOinstant = 1.0 / KOinstant(q, para)
     # Rs = 1.0 ./ (invKOinstant - Pi) - para.fs
@@ -406,9 +407,9 @@ function counterR(p::ParaMC, qd, τIn, τOut, order)
 
     # if qd <= qgrid.grid[1]
     # the current interpolation vanishes at q=0, which needs to be corrected!
-    if qd <= 1e-6 * kF
+    if qd <= Q_CUTOFF * kF
         # q = qgrid.grid[1] + 1.0e-6
-        qd = 1e-6 * kF
+        qd = Q_CUTOFF * kF
     end
 
     if order <= p.order
@@ -429,9 +430,9 @@ function counterR_df(p::ParaMC, qd, τIn, τOut, order)
 
     # if qd <= qgrid.grid[1]
     # the current interpolation vanishes at q=0, which needs to be corrected!
-    if qd <= 1e-6 * kF
+    if qd <= Q_CUTOFF * kF
         # q = qgrid.grid[1] + 1.0e-6
-        qd = 1e-6 * kF
+        qd = Q_CUTOFF * kF
     end
 
     if order <= p.order
@@ -471,8 +472,8 @@ where Π₀ is the polarization of free electrons.
 
     # if qd <= qgrid.grid[1]
     # the current interpolation vanishes at q=0, which needs to be corrected!
-    if qd <= 1e-6 * kF
-        qd = 1e-6 * kF
+    if qd <= Q_CUTOFF * kF
+        qd = Q_CUTOFF * kF
     end
 
     vd = KOinstant(qd, p)
@@ -519,10 +520,10 @@ kostatic - dynamic = v_q
     kF, maxK = p.kF, p.maxK
 
     if qd > maxK
-        return (KOinstant(qd, p)-p.fs) / p.β #if qd is very large, the interactin is reduced to the bare one
+        return (KOinstant(qd, p) - p.fs) / p.β #if qd is very large, the interactin is reduced to the bare one
     end
-    if qd <= 1e-6 * kF
-        qd = 1e-6 * kF
+    if qd <= Q_CUTOFF * kF
+        qd = Q_CUTOFF * kF
     end
     # if there is no dynamic interactoin
     # return KOinstant(qd) - p.fs
@@ -570,8 +571,8 @@ where Π₀ is the polarization of free electrons.
 
     # if qd <= qgrid.grid[1]
     # the current interpolation vanishes at q=0, which needs to be corrected!
-    if qd <= 1e-6 * kF
-        qd = 1e-6 * kF
+    if qd <= Q_CUTOFF * kF
+        qd = Q_CUTOFF * kF
     end
 
     return linear2D(p.dW0_f, p.qgrid, p.τgrid, qd, dτ) # dynamic interaction, don't forget the singular factor vq

src/common/para_builder.jl

@@ -85,6 +85,8 @@ function MCinitialize!(para::ParaMC, bubble::Bool=true)
     !para.isDynamic && return
     para.dW0 .= KOdynamic_T(para)
     para.dW0_f .= KOdynamic_T_df(para)
+    para.cRs::Vector{Matrix{Float64}} = []
+    para.cRs_f::Vector{Matrix{Float64}} = []
     for o in 1:para.order-1
         push!(para.cRs, counterKO_T(para; order=o, bubble=bubble))
         push!(para.cRs_f, counterKO_T_df(para; order=o, bubble=bubble))

src/vertex3/ver3KW.jl

@@ -86,23 +86,21 @@ function KW(para::ParaMC, diagram;
     obs = [zeros(ComplexF64, 2, Nkin, Nqout, Nwin, Nwqout) for p in diagpara] # observable for the Fock diagram 
 
     if isnothing(config)
-        config = MCIntegration.Configuration(
-            var=(K, T, vKin, vQout, vWin, vWqout),
+        config = MCIntegration.Configuration(; var=(K, T, vKin, vQout, vWin, vWqout),
             dof=dof,
             obs=obs,
             type=Weight,
             userdata=(para, diag, root, extT, kin, qout, nkin, nqout),
-            kwargs...
-        )
+            kwargs...)
     end
     result = integrate(integrandKW; config=config, measure=measureKW, solver=:mcmc, neval=neval, print=print, kwargs...)
 
     if isnothing(result) == false
-        if print >= 0
-            report(result.config)
-            report(result; pick=o -> (real(o[1, 1, 1, 1, 1])), name="uu")
-            report(result; pick=o -> (real(o[2, 1, 1, 1, 1])), name="ud")
-        end
+        # if print >= 0
+        #     report(result.config)
+        #     report(result; pick=o -> (real(o[1, 1, 1, 1, 1])), name="uu")
+        #     report(result; pick=o -> (real(o[2, 1, 1, 1, 1])), name="ud")
+        # end
 
         datadict = Dict{eltype(partition),Any}()
         for k in 1:length(dof)

src/vertex3/ver3angleavg.jl

@@ -0,0 +1,119 @@
+function integrandAA(idx, var, config)
+    para, diag, root, extT, kamp, kamp2, nkin, nqout = config.userdata
+
+    kF, β = para.kF, para.β
+    varK, varT, varX = var[1], var[2], var[7]
+    loopNum = config.dof[idx][1]
+    # error(loopNum)
+    _kin, _kout = kamp[var[3][1]], kamp2[var[4][1]]
+    _nkin, _nqout = nkin[var[5][1]], nqout[var[6][1]]
+
+    # println(kinL, ", ", koutL, ", ", kinR)
+    x = varX[1][1]
+    varK.data[1, 1] = _kout * x - _kin
+    varK.data[2, 1] = _kout * sqrt(1 - x^2)
+    varK.data[1, 2] = _kin
+
+    diagram = diag[idx]
+    weight = diagram.node.current
+    rootuu, rootud = root[1][idx], root[2][idx]
+    extTuu, extTud = extT[1][idx], extT[2][idx]
+
+    # varK.data[:, 2] = [kF * x, kF * sqrt(1 - x^2), 0.0]
+
+    ExprTree.evalKT!(diagram, varK.data, varT.data, para)
+    if !isempty(rootuu)
+        wuu = sum(weight[root] * phaseC(varT, extTuu[ri], _nkin, _nqout, β) for (ri, root) in enumerate(rootuu))
+    else
+        wuu = zero(ComplexF64)
+    end
+    if !isempty(rootud)
+        wud = sum(weight[root] * phaseC(varT, extTud[ri], _nkin, _nqout, β) for (ri, root) in enumerate(rootud))
+    else
+        wud = zero(ComplexF64)
+    end
+    # factor = para.NF / (2π)^(para.dim * loopNum)
+    factor = 1.0 / (2π)^(para.dim * loopNum)
+
+    # if isdefined(Main, :Infiltrator)
+    #     Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__)
+    # end
+
+    return Weight(wuu * factor, wud * factor)
+    # return Weight(zero(ComplexF64), wud * para.NF)
+end
+
+function measureAA(idx, var, obs, weight, config)
+    kin, kout = var[3][1], var[4][1]
+    nkin, nqout = var[5][1], var[6][1]
+    obs[idx][1, kin, kout, nkin, nqout] += weight.d
+    obs[idx][2, kin, kout, nkin, nqout] += weight.e
+end
+
+function AA(para::ParaMC, diagram;
+    kamp=[para.kF,], kamp1=[kamp[1],],
+    nkin=[0,],
+    nqout=[0,],
+    neval=1e6, #number of evaluations
+    print=0,
+    alpha=3.0, #learning ratio
+    config=nothing,
+    kwargs...
+)
+    UEG.MCinitialize!(para)
+
+    dim, β, kF, NF = para.dim, para.β, para.kF, para.NF
+
+    partition, diagpara, diag, root, extT = diagram
+    @assert length(diagpara) == length(diag) == length(root[1]) == length(extT[1])
+    @assert length(root[1]) == length(root[2])
+    @assert length(extT[1]) == length(extT[2])
+
+    K = MCIntegration.FermiK(para.dim, kF, 0.2 * kF, 10.0 * kF, offset=2)
+    T = MCIntegration.Continuous(0.0, β, offset=1, alpha=alpha)
+    T.data[1] = 0.0
+    X = MCIntegration.Continuous(-1.0, 1.0) #x=cos(θ)
+
+    Nkin, Nkout = length(kamp), length(kamp1)
+    Nwin, Nwqout = length(nkin), length(nqout)
+
+    vKin = MCIntegration.Discrete(1, Nkin, alpha=alpha)
+    vKout = MCIntegration.Discrete(1, Nkout, alpha=alpha)
+
+    vWin = MCIntegration.Discrete(1, Nwin, alpha=alpha)
+    vWqout = MCIntegration.Discrete(1, Nwqout, alpha=alpha)
+
+    dof = [[p.innerLoopNum, p.totalTauNum - 1, 1, 1, 1, 1, 1] for p in diagpara] # K, T, ExtKidx
+    obs = [zeros(ComplexF64, 2, Nkin, Nkout, Nwin, Nwqout) for p in diagpara] # observable for the Fock diagram 
+
+    if isnothing(config)
+        config = MCIntegration.Configuration(; var=(K, T, vKin, vKout, vWin, vWqout, X),
+            dof=dof,
+            obs=obs,
+            type=Weight,
+            userdata=(para, diag, root, extT, kamp, kamp1, nkin, nqout),
+            kwargs...)
+    end
+    result = integrate(integrandAA; config=config, measure=measureAA, solver=:mcmc, neval=neval, print=print, kwargs...)
+
+    if isnothing(result) == false
+        # if print >= 0
+        #     report(result.config)
+        #     report(result; pick=o -> (real(o[1, 1, 1, 1, 1])), name="uu")
+        #     report(result; pick=o -> (real(o[2, 1, 1, 1, 1])), name="ud")
+        # end
+
+        datadict = Dict{eltype(partition),Any}()
+        for k in 1:length(dof)
+            avg, std = result.mean[k], result.stdev[k]
+            r = measurement.(real(avg), real(std))
+            i = measurement.(imag(avg), imag(std))
+            data = Complex.(r, i)
+            datadict[partition[k]] = data
+        end
+        return datadict, result
+    else
+        return nothing, nothing
+    end
+
+end

src/vertex3/vertex3.jl

@@ -109,5 +109,6 @@ end
 # end
 
 include("ver3KW.jl")
+include("ver3angleavg.jl")
 
 end