Merge pull request #20 from numericalEFT/dev

Dev

Author: fsxbhyy

Project.toml

@@ -1,7 +1,7 @@
 name = "Lehmann"
 uuid = "95bf888a-8996-4655-9f35-1c0506bdfefe"
 authors = ["Kun Chen", "Tao Wang", "Xiansheng Cai"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

docs/src/manual/kernel.md

@@ -23,11 +23,11 @@ We use the following conventions:
 - Fourier transform follows the convention in the book "Quantum Many-particle Systems" by J. Negele and H. Orland, Page 95,
 
 ```math
-G(\tau) = \frac{1}{\beta} \sum_n G(i\omega_n) \text{e}^{i\omega_n \tau}
+G(\tau) = \frac{1}{\beta} \sum_n G(i\omega_n) \text{e}^{-i\omega_n \tau}
 ```
 
 ```math
-G(i\omega_n) = \int_0^\beta G(\tau) \text{e}^{-i\omega_n \tau} d\tau
+G(i\omega_n) = \int_0^\beta G(\tau) \text{e}^{i\omega_n \tau} d\tau
 ```
 
 # Fermion without Symmetry 
@@ -65,7 +65,7 @@ K(τ, ω) = e^{-ω|τ|}+e^{-ω(β-|τ|)}
 ```
 - Matusbara frequency
 ```math
-K(iω_n, ω) = -\frac{2iω_n}{ω^2+ω_n^2}(1+e^{-ωβ}),
+K(iω_n, ω) = \frac{2iω_n}{ω^2+ω_n^2}(1+e^{-ωβ}),
 ```
 
 # Boson with the Particle-hole Symmetry 
@@ -104,5 +104,5 @@ K(τ, ω) = e^{-ω|τ|}-e^{-ω(β-|τ|)}
 ```
 - Matusbara frequency
 ```math
-K(iω_n, ω) = -\frac{2iω_n}{ω^2+ω_n^2}(1-e^{-ωβ}),
+K(iω_n, ω) = \frac{2iω_n}{ω^2+ω_n^2}(1-e^{-ωβ}),
 ```

example/SYK.jl

@@ -42,7 +42,7 @@ end
 function dyson(d, sigma_q, mu)
     if d.symmetry == :ph #symmetrized G
         @assert mu ≈ 0.0 "Only the case μ=0 enjoys the particle-hole symmetry."
-        return 1im * imag.(1 ./ (d.ωn * 1im .+ mu .- sigma_q))
+        return 1im * imag.(-1 ./ (d.ωn * 1im .- mu .+ sigma_q))
     elseif d.symmetry == :none
         return -1 ./ (d.ωn * 1im .- mu .+ sigma_q)
     else
@@ -102,13 +102,13 @@ for Euv in LinRange(5.0, 10.0, 50)
     # printstyled("=====     Symmetrized DLR solver for SYK model     =======\n", color = :yellow)
     mix = 0.01
     dsym = DLRGrid(Euv = Euv, β = β, isFermi = true, rtol = rtol, symmetry = :ph, rebuild = true, verbose = false) # Initialize DLR object
-    G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix, verbose = verbose)
+    @time G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix, verbose = verbose)
     # printG(dsym, G_x_ph)
 
     # printstyled("=====     Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)
     mix = 0.01
     dnone = DLRGrid(Euv = Euv, β = β, isFermi = true, rtol = rtol, symmetry = :none, rebuild = true, verbose = false) # Initialize DLR object
-    G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, verbose = verbose)
+    @time G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, verbose = verbose)
     # printG(dnone, G_x_none)
 
     # printstyled("=====     Unsymmetrized versus Symmetrized DLR solver    =======\n", color = :yellow)

src/discrete/kernel.jl

@@ -103,11 +103,11 @@ function preciseKernelT(dlrGrid, τ, ω, print::Bool = true)
         #symmetrize K(τ, ω)=K(β-τ, -ω) for τ>0 
         @assert isodd(τ.np) #symmetrization is only possible for odd τ panels
         halfτ = ((τ.np - 1) ÷ 2) * τ.degree
-        kernel[1:halfτ, :] = Spectral.kernelT(true, symmetry, τ.grid[1:halfτ], ω.grid, 1.0, true)
-        kernel[end:-1:(halfτ+1), :] = Spectral.kernelT(true, symmetry, τ.grid[1:halfτ], ω.grid[end:-1:1], 1.0, true)
+        kernel[1:halfτ, :] = Spectral.kernelT(Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid, 1.0, true)
+        kernel[end:-1:(halfτ+1), :] = Spectral.kernelT(Val(true), Val(symmetry), τ.grid[1:halfτ], ω.grid[end:-1:1], 1.0, true)
         # use the fermionic kernel for both the fermionic and bosonic propagators
     else
-        kernel = Spectral.kernelT(dlrGrid.isFermi, symmetry, τGrid, ωGrid, 1.0, true)
+        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(symmetry), τGrid, ωGrid, 1.0, true)
     end
 
     # print && println("=====  Kernel Discretization =====")
@@ -191,8 +191,8 @@ function preciseKernelΩn(dlrGrid, ω, print::Bool = true)
     symmetry = dlrGrid.symmetry
     n = nGrid(dlrGrid.isFermi, symmetry, dlrGrid.Λ)
 
-    nkernelFermi = Spectral.kernelΩ(true, symmetry, n, Float64.(ωGrid), 1.0, true)
-    nkernelBose = Spectral.kernelΩ(false, symmetry, n, Float64.(ωGrid), 1.0, true)
+    nkernelFermi = Spectral.kernelΩ(Val(true), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
+    nkernelBose = Spectral.kernelΩ(Val(false), Val(symmetry), n, Float64.(ωGrid), 1.0, true)
 
     return n, nkernelFermi, nkernelBose
 end

src/dlr.jl

@@ -26,7 +26,7 @@ struct DLRGrid
 - `ωn` : (2n+1)π/β
 - `τ` : selected representative imaginary-time grid
 """
-struct DLRGrid
+mutable struct DLRGrid
     isFermi::Bool
     symmetry::Symbol
     Euv::Float64
@@ -41,6 +41,9 @@ struct DLRGrid
     ωn::Vector{Float64} # (2n+1)π/β
     τ::Vector{Float64}
 
+    kernel_τ::Any
+    kernel_n::Any
+
     """
     function DLRGrid(Euv, β, rtol, isFermi::Bool; symmetry::Symbol = :none, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
     function DLRGrid(; Euv, β, isFermi::Bool, symmetry::Symbol = :none, rtol = 1e-8, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
@@ -110,41 +113,34 @@ struct DLRGrid
 
         if rebuild == false
             if isnothing(folder)
-                Λfloor = Λ < 100 ? Int(100) : 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n
-
+                Λ = Λ < 100 ? Int(100) : 10^(Int(ceil(log10(Λ)))) # get smallest n so that Λ<10^n
                 folderList = [string(@__DIR__, "/../basis/"),]
-                file = filename(Λfloor, rtolpower)
-
-                dlrfile = finddlr(folderList, file)
-
-                if isnothing(dlrfile) == false
-                    dlr = new(isFermi, symmetry, Euv, β, Λfloor, 10.0^(float(rtolpower)), [], [], [], [])
-                    _load!(dlr, dlrfile, verbose)
-                    return dlr
-                else
-                    @warn("No DLR is found in the folder $folder, try to rebuild instead.")
-                end
             else
-                file = filename(Euv * β, rtolpower)
                 folderList = [folder,]
+            end
 
-                dlrfile = finddlr(folderList, file)
+            file = filename(Λ, rtolpower)
+            dlrfile = finddlr(folderList, file)
 
-                if isnothing(dlrfile) == false
-                    dlr = new(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], [])
-                    _load!(dlr, dlrfile, verbose)
-                    return dlr
-                else
-                    @warn("No DLR is found in the folder $folder, try to rebuild instead.")
-                end
+            if isnothing(dlrfile) == false
+                dlr = new(isFermi, symmetry, Euv, β, Λ, rtol, [], [], [], [], nothing, nothing)
+                _load!(dlr, dlrfile, verbose)
+                dlr.kernel_τ = Spectral.kernelT(Val(dlr.isFermi), Val(dlr.symmetry), dlr.τ, dlr.ω, dlr.β, true)
+                dlr.kernel_n = Spectral.kernelΩ(Val(dlr.isFermi), Val(dlr.symmetry), dlr.n, dlr.ω, dlr.β, true)
+                return dlr
+            else
+                @warn("No DLR is found in the folder $folder, try to rebuild instead.")
             end
 
         end
 
         # try to rebuild the dlrGrid
-        dlr = new(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], [])
+        dlr = new(isFermi, symmetry, Euv, β, Euv * β, rtol, [], [], [], [], nothing, nothing)
         file2save = filename(Euv * β, rtolpower)
         _build!(dlr, folder, file2save, algorithm, verbose)
+
+        dlr.kernel_τ = Spectral.kernelT(Val(dlr.isFermi), Val(dlr.symmetry), dlr.τ, dlr.ω, dlr.β, true)
+        dlr.kernel_n = Spectral.kernelΩ(Val(dlr.isFermi), Val(dlr.symmetry), dlr.n, dlr.ω, dlr.β, true)
         return dlr
     end
 

src/operation.jl

@@ -89,10 +89,20 @@ function tau2dlr(dlrGrid::DLRGrid, green, τGrid = dlrGrid.τ; error = nothing,
     @assert size(green)[axis] == length(τGrid)
     ωGrid = dlrGrid.ω
 
-    kernel = Spectral.kernelT(dlrGrid.isFermi, dlrGrid.symmetry, τGrid, ωGrid, dlrGrid.β, true)
-    typ = promote_type(eltype(kernel), eltype(green))
-    kernel = convert.(typ, kernel)
-    green = convert.(typ, green)
+    typ = promote_type(eltype(dlrGrid.kernel_τ), eltype(green))
+
+    if length(τGrid) == dlrGrid.size && isapprox(τGrid, dlrGrid.τ; rtol = 1e-14)
+        kernel = dlrGrid.kernel_τ
+    else
+        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), τGrid, ωGrid, dlrGrid.β, true; type = typ)
+    end
+
+    if typ != eltype(kernel)
+        kernel = convert.(typ, kernel)
+    end
+    if typ != eltype(green)
+        green = convert.(typ, green)
+    end
 
     g, partialsize = _tensor2matrix(green, axis)
 
@@ -102,7 +112,7 @@ function tau2dlr(dlrGrid::DLRGrid, green, τGrid = dlrGrid.τ; error = nothing,
     end
     coeff = _weightedLeastSqureFit(g, error, kernel)
 
-    if all(abs.(coeff) .< 1e16) == false
+    if all(x -> abs(x) < 1e16, coeff) == false
         @warn("Some of the DLR coefficients are larger than 1e16. The quality of DLR fitting could be bad.")
     end
 
@@ -127,7 +137,11 @@ function dlr2tau(dlrGrid::DLRGrid, dlrcoeff, τGrid = dlrGrid.τ; axis = 1)
     β = dlrGrid.β
     ωGrid = dlrGrid.ω
 
-    kernel = Spectral.kernelT(dlrGrid.isFermi, dlrGrid.symmetry, τGrid, ωGrid, β, true)
+    if length(τGrid) == dlrGrid.size && isapprox(τGrid, dlrGrid.τ; rtol = 1e-14)
+        kernel = dlrGrid.kernel_τ
+    else
+        kernel = Spectral.kernelT(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), τGrid, ωGrid, dlrGrid.β, true)
+    end
 
     coeff, partialsize = _tensor2matrix(dlrcoeff, axis)
 
@@ -154,10 +168,21 @@ function matfreq2dlr(dlrGrid::DLRGrid, green, nGrid = dlrGrid.n; error = nothing
     @assert eltype(nGrid) <: Integer
     ωGrid = dlrGrid.ω
 
-    kernel = Spectral.kernelΩ(dlrGrid.isFermi, dlrGrid.symmetry, nGrid, ωGrid, dlrGrid.β, true)
-    typ = promote_type(eltype(kernel), eltype(green))
-    kernel = convert.(typ, kernel)
-    green = convert.(typ, green)
+    typ = promote_type(eltype(dlrGrid.kernel_n), eltype(green))
+
+    if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol = 1e-14)
+        kernel = dlrGrid.kernel_n
+    else
+        kernel = Spectral.kernelΩ(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), nGrid, ωGrid, dlrGrid.β, true; type = typ)
+    end
+
+    if typ != eltype(green)
+        green = convert.(typ, green)
+    end
+
+    if typ != eltype(kernel)
+        dlrGrid.kernel_n = convert.(typ, dlrGrid.kernel_n)
+    end
 
     g, partialsize = _tensor2matrix(green, axis)
 
@@ -166,7 +191,7 @@ function matfreq2dlr(dlrGrid::DLRGrid, green, nGrid = dlrGrid.n; error = nothing
         error, psize = _tensor2matrix(error, axis)
     end
     coeff = _weightedLeastSqureFit(g, error, kernel)
-    if all(abs.(coeff) .< 1e16) == false
+    if all(x -> abs(x) < 1e16, coeff) == false
         @warn("Some of the DLR coefficients are larger than 1e16. The quality of DLR fitting could be bad.")
     end
     return _matrix2tensor(coeff, partialsize, axis)
@@ -189,7 +214,11 @@ function dlr2matfreq(dlrGrid::DLRGrid, dlrcoeff, nGrid = dlrGrid.n; axis = 1)
     @assert eltype(nGrid) <: Integer
     ωGrid = dlrGrid.ω
 
-    kernel = Spectral.kernelΩ(dlrGrid.isFermi, dlrGrid.symmetry, nGrid, ωGrid, dlrGrid.β, true)
+    if length(nGrid) == dlrGrid.size && isapprox(nGrid, dlrGrid.n; rtol = 1e-14)
+        kernel = dlrGrid.kernel_n
+    else
+        kernel = Spectral.kernelΩ(Val(dlrGrid.isFermi), Val(dlrGrid.symmetry), nGrid, ωGrid, dlrGrid.β, true)
+    end
 
     coeff, partialsize = _tensor2matrix(dlrcoeff, axis)
 

src/sample/sample.jl

@@ -89,15 +89,19 @@ function _Green(::Val{IsMatFreq}, Euv, β, isFermi, Grid, symmetry, n, pbp, npo,
                     #spectral density is defined for positivie frequency only for correlation functions
                     continue
                 end
-                ker = IsMatFreq ? Spectral.kernelΩ(isFermi, symmetry, τ, Euv * x, β, regularized) : Spectral.kernelT(isFermi, symmetry, τ, Euv * x, β, regularized)
+                ker = IsMatFreq ?
+                      Spectral.kernelΩ(Val(isFermi), Val(symmetry), τ, Euv * x, β, regularized) :
+                      Spectral.kernelT(Val(isFermi), Val(symmetry), τ, Euv * x, β, regularized)
                 G[τi] += (b - a) / 2 * wl[jj] * ker * sqrt(1 - x^2)
             end
         end
 
         a, b = 1.0 / 2, 1.0
         for jj = 1:n
             x = (a + b) / 2 + (b - a) / 2 * xj[jj]
-            ker = IsMatFreq ? Spectral.kernelΩ(isFermi, symmetry, τ, Euv * x, β, regularized) : Spectral.kernelT(isFermi, symmetry, τ, Euv * x, β, regularized)
+            ker = IsMatFreq ?
+                  Spectral.kernelΩ(Val(isFermi), Val(symmetry), τ, Euv * x, β, regularized) :
+                  Spectral.kernelT(Val(isFermi), Val(symmetry), τ, Euv * x, β, regularized)
             G[τi] += ((b - a) / 2)^1.5 * wj[jj] * ker * sqrt(1 + x)
         end
 
@@ -106,7 +110,9 @@ function _Green(::Val{IsMatFreq}, Euv, β, isFermi, Grid, symmetry, n, pbp, npo,
             a, b = -1.0, -1.0 / 2
             for jj = 1:n
                 x = (a + b) / 2 + (b - a) / 2 * (-xj[n-jj+1])
-                ker = IsMatFreq ? Spectral.kernelΩ(isFermi, symmetry, τ, Euv * x, β, regularized) : Spectral.kernelT(isFermi, symmetry, τ, Euv * x, β, regularized)
+                ker = IsMatFreq ?
+                      Spectral.kernelΩ(Val(isFermi), Val(symmetry), τ, Euv * x, β, regularized) :
+                      Spectral.kernelT(Val(isFermi), Val(symmetry), τ, Euv * x, β, regularized)
                 G[τi] += ((b - a) / 2)^1.5 * wj[n-jj+1] * ker * sqrt(1 - x)
             end
         end
@@ -153,9 +159,9 @@ function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbo
             end
 
             if IsMatFreq == false
-                g[τi] += Spectral.kernelT(isFermi, symmetry, τ, ω, β, regularized)
+                g[τi] += Spectral.kernelT(Val(isFermi), Val(symmetry), τ, ω, β, regularized)
             else
-                g[τi] += Spectral.kernelΩ(isFermi, symmetry, τ, ω, β, regularized)
+                g[τi] += Spectral.kernelΩ(Val(isFermi), Val(symmetry), τ, ω, β, regularized)
             end
         end
     end