Merge pull request #103 from numericalEFT/xcai-dev

Xcai dev

Project.toml

@@ -14,11 +14,11 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-BrillouinZoneMeshes = "0.1"
+BrillouinZoneMeshes = "0.2.2"
 CompositeGrids = "0.1"
 Lehmann = "0.2"
 StaticArrays = "1"
-julia = "1.6, 1.7, 1.8"
+julia = "1.6, 1.7, 1.8, 1.9"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -75,8 +75,9 @@ The momentum is defined on the first Brillouin zone captured by a 2D k-mesh.
     using GreenFunc: BrillouinZoneMeshes
 
     DIM, nk = 2, 8
-    latvec = [1.0 0.0; 0.0 1.0] .* 2π
-    bzmesh = BrillouinZoneMeshes.BaseMesh.UniformMesh{DIM, nk}([0.0, 0.0], latvec)
+    lattice = Matrix([1.0 0; 0 1]')
+    br = BrillouinZoneMeshes.BZMeshes.Cell(lattice=lattice)
+    bzmesh = BrillouinZoneMeshes.BZMeshes.UniformBZMesh(cell=br, size=(nk, nk))
     ωₙmesh = ImFreq(10.0, FERMION)
     g_freq =  MeshArray(bzmesh, ωₙmesh; dtype=ComplexF64)
 

example/hubbard_response/HubbardRPA.jl

@@ -13,6 +13,8 @@ export hubbard_rpa, save_hubbard_rpa_list, load_hubbard_rpa_list
 
 # load python script(import from file not allowed in PythonCall)
 
+println("HubbardRPA.jl start")
+
 f = open("./example/hubbard_response/tprf_rpa.py", "r")
 # exec python script
 pyexec(read(f, String), Main)
@@ -21,10 +23,10 @@ gamma_rpa = pyeval(Py, "gamma_rpa", Main)
 @with_kw struct Para
     norb::Int = 1
     t::Float64 = 1.0
-    nk::Int = 32
+    nk::Int = 16
     dim::Int = 2
-    beta::Float64 = 10.0
-    nw::Int = 100
+    beta::Float64 = 5.0
+    nw::Int = 50
     mu::Float64 = 0.0
     U::Float64 = 1.0
 end
@@ -68,4 +70,5 @@ function load_hubbard_rpa_list(;
     end
 end
 
+println("HubbardRPA.jl end")
 end

example/hubbard_response/propagator.jl

@@ -22,25 +22,36 @@ catch
     global paras, greens, gammas = save_hubbard_rpa_list(betas=betas, fname=fname)
 end
 
-gr_w = greens[pid]
-ga_w = gammas[pid]
-para = paras[pid]
+const gr_w = greens[pid]
+const ga_w = gammas[pid]
+const para = paras[pid]
 
+const beta = para.beta
+const Euv = 100
+const extTgrid = CompositeGrid.LogDensedGrid(:uniform, [0.0, beta], [0.0, beta], 5, 1 / Euv, 5)
+# const extTgrid = CompositeGrids.SimpleG.Uniform{Float64}([0, para.beta], para.nw)
 
-gr_dlr = gr_w |> to_dlr
-gr_imt = dlr_to_imtime(gr_dlr, CompositeGrids.SimpleG.Uniform{Float64}([0, para.beta], para.nw))
+const gr_dlr = gr_w |> to_dlr
+const gr_imt = dlr_to_imtime(gr_dlr, extTgrid)
 
-ga0 = ga_w[1, 1, 1, 1, 1, 1]
+const ga0 = ga_w[1, 1, 1, 1, 1, 1]
 ga_w .-= ga0
-ga_dlr = ga_w |> to_dlr
-ga_imt = dlr_to_imtime(ga_dlr, CompositeGrids.SimpleG.Uniform{Float64}([0, para.beta], para.nw))
+const ga_dlr = ga_w |> to_dlr
+const ga_imt = dlr_to_imtime(ga_dlr, extTgrid)
 # use default here cause ERROR:
 # when dlr.τ is converted to CompositeGrid.SimpleG.Arbitrary, the bound is not [0, β]
 
+println(ga_imt.data[1, 1, 1, 1, 1, :])
+
+const rdyn = MeshArray(gr_imt.mesh[3:4]...)
+rdyn.data .= 0.0
+const r0 = MeshArray(gr_imt.mesh[3])
+r0.data .= 0.0
+
 function green(k, t1, t2)
     t = t2 - t1
     factor = 1.0
-    if t < 0
+    if t <= 0
         t = t + para.beta
         factor = -1.0
     end
@@ -56,7 +67,7 @@ end
 
 function dynamicW0(k, t1, t2)
     t = t2 - t1
-    if t < 0
+    if t <= 0
         t = t + para.beta
     end
 
@@ -65,10 +76,22 @@ function dynamicW0(k, t1, t2)
     return real(ga_imt[1, 1, 1, 1, ik, it])
 end
 
-function bareR(k)
-    return 1.0
+function bareR(k; data=r0.data)
+    # return 0.0
+    ik = locate(r0.mesh[1], k)
+    return real(data[ik])
 end
 
-function dynamicR(k, t1, t2)
-    return 0.0
+function dynamicR(k, t1, t2; data=rdyn.data)
+    # return 0.0
+    t = t2 - t1
+    factor = 1.0
+    if t <= 0
+        t = t + para.beta
+        factor = -1.0
+    end
+
+    ik, it = locate(rdyn.mesh[1], k), locate(rdyn.mesh[2], t)
+
+    return real(data[ik, it])
 end

example/hubbard_response/sc_response_mc.jl

@@ -21,12 +21,14 @@ function integrand(vars, config)
     #### instant R = instant W0 * G * G * (instant R + dynamic R) ############
     instant = green(k, t1, t3) * dynamicR(k, t3, t4) * green(k, t4, t1)
     instant += green(k, t1, t3) * bareR(k) * green(k, t3, t1) / beta
+    instant = instant + green(k, t1, t3) * green(k, t3, t1) / beta
     # both t3 and t4 will be integrated over, but bareR doesn't depend on t3, t4. So one must divide by beta
     instant *= bareW0(q)
 
     #### dynamic R = dynamic W0 * G * G * (instant R + dynamic R) ############
     dynamic = green(k, t1, t3) * dynamicR(k, t3, t4) * green(k, t4, t2)
     dynamic += green(k, t1, t3) * bareR(k) * green(k, t3, t2) / beta
+    dynamic = dynamic + green(k, t1, t3) * green(k, t3, t2) / beta
     dynamic *= dynamicW0(q, t1, t2)
 
     return instant, dynamic
@@ -47,15 +49,17 @@ function PPver(;
     para::HubbardRPA.Para=para,
     kwargs...)
     HubbardRPA.@unpack t, nk, dim, beta, U = para
+    nkf = Base.floor(Int, nk / 2 + 1)
 
     Euv = 4t
 
-    kmesh = ga_w.mesh[5]
+    kmesh = rdyn.mesh[1]
 
     ##### prepare external K grid and tau grid ##############
     extKgrid = [(k[1], k[2]) for k in kmesh]
-    extTgrid = CompositeGrid.LogDensedGrid(:uniform, [0.0, beta], [0.0, beta], 8, 1 / Euv, 8) #roughly ~100 points if resolution = β/128
-
+    # extTgrid = CompositeGrid.LogDensedGrid(:uniform, [0.0, beta], [0.0, beta], 8, 1 / Euv, 8) #roughly ~100 points if resolution = β/128
+    extTgrid = rdyn.mesh[2]
+    println(extTgrid)
     nt = length(extTgrid)
 
     latvec = 1.0 * π
@@ -72,6 +76,7 @@ function PPver(;
     # there are only 2 time variables (T[3], T[4]), while T[1] and T[2] are fixed to 0.0 and the external time
 
     obs = [zeros(Float64, nk * nk), zeros(Float64, nk * nk, nt)] # instant R and dynamic R
+    # obs = [r0.data, rdyn.data]
 
     if isnothing(config)
         config = Configuration(;
@@ -86,14 +91,19 @@ function PPver(;
     end
 
     result = integrate(integrand; config=config, measure=measure, print=print, neval=neval, kwargs...)
+    # for i in 1:9
+    #     r0.data .= result.mean[1]
+    #     rdyn.data .= result.mean[2]
+    #     result = integrate(integrand; config=config, measure=measure, print=print, neval=neval, kwargs...)
+    # end
     # # result = integrate(integrandKW; config=config, print=print, neval=neval, kwargs...)
     # # niter=niter, print=print, block=block, kwargs...)
 
     if isnothing(result) == false
 
         if print >= -1
             report(result.config)
-            println(report(result, pick=o -> first(o)))
+            println(report(result, pick=o -> o[nkf]))#irst(o)))
             println(result)
         end
         if print >= -2
@@ -110,6 +120,8 @@ function PPver(;
             # datadict[partition[o]] = data
             datadict[o] = data
         end
+        r0.data .= result.mean[1]
+        rdyn.data .= result.mean[2]
         return datadict, result
     else
         return nothing, nothing
@@ -119,4 +131,11 @@ end
 
 # solve linear response function and compute Tc 
 
-PPver(neval=1e6)
+datadict, result = PPver(neval=1e6)
+rdyn_freq = rdyn |> to_dlr |> dlr_to_imfreq
+println(r0.data[1:16])
+println(r0.data[129:144])
+println(rdyn.data[1, :])
+println(rdyn.data[9, :])
+# println(rdyn_freq.data[9, :])
+println("1/Δ0 = ", 1 / (1.0 + r0.data[9] + rdyn_freq.data[9, 1]))

src/meshgrids/MeshGrids.jl

@@ -9,11 +9,9 @@ locate(m::AbstractGrid, pos) = CompositeGrids.Interp.locate(m, pos)
 volume(m::AbstractGrid, index) = CompositeGrids.Interp.volume(m, index)
 volume(m::AbstractGrid) = CompositeGrids.Interp.volume(m)
 
-locate(m::BrillouinZoneMeshes.BaseMesh.AbstractMesh, pos) = BrillouinZoneMeshes.BaseMesh.locate(m, pos)
-volume(m::BrillouinZoneMeshes.BaseMesh.AbstractMesh, index) = BrillouinZoneMeshes.BaseMesh.volume(m, index)
-volume(m::BrillouinZoneMeshes.BaseMesh.AbstractMesh) = BrillouinZoneMeshes.BaseMesh.volume(m, pos)
-
-
+locate(m::BrillouinZoneMeshes.AbstractMesh, pos) = BrillouinZoneMeshes.AbstractMeshes.locate(m, pos)
+volume(m::BrillouinZoneMeshes.AbstractMesh, index) = BrillouinZoneMeshes.AbstractMeshes.volume(m, index)
+volume(m::BrillouinZoneMeshes.AbstractMesh) = BrillouinZoneMeshes.AbstractMeshes.volume(m)
 
 export locate, volume
 
@@ -26,7 +24,6 @@ export FERMION, BOSON
 export TemporalGrid
 include("common.jl")
 
-
 include("dlrfreq.jl")
 export DLRFreq
 

src/meshgrids/common.jl

@@ -72,4 +72,15 @@ function _grid(grid, n=3)
 end
 
 Base.show(io::IO, ::MIME"text/plain", tg::TemporalGrid) = Base.show(io, tg)
-Base.show(io::IO, ::MIME"text/html", tg::TemporalGrid) = Base.show(io, tg)
+Base.show(io::IO, ::MIME"text/html", tg::TemporalGrid) = Base.show(io, tg)
+
+_to_AbstractGrid(grid::AbstractGrid, dtype) = grid, false # do nothing if already AbstractGrid
+function _to_AbstractGrid(grid::AbstractVector, dtype)
+    # convert grid into AbstractGrid
+    rev = issorted(grid, rev=true)
+    if rev
+        grid = reverse(grid)
+    end
+    grid = SimpleG.Arbitrary{dtype}(dtype.(grid))
+    return grid, rev
+end

src/meshgrids/imfreq.jl

@@ -62,11 +62,14 @@ function ImFreq(β, isFermi::Bool=false;
         grid = SimpleG.Arbitrary{Int}(dlr.n)
         rev = false
     elseif (grid isa AbstractVector)
-        rev = issorted(grid, rev=true)
-        if rev
-            grid = reverse(grid)
-        end
-        grid = SimpleG.Arbitrary{Int}(Int.(grid))
+        # if !(grid isa AbstractGrid)
+        #     rev = issorted(grid, rev=true)
+        #     if rev
+        #         grid = reverse(grid)
+        #     end
+        #     grid = SimpleG.Arbitrary{Int}(Int.(grid))
+        # end
+        grid, rev = _to_AbstractGrid(grid, Int)
     else
         error("Proper grid or basis are required.")
     end

src/meshgrids/imtime.jl

@@ -63,12 +63,14 @@ function ImTime(β, isFermi::Bool=false;
         # grid = SimpleG.Uniform{dtype}([0, β], Int(round(β / resolution)))
         # grid = CompositeGrid.LogDensedGrid(:uniform, [0.0, β], [0.0, β], 8, 1 / Euv, 8) #roughly ~100 points if resolution = β/128
     elseif (grid isa AbstractVector)
-        rev = issorted(grid, rev=true)
-        if rev
-            grid = reverse(grid)
-        end
-
-        grid = SimpleG.Arbitrary{dtype}(grid)
+        # if !(grid isa AbstractGrid)
+        #     rev = issorted(grid, rev=true)
+        #     if rev
+        #         grid = reverse(grid)
+        #     end
+        #     grid = SimpleG.Arbitrary{dtype}(grid)
+        # end
+        grid, rev = _to_AbstractGrid(grid, dtype)
     else
         error("Proper grid and basis are required!")
     end

test/test_MeshGrids.jl

@@ -177,4 +177,46 @@
         end
         @test volume(tg2) ≈ sum(volume(tg2, i) for i in 1:length(tg2))
     end
+
+    @testset "CompositeGrids and BrillouinZoneMeshes" begin
+        # test availability of imported package
+
+        # CompositeGrids
+        N1 = 7
+        mesh1 = MeshGrids.SimpleGrid.Uniform{Float64}([0.0, 1.0], N1)
+        for (i, x) in enumerate(mesh1)
+            @test MeshGrids.locate(mesh1, x) == i
+        end
+        @test MeshGrids.volume(mesh1) ≈ sum(MeshGrids.volume(mesh1, i) for i in 1:length(mesh1))
+
+        # test if CompositeGrid could be handled correctly
+        β = 1.0
+        # ImTime
+        mesh1 = [0.0, 0.5β, β]
+        mesh2 = MeshGrids.CompositeGrid.LogDensedGrid(:cheb, [0.0, β], [0.0, β], 4, 0.01, 4)
+        itm1 = MeshGrids.ImTime(β; grid=mesh1)
+        itm2 = MeshGrids.ImTime(β; grid=mesh2)
+        @test itm1.grid isa AbstractGrid
+        @test itm2.grid isa AbstractGrid
+        @test itm2.grid == mesh2 # type of mesh2 should be preserved
+
+        # ImFreq
+        mesh1 = [1, 2, 3]
+        mesh2 = MeshGrids.SimpleGrid.Arbitrary{Int}([1, 2, 3, 4])
+        ifm1 = MeshGrids.ImFreq(β; grid=mesh1)
+        ifm2 = MeshGrids.ImFreq(β; grid=mesh2)
+        @test ifm1.grid isa AbstractGrid
+        @test ifm2.grid isa AbstractGrid
+        @test ifm2.grid == mesh2 # type of mesh2 should be preserved
+
+        # BrillouinZoneMeshes
+        lattice = Matrix([2.0 0 0; 1 sqrt(3) 0; 7 11 19]')
+        msize = (3, 5, 7)
+        br = MeshGrids.BZMeshes.Cell(lattice=lattice)
+        mesh2 = MeshGrids.BZMeshes.UniformBZMesh(cell=br, size=msize)
+        for (i, x) in enumerate(mesh2)
+            @test MeshGrids.locate(mesh2, x) == i
+        end
+        @test MeshGrids.volume(mesh2) ≈ sum(MeshGrids.volume(mesh2, i) for i in 1:length(mesh2))
+    end
 end