Merge pull request #16 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.0"
+version = "0.2.1"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

README.md

@@ -33,6 +33,7 @@ This package has been registered. So, simply type `import Pkg; Pkg.add("Lehmann"
 In the following [demo](example/demo.jl), we will show how to compress a Green's function of ~10000 data point into ~20 DLR coefficients, and perform fast interpolation and fourier transform up to the accuracy ~1e-10.
 
 ```julia
+using Lehmann
 β = 100.0 # inverse temperature
 Euv = 1.0 # ultraviolt energy cutoff of the Green's function
 rtol = 1e-8 # accuracy of the representation
@@ -41,9 +42,6 @@ symmetry = :none # :ph if particle-hole symmetric, :pha is antisymmetric, :none
 
 diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 
-# Use semicircle spectral density to generate the sample Green's function
-sample(grid, type) = Sample.SemiCircle(Euv, β, isFermi, grid, type, symmetry)
-
 dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters and basis
 # A set of most representative grid points are generated:
 # dlr.ω gives the real-frequency grids
@@ -53,9 +51,9 @@ dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters a
 println("Prepare the Green's function sample ...")
 Nτ, Nωn = 10000, 10000 # many τ and n points are needed because Gτ is quite singular near the boundary
 τgrid = collect(LinRange(0.0, β, Nτ))  # create a τ grid
-Gτ = sample(τgrid, :τ)
+Gτ = Sample.SemiCircle(dlr, :τ, τgrid) # Use semicircle spectral density to generate the sample Green's function in τ
 ngrid = collect(-Nωn:Nωn)  # create a set of Matsubara-frequency points
-Gn = sample(ngrid, :ωn)
+Gn = Sample.SemiCircle(dlr, :n, ngrid) # Use semicircle spectral density to generate the sample Green's function in ωn
 
 println("Compress Green's function into ~20 coefficients ...")
 spectral_from_Gτ = tau2dlr(dlr, Gτ, τgrid)
@@ -64,9 +62,9 @@ spectral_from_Gω = matfreq2dlr(dlr, Gn, ngrid)
 
 println("Prepare the target Green's functions to benchmark with ...")
 τ = collect(LinRange(0.0, β, Nτ * 2))  # create a dense τ grid to interpolate
-Gτ_target = sample(τ, :τ)
+Gτ_target = Sample.SemiCircle(dlr, :τ, τ)
 n = collect(-2Nωn:2Nωn)  # create a set of Matsubara-frequency points
-Gn_target = sample(n, :ωn)
+Gn_target = Sample.SemiCircle(dlr, :n, n)
 
 println("Interpolation benchmark ...")
 Gτ_interp = dlr2tau(dlr, spectral_from_Gτ, τ)

example/Untitled.ipynb

@@ -7,9 +7,6 @@ symmetry = :none # :ph if particle-hole symmetric, :pha is antisymmetric, :none
 
 diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 
-# Use semicircle spectral density to generate the sample Green's function
-sample(grid, type) = Sample.SemiCircle(Euv, β, isFermi, grid, type, symmetry)
-
 dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters and basis
 # A set of most representative grid points are generated:
 # dlr.ω gives the real-frequency grids
@@ -19,9 +16,9 @@ dlr = DLRGrid(Euv, β, rtol, isFermi, symmetry) #initialize the DLR parameters a
 println("Prepare the Green's function sample ...")
 Nτ, Nωn = 10000, 10000 # many τ and n points are needed because Gτ is quite singular near the boundary
 τgrid = collect(LinRange(0.0, β, Nτ))  # create a τ grid
-Gτ = sample(τgrid, :τ)
+Gτ = Sample.SemiCircle(dlr, :τ, τgrid) # Use semicircle spectral density to generate the sample Green's function in τ
 ngrid = collect(-Nωn:Nωn)  # create a set of Matsubara-frequency points
-Gn = sample(ngrid, :ωn)
+Gn = Sample.SemiCircle(dlr, :n, ngrid) # Use semicircle spectral density to generate the sample Green's function in ωn
 
 println("Compress Green's function into ~20 coefficients ...")
 spectral_from_Gτ = tau2dlr(dlr, Gτ, τgrid)
@@ -30,9 +27,9 @@ spectral_from_Gω = matfreq2dlr(dlr, Gn, ngrid)
 
 println("Prepare the target Green's functions to benchmark with ...")
 τ = collect(LinRange(0.0, β, Nτ * 2))  # create a dense τ grid to interpolate
-Gτ_target = sample(τ, :τ)
+Gτ_target = Sample.SemiCircle(dlr, :τ, τ)
 n = collect(-2Nωn:2Nωn)  # create a set of Matsubara-frequency points
-Gn_target = sample(n, :ωn)
+Gn_target = Sample.SemiCircle(dlr, :n, n)
 
 println("Interpolation benchmark ...")
 Gτ_interp = dlr2tau(dlr, spectral_from_Gτ, τ)

example/demo.jl

@@ -6,14 +6,14 @@ d = DLRGrid(Euv = 1.0, β = 1000.0, rtol = 1e-14, isFermi = false) # Initialize
 tau_k = d.τ  # DLR imaginary time points
 
 println("Prepare the Green's function sample ...")
-G_k = Sample.SemiCircle(d.Euv, d.β, d.isFermi, tau_k, :τ)  # Evaluate known G at tau_k
+G_k = Sample.SemiCircle(d, :τ, tau_k)  # Evaluate known G at tau_k
 G_x = tau2dlr(d, G_k) # DLR coeffs from G_k
 
 println("Interpolate imaginary-time Green's function ...")
 tau_i = collect(LinRange(0, d.β, 40)) # Equidistant tau grid
 G_i = dlr2tau(d, G_x, tau_i) # Evaluate DLR at tau_i
 
-G_s = Sample.SemiCircle(d.Euv, d.β, d.isFermi, tau_i, :τ)  # Evaluate known G at tau_k
+G_s = Sample.SemiCircle(d, :τ, tau_i)  # Evaluate known G at tau_k
 
 diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 println("Maximum difference: ", diff(G_i, G_s))

example/demo1.jl

@@ -10,7 +10,7 @@ d = DLRGrid(Euv = 1.0, β = 1000.0, isFermi = false) # Initialize DLR object
 tau_i = collect(LinRange(0, d.β, 100))
 
 println("Prepare the Green's function sample ...")
-G_i = Sample.SemiCircle(d.Euv, d.β, d.isFermi, tau_i, :τ)  # Evaluate known G at tau_k
+G_i = Sample.SemiCircle(d, :τ, tau_i)  # Evaluate known G at tau_k
 # G_i = true_G_tau(tau_i, beta)
 
 G_i_noisy = G_i .+ eta .* (rand(length(G_i)) .- 0.5)

example/demo2.jl

@@ -47,7 +47,7 @@ function syk_sigma_dlr(d, G_x, J = 1.0)
     return Sigma_x
 end
 
-function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol; mix = 0.3, maxiter = 5, G_x = zeros(ComplexF64, length(d)))
+function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol; mix = 0.1, maxiter = 1000, G_x = zeros(ComplexF64, length(d)))
 
     for iter in 1:maxiter
 
@@ -62,7 +62,7 @@ function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol; mix = 0.3, maxiter =
         G_q_new = -1 ./ (d.ωn * 1im .- mu .+ tau2matfreq(d, Sigma_x))
         +1 ./ (-d.ωn * 1im .- mu .+ tau2matfreq(d, Sigma_x)) # Solve Dyson
 
-        G_x_new = (matfreq2tau(d, G_q_new))
+        G_x_new = matfreq2tau(d, G_q_new)
 
         println("imag", maximum(abs.(imag.(G_x_new))))
 
@@ -77,7 +77,7 @@ function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol; mix = 0.3, maxiter =
 end
 
 d = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-10, symmetry = :ph) # Initialize DLR object
-G_q = solve_syk_with_fixpoint_iter(d, 0.0)
+G_q = solve_syk_with_fixpoint_iter(d, 0.00)
 # G_q = solve_syk_with_fixpoint_iter(d, 0.15, G_x = G_q)
 # G_q = solve_syk_with_fixpoint_iter(d, 0.12, G_x = G_q)
 # G_q = solve_syk_with_fixpoint_iter(d, 0.11, G_x = G_q)

example/demo3.jl

@@ -52,14 +52,16 @@ function _weightedLeastSqureFit(Gτ, error, kernel)
         C = Gτ
     else
         @assert size(error) == size(Gτ)
-        for i = 1:size(error)[1]
-            error[i, :] /= sum(error[i, :]) / length(error[i, :])
+        w = 1.0 ./ (error .+ 1e-16)
+
+        for i = 1:size(w)[1]
+            w[i, :] /= sum(w[i, :]) / length(w[i, :])
         end
 
         # W = Diagonal(weight)
         # B = transpose(kernel) * W * kernel
         # C = transpose(kernel) * W * Gτ
-        w = 1.0 ./ error
+        # w = 1.0 ./ (error .+ 1e-16)
         B = w .* kernel
         # B = Diagonal(w) * kernel
         C = w .* Gτ
@@ -99,6 +101,11 @@ function tau2dlr(dlrGrid::DLRGrid, green, τGrid = dlrGrid.τ; error = nothing,
         error, psize = _tensor2matrix(error, axis)
     end
     coeff = _weightedLeastSqureFit(g, error, kernel)
+
+    if all(abs.(coeff) .< 1e16) == 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)
 end
 
@@ -159,6 +166,9 @@ 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
+        @warn("Some of the DLR coefficients are larger than 1e16. The quality of DLR fitting could be bad.")
+    end
     return _matrix2tensor(coeff, partialsize, axis)
 end
 

src/operation.jl

@@ -3,16 +3,18 @@ using FastGaussQuadrature
 using ..Spectral
 """
     SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol = :none; rtol = nothing, degree = 24, regularized::Bool = true)
+    SemiCircle(dlr, type::Symbol, Grid = dlrGrid(dlr, type); degree = 24, regularized::Bool = true)
 
 Generate Green's function from a semicircle spectral density. 
 Return the function on Grid and the systematic error.
 
 #Arguments
+- `dlr`: dlrGrid struct
 - `Euv` : ultraviolet energy cutoff
 - `β` : inverse temperature
 - `isFermi`: is fermionic or bosonic
 - `Grid`: grid to evalute on
-- `type`: imaginary-time with :τ, or Matsubara-frequency with :ωn
+- `type`: imaginary-time with :τ, or Matsubara-frequency with :n
 - `symmetry`: particle-hole symmetric :ph, particle-hole antisymmetric :pha, or :none
 - `rtol`: accuracy to achieve
 - `degree`: polynomial degree for integral
@@ -23,10 +25,10 @@ function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol
     # S(ω) = sqrt(1 - (ω / Euv)^2) / Euv # semicircle -1<ω<1
     if type == :τ
         IsMatFreq = false
-    elseif type == :ωn
+    elseif type == :n
         IsMatFreq = true
     else
-        error("not implemented!")
+        error("$type is not implemented!")
     end
 
     ##### Panels endpoints for composite quadrature rule ###
@@ -48,6 +50,20 @@ function SemiCircle(Euv, β, isFermi::Bool, Grid, type::Symbol, symmetry::Symbol
     end
     return g1
 end
+function SemiCircle(dlr, type::Symbol, Grid = dlrGrid(dlr, type); degree = 24, regularized::Bool = true)
+    return SemiCircle(dlr.Euv, dlr.β, dlr.isFermi, Grid, type, dlr.symmetry; rtol = dlr.rtol, degree = degree, regularized = regularized)
+end
+
+function dlrGrid(dlr, type::Symbol)
+    if type == :τ
+        return dlr.τ
+    elseif type == :n
+        return dlr.n
+    else
+        error("$type not implemented!")
+    end
+end
+
 
 # function getG(::Val{true}, Grid)
 #     return zeros(ComplexF64, length(Grid))
@@ -100,16 +116,18 @@ end
 
 """
     MultiPole(β, isFermi::Bool, symmetry::Symbol, Grid, type::Symbol, poles, regularized::Bool = true)
+    MultiPole(dlr, type::Symbol, poles, Grid = dlrGrid(dlr, type); regularized::Bool = true)
 
 Generate Green's function from a spectral density with delta peaks specified by the argument ``poles``. 
 Return the function on Grid and the systematic error.
 
 #Arguments
+- `dlr`: dlrGrid struct
 - `β` : inverse temperature
 - `isFermi`: is fermionic or bosonic
 - `symmetry`: particle-hole symmetric :ph, particle-hole antisymmetric :pha, or :none
 - `Grid`: grid to evalute on
-- `type`: imaginary-time with :τ, or Matsubara-frequency with :ωn
+- `type`: imaginary-time with :τ, or Matsubara-frequency with :n
 - `poles`: a list of frequencies for the delta functions
 - `regularized`: use regularized bosonic kernel if symmetry = :none
 """
@@ -119,10 +137,10 @@ function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbo
     # poles = [0.0]
     if type == :τ
         IsMatFreq = false
-    elseif type == :ωn
+    elseif type == :n
         IsMatFreq = true
     else
-        error("not implemented!")
+        error("$type is not implemented!")
     end
 
     g = IsMatFreq ? zeros(ComplexF64, length(Grid)) : zeros(Float64, length(Grid))
@@ -143,5 +161,8 @@ function MultiPole(β, isFermi::Bool, Grid, type::Symbol, poles, symmetry::Symbo
     end
     return g
 end
+function MultiPole(dlr, type::Symbol, poles, Grid = dlrGrid(dlr, type); regularized::Bool = true)
+    return MultiPole(dlr.β, dlr.isFermi, Grid, type, poles, dlr.symmetry; regularized = regularized)
+end
 
 end

src/sample/sample.jl

@@ -31,7 +31,7 @@ Compute the imaginary-time kernel of different type.
         return isFermi ? kernelFermiT_PHA(τ, ω, β) : kernelBoseT_PHA(τ, ω, β)
     else
         @error "Symmetry $symmetry is not implemented!"
-        return 0.0
+        return T(0)
     end
 end
 """
@@ -259,7 +259,7 @@ Compute the imaginary-time kernel of different type. Assume ``k_B T/\\hbar=1``
         return isFermi ? kernelFermiΩ_PHA(n, ω, β) : kernelBoseΩ_PHA(n, ω, β)
     else
         @error "Symmetry $symmetry  is not implemented!"
-        return 0.0
+        return T(0)
     end
 end