Merge pull request #39 from numericalEFT/bugfix

Bugfix

Project.toml

@@ -1,7 +1,7 @@
 name = "MCIntegration"
 uuid = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167"
 authors = ["Kun Chen", "Xiansheng Cai", "Pengcheng Hou"]
-version = "0.3.3"
+version = "0.3.4"
 
 [deps]
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"

example/bubble.jl

@@ -1,89 +1,124 @@
 # This example demonstrated how to calculate the bubble diagram of free electrons using the Monte Carlo module
 
-using LinearAlgebra, Random, Printf, BenchmarkTools, InteractiveUtils, Parameters
-using ElectronGas
+using LinearAlgebra, Random, Printf
 using StaticArrays
-using Lehmann
 using MCIntegration
 # using ProfileView
 
-const Steps = 1e6
-
-# include("parameter.jl")
-beta = 25.0
-rs = 1.0
-const basic = Parameter.rydbergUnit(1 / beta, rs, 3)
-const β = basic.β
-const kF = basic.kF
-const me = basic.me
-const spin = basic.spin
-
-@with_kw struct Para
-    n::Int = 0 # external Matsubara frequency
-    Qsize::Int = 8
-    extQ::Vector{SVector{3,Float64}} = [@SVector [q, 0.0, 0.0] for q in LinRange(0.0 * kF, 2.0 * kF, Qsize)]
+Steps = 1e6
+
+Base.@kwdef struct Para
+    rs::Float64 = 1.0
+    beta::Float64 = 25.0
+    spin::Int = 2
+    Qsize::Int = 4
+    # n::Int = 0 # external Matsubara frequency
+    dim::Int = 3
+    me::Float64 = 0.5
+
+    kF::Float64 = (dim == 3) ? (9π / (2spin))^(1 / 3) / rs : sqrt(4 / spin) / rs
+    extQ::Vector{SVector{3,Float64}} = [@SVector [q, 0.0, 0.0] for q in LinRange(0.0 * kF, 1.5 * kF, Qsize)]
+    β::Float64 = beta / (kF^2 / 2me)
 end
 
-function integrand(T, K, Ext; userdata)
+function lindhard(q, para) #free electron polarization
+    me, kF, β, spin = para.me, para.kF, para.β, para.spin
+    density = me * kF / (2π^2)
+    # check sign of q, use -q if negative
+    (q < 1e-6) && (q = 1e-6)
+    x = q / 2 / kF
+    if abs(q - 2 * kF) > 1e-6
+        Π = (1 + (1 - x^2) * log1p(4 * x / ((1 - x)^2)) / 4 / x)
+    else
+        Π = 1.0
+    end
+    return -Π * density * spin / 2
+end
+
+function green(τ::T, ω::T, β::T) where {T}
+    if τ >= T(0.0)
+        return ω > T(0.0) ?
+               exp(-ω * τ) / (1 + exp(-ω * β)) :
+               exp(ω * (β - τ)) / (1 + exp(ω * β))
+    else
+        return ω > T(0.0) ?
+               -exp(-ω * (τ + β)) / (1 + exp(-ω * β)) :
+               -exp(-ω * τ) / (1 + exp(ω * β))
+    end
+end
+
+function integrand(vars, config) #for the vegas and vegasmc algorithms
     # @assert idx == 1 "$(idx) is not a valid integrand"
-    para, _Ext = userdata
-    k = K[1]
-    # Tin, Tout = T[1], T[2]
+    R, Theta, Phi, T, Ext = vars
+    para = config.userdata
+    kF, β, me = para.kF, para.β, para.me
+
+    r = R[1] / (1 - R[1])
+    θ = Theta[1]
+    ϕ = Phi[1]
+    # varK[:, i+1] .= [r * sin(θ) * cos(ϕ), r * sin(θ) * sin(ϕ), r * cos(θ)]
+    k = [r * sin(θ) * cos(ϕ), r * sin(θ) * sin(ϕ), r * cos(θ)]
+    factor = 1.0 / (2π)^(para.dim)  #each momentum loop is ∫dkxdkydkz/(2π)^3
+    factor *= r^2 / (1 - R[1])^2 * sin(θ)
+
     Tin, Tout = 0.0, T[1]
     extidx = Ext[1]
     q = para.extQ[extidx] # external momentum
     kq = k + q
     τ = (Tout - Tin)
     ω1 = (dot(k, k) - kF^2) / (2me)
-    g1 = Spectral.kernelFermiT(τ, ω1, β)
+    g1 = green(τ, ω1, β)
     ω2 = (dot(kq, kq) - kF^2) / (2me)
-    g2 = Spectral.kernelFermiT(-τ, ω2, β)
-    phase = 1.0 / (2π)^3
-    return g1 * g2 * spin * phase * cos(2π * para.n * τ / β)
+    g2 = green(-τ, ω2, β)
+    n = 0 # external Matsubara frequency
+    return g1 * g2 * para.spin * factor * cos(2π * n * τ / β)
 end
 
-function measure(obs, weight; userdata)
-    # @assert idx == 1 "$(idx) is not a valid integrand"
-    para, Ext = userdata
-    obs[Ext[1]] += weight[1]
+function integrand(idx, vars, config) #for the mcmc algorithm
+    return integrand(vars, config)::Float64
 end
 
-function run(steps)
+function measure(vars, obs, weight, config) # for vegas and vegasmc algorithms
+    Ext = vars[end]
+    obs[1][Ext[1]] += weight[1]
+end
+function measure(idx, vars, obs, weight, config) # for the mcmc algorithm
+    measure(vars, obs, weight, config)
+end
 
+function run(steps, alg)
     para = Para()
-    @unpack extQ, Qsize = para
+    extQ, Qsize = para.extQ, para.Qsize
+    kF, β = para.kF, para.β
 
-    # T = MCIntegration.Tau(β, β / 2.0)
-    T = MCIntegration.Continuous(0.0, β; alpha=3.0, adapt=true)
-    K = MCIntegration.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
-    Ext = MCIntegration.Discrete(1, length(extQ); adapt=true) # external variable is specified
+    T = Continuous(0.0, β; alpha=3.0, adapt=true)
+    R = Continuous(0.0, 1.0; alpha=3.0, adapt=true)
+    θ = Continuous(0.0, 1π; alpha=3.0, adapt=true)
+    ϕ = Continuous(0.0, 2π; alpha=3.0, adapt=true)
+    # K = MCIntegration.FermiK(3, kF, 0.2 * kF, 10.0 * kF)
+    Ext = Discrete(1, length(extQ); adapt=false) # external variable is specified
 
-    dof = [[1, 1, 1],] # degrees of freedom of the normalization diagram and the bubble
-    obs = zeros(Float64, Qsize) # observable for the normalization diagram and the bubble
+    dof = [[1, 1, 1, 1, 1],] # degrees of freedom of the normalization diagram and the bubble
+    obs = [zeros(Float64, Qsize),] # observable for the normalization diagram and the bubble
 
-    # config = MCIntegration.Configuration(var=(T, K, Ext), dof=dof, obs=obs, para=para)
-    result = MCIntegration.integrate(integrand; measure=measure, userdata=(para, Ext),
-        var=(T, K, Ext), dof=dof, obs=obs, solver=:vegas,
-        neval=steps, print=0, block=16)
+    result = integrate(integrand; measure=measure, userdata=para,
+        var=(R, θ, ϕ, T, Ext), dof=dof, obs=obs, solver=alg,
+        neval=steps, print=0, debug=true)
 
     if isnothing(result) == false
-        @unpack n, extQ = Para()
         avg, std = result.mean, result.stdev
 
+        println("Algorithm : $(alg)")
         @printf("%10s  %10s   %10s  %10s\n", "q/kF", "avg", "err", "exact")
         for (idx, q) in enumerate(extQ)
             q = q[1]
-            p = Polarization.Polarization0_ZeroTemp(q, para.n, basic) * spin
-            @printf("%10.6f  %10.6f ± %10.6f  %10.6f\n", q / basic.kF, avg[idx], std[idx], p)
+            p = lindhard(q, para)
+            @printf("%10.6f  %10.6f ± %10.6f  %10.6f\n", q / kF, avg[idx], std[idx], p)
         end
-        # println(MCIntegration.summary(result))
-        # i = 1
-        # println(result.config.var[i].histogram)
-        # println(sum(result.config.var[i].histogram))
-        # println(result.config.var[i].accumulation)
-        # println(result.config.var[i].distribution)
+        report(result)
     end
 end
 
-run(Steps)
-# @time run(Steps)
+run(Steps, :mcmc)
+run(Steps, :vegas)
+run(Steps, :vegasmc)