Merge pull request #22 from numericalEFT/dev

Dev

Author: fsxbhyy

experiment/yuanyao.txt renamed to example/QMC.txt

@@ -27,12 +27,12 @@ diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
 
 conformal_tau(τ, β) = π^(1 / 4) / sqrt(2β) * 1 / sqrt(sin(π * τ / β))
 
-function syk_sigma_dlr(d, G_x, J = 1.0)
+function syk_sigma_dlr(d, G_x, J = 1.0; sumrule = nothing, verbose = false)
 
     tau_k = d.τ # DLR imaginary time nodes
     tau_k_rev = d.β .- tau_k # Reversed imaginary time nodes
 
-    G_x_rev = tau2tau(d, G_x, tau_k_rev) # G at beta - tau_k
+    G_x_rev = tau2tau(d, G_x, tau_k_rev, sumrule = sumrule, verbose = verbose) # G at beta - tau_k
 
     Sigma_x = J .^ 2 .* G_x .^ 2 .* G_x_rev # SYK self-energy in imaginary time
 
@@ -50,15 +50,15 @@ function dyson(d, sigma_q, mu)
     end
 end
 
-function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 5000, G_x = zeros(ComplexF64, length(d)), verbose = true)
+function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 5000, G_x = zeros(ComplexF64, length(d)), sumrule = nothing, verbose = true)
 
     for iter in 1:maxiter
 
-        Sigma_x = syk_sigma_dlr(d, G_x)
+        Sigma_x = syk_sigma_dlr(d, G_x, sumrule = sumrule, verbose = verbose)
 
         G_q_new = dyson(d, tau2matfreq(d, Sigma_x), mu)
 
-        G_x_new = matfreq2tau(d, G_q_new)
+        G_x_new = matfreq2tau(d, G_q_new, sumrule = sumrule, verbose = verbose)
 
 
         if verbose
@@ -94,21 +94,24 @@ printG(dsym_correct, G_x_correct)
 
 printstyled("=====    Test Symmetrized and Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)
 
-@printf("%30s%30s%30s%15s\n", "Euv", "symmetrized solver error", "unsymmetrized solver error", "good or bad")
+@printf("%30s%30s%30s%30s%20s\n", "Euv", "sym_solver", "unsym_solver", "unsym_solver+sum_rule", "good or bad")
 for Euv in LinRange(5.0, 10.0, 50)
 
     rtol = 1e-10
     β = 10000.0
     # 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
-    @time G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix, verbose = verbose)
-    # printG(dsym, G_x_ph)
+    G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix, sumrule = nothing, verbose = verbose)
 
     # 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
-    @time G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, verbose = verbose)
+    G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, sumrule = nothing, verbose = verbose)
+
+    # printstyled("=====     Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)
+    mix = 0.01
+    G_x_none_sumrule = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix, sumrule = 1.0, verbose = verbose)
     # printG(dnone, G_x_none)
 
     # printstyled("=====     Unsymmetrized versus Symmetrized DLR solver    =======\n", color = :yellow)
@@ -122,11 +125,13 @@ for Euv in LinRange(5.0, 10.0, 50)
 
     G_x_interp_ph = tau2tau(dsym_correct, G_x_correct, dsym.τ)
     G_x_interp_none = tau2tau(dsym_correct, G_x_correct, dnone.τ)
+    G_x_interp_none_sumrule = tau2tau(dsym_correct, G_x_correct, dnone.τ)
     d_ph = diff(G_x_interp_ph, G_x_ph)
     d_none = diff(G_x_interp_none, G_x_none)
-    flag = (d_ph < 100rtol) && (d_none < 100rtol) ? "good" : "bad"
+    d_none_sumrule = diff(G_x_interp_none_sumrule, G_x_none_sumrule)
+    flag = (d_ph < 100rtol) && (d_none < 100rtol) && (d_none_sumrule < 100rtol) ? "good" : "bad"
 
-    @printf("%30.15f%30.15e%30.15e%20s\n", Euv, d_ph, d_none, flag)
+    @printf("%30.15f%30.15e%30.15e%30.15e%20s\n", Euv, d_ph, d_none, d_none_sumrule, flag)
     # println("symmetric Euv = $Euv maximumal difference: ", diff(G_x_interp, G_x_ph))
     # println("non symmetric Euv = $Euv maximumal difference: ", diff(G_x_interp, G_x_none))
 

example/SYK.jl

@@ -0,0 +1,66 @@
+using Lehmann
+using DelimitedFiles
+using Printf
+using Gaston  # for ploting
+
+β = 10.0 # inverse temperature
+Euv = 4.0 # UV energy cutoff for the spectral density
+eps = 1.0e-8 # accuracy for DLR basis
+
+# load Auxiliary field quantum Monte Carlo data for a fermionic Green's function from Yuanyao He
+data = DelimitedFiles.readdlm("QMC.txt")
+τ = data[:, 1] * β
+y = data[:, 2]
+err = data[:, 3]
+
+# construct DLR basis
+dlr = DLRGrid(Euv, β, eps, true)
+
+######################### G_dlr(τ) ################################
+printstyled("Fit without Sum Rule\n", color = :green)
+coeff = tau2dlr(dlr, y, τ; error = err)
+yp = dlr2tau(dlr, coeff, τ)
+# print comparsion for selected τ point 
+printstyled("Compare G(τ) <-->  G_dlr(τ)\n", color = :green)
+printstyled("G(0⁺)+G(β⁻)-1: ", y[1] + y[end] - 1, " vs ", yp[1] + yp[end] - 1, "\n", color = :red) #check sum rule
+printstyled("maximum dlr fitting error: ", maximum(abs.(y - yp)), "\n", color = :red)
+printstyled(" τi       G(τ)        G_dlr(τ)      δG\n", color = :yellow)
+
+printstyled("Fit with Sum Rule\n", color = :green)
+coeff_sumrule = tau2dlr(dlr, y, τ; error = err, sumrule = 1.0)
+yp = dlr2tau(dlr, coeff_sumrule, τ)
+# print comparsion for selected τ point 
+printstyled("Compare G(τ) <-->  G_dlr(τ)\n", color = :green)
+printstyled("G(0⁺)+G(β⁻)-1: ", y[1] + y[end] - 1, " vs ", yp[1] + yp[end] - 1, "\n", color = :red) #check sum rule
+printstyled("maximum dlr fitting error: ", maximum(abs.(y - yp)), "\n", color = :red)
+printstyled(" τi       G(τ)        G_dlr(τ)      δG\n", color = :yellow)
+
+for i in 1:25:length(y)
+    @printf("%3i    %10.8f    %10.8f   %10.8f\n", i, y[i], yp[i], abs(y[i] - yp[i]))
+end
+
+######################### G_dlr(iωn) ################################
+n = collect(0:50)
+Gn = dlr2matfreq(dlr, coeff, n)
+
+printstyled("Matsubara frequency G_dlr(ωn)\n", color = :green)
+printstyled(" n          Re           Im\n", color = :yellow)
+for i in 1:length(n)
+    @printf("%3i    %10.8f   %10.8f\n", n[i], real(Gn[i]), -imag(Gn[i]))
+end
+
+######################### Σ_dlr(iωn) ################################
+g0inv = @. -1im * (2 * n + 1) * π / β + 0.4
+Ginv = 1.0 ./ Gn
+Σ = g0inv - Ginv
+
+printstyled("Matsubara frequency Σ_dlr(ωn)\n", color = :green)
+printstyled(" n          Re           Im\n", color = :yellow)
+for i in 1:length(n)
+    @printf("%3i    %10.8f   %10.8f\n", n[i], real(Σ[i]), -imag(Σ[i]))
+end
+
+# set(term = "qt")
+# p = plot(n, -imag.(Σ))
+# display(p)
+# readline()

example/demoMC.jl

@@ -0,0 +1,36 @@
+using Lehmann
+using Printf
+using Gaston
+β = 40.0
+rtol = 1e-8
+eta = 1e-2
+N = 128
+
+dlr = DLRGrid(β = β, Euv = 1.0, isFermi = true, rtol = rtol, symmetry = :none)
+τgrid = collect(LinRange(0.0, β, N))
+println(τgrid)
+G = Sample.SemiCircle(dlr, :τ, τgrid)
+G += rand(length(G)) / 2.0 * eta
+coeff = tau2dlr(dlr, G, τgrid)
+coeff_sumrule = tau2dlr(dlr, G, τgrid, sumrule = π / 2)
+println("sum rule: ", abs(sum(coeff) - π / 2), " versus ", abs(sum(coeff_sumrule) - π / 2))
+
+
+dlr10 = DLRGrid(β = β, Euv = 100.0, isFermi = true, rtol = rtol, symmetry = :none)
+
+Gtrue = Sample.SemiCircle(dlr, :τ, dlr10.τ)
+Gfit = real.(dlr2tau(dlr, coeff, dlr10.τ))
+Gfit_sumrule = real.(dlr2tau(dlr, coeff_sumrule, dlr10.τ))
+
+@printf("%15s%30s%30s%30s\n", "tau", "true", "no sumrule", "with sumrule")
+for i in 1:dlr10.size
+    @printf("%15.6f%30.15f%30.15e%30.15e\n", dlr10.τ[i], Gtrue[i], abs(Gfit[i] - Gtrue[i]), abs(Gfit_sumrule[i] - Gtrue[i]))
+end
+
+set(term = "qt")
+# p = plot(dlr10.τ, Gtrue, label = "original", axis = "semilogy")
+xrange = (0.0, 40)
+p = plot(dlr10.τ, Gfit - Gtrue, plotstyle = :linespoints, leg = Symbol("DLR no sum rule"), Axes(xrange = xrange))
+plot!(dlr10.τ, Gfit_sumrule - Gtrue, plotstyle = :linespoints, leg = Symbol("DLR with sum rule"))
+display(p)
+readline()

example/demoSumRule.jl

@@ -46,35 +46,38 @@ mutable struct DLRGrid
 
     """
     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)
+    function DLRGrid(; isFermi::Bool, β = -1.0, beta = -1.0, Euv = 1.0, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
 
     Create DLR grids
 
     #Arguments:
-    - `Euv`       : the UV energy scale of the spectral density 
-    - `β`         : inverse temeprature
-    - `isFermi`   : bool is fermionic or bosonic
-    - `symmetry`  : particle-hole symmetric :ph, or particle-hole asymmetric :pha, or :none
-    - `rtol`      : tolerance absolute error
-    - `rebuild`   : set false to load DLR basis from the file, set true to recalculate the DLR basis on the fly
-    - `folder`    : if rebuild is true and folder is set, then dlrGrid will be rebuilt and saved to the specified folder
-                    if rebuild is false and folder is set, then dlrGrid will be loaded from the specified folder
-    - `algorithm` : if rebuild = true, then set :functional to use the functional algorithm to generate the DLR basis, or set :discrete to use the matrix algorithm.
-    - `verbose`   : false not to print DLRGrid to terminal, or true to print
+    - `Euv`         : the UV energy scale of the spectral density 
+    - `β` or `beta` : inverse temeprature
+    - `isFermi`     : bool is fermionic or bosonic
+    - `symmetry`    : particle-hole symmetric :ph, or particle-hole asymmetric :pha, or :none
+    - `rtol`        : tolerance absolute error
+    - `rebuild`     : set false to load DLR basis from the file, set true to recalculate the DLR basis on the fly
+    - `folder`      : if rebuild is true and folder is set, then dlrGrid will be rebuilt and saved to the specified folder
+                      if rebuild is false and folder is set, then dlrGrid will be loaded from the specified folder
+    - `algorithm`   : if rebuild = true, then set :functional to use the functional algorithm to generate the DLR basis, or set :discrete to use the matrix algorithm.
+    - `verbose`     : false not to print DLRGrid to terminal, or true to print
     """
     function DLRGrid(Euv, β, rtol, isFermi::Bool, symmetry::Symbol = :none; rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
         Λ = Euv * β # dlr only depends on this dimensionless scale
         # println("Get $Λ")
         @assert rtol > 0.0 "rtol=$rtol is not positive and nonzero!"
         @assert Λ > 0 "Energy scale $Λ must be positive!"
         @assert symmetry == :ph || symmetry == :pha || symmetry == :none "symmetry must be :ph, :pha or nothing"
+        @assert algorithm == :functional || algorithm == :discrete "Algorithm is either :functional or :discrete"
+        @assert β > 0.0 "Inverse temperature must be temperature."
+        @assert Euv > 0.0 "Energy cutoff must be positive."
 
         if Λ > 1e8 && symmetry == :none
             @warn("Current DLR without symmetry may cause ~ 3-4 digits loss for Λ ≥ 1e8!")
         end
 
-        if rtol >= 1e-6
-            @warn("Current implementation may cause ~ 3-4 digits loss for rtol ≥ 1e-6!")
+        if rtol > 1e-6
+            @warn("Current implementation may cause ~ 3-4 digits loss for rtol > 1e-6!")
         end
 
         rtolpower = Int(floor(log10(rtol))) # get the biggest n so that rtol>1e-n
@@ -144,7 +147,15 @@ mutable struct DLRGrid
         return dlr
     end
 
-    function DLRGrid(; Euv, β, isFermi::Bool, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
+    function DLRGrid(; isFermi::Bool, β = -1.0, beta = -1.0, Euv = 1.0, symmetry::Symbol = :none, rtol = 1e-14, rebuild = false, folder = nothing, algorithm = :functional, verbose = false)
+        if β <= 0.0 && beta > 0.0
+            β = beta
+        elseif β > 0.0 && beta <= 0.0
+            beta = β
+        elseif β < 0.0 && beta < 0.0
+            error("Either β or beta needs to be initialized with a positive value!")
+        end
+        @assert β ≈ beta
         return DLRGrid(Euv, β, rtol, isFermi, symmetry; rebuild = rebuild, folder = folder, algorithm = algorithm, verbose = verbose)
     end
 end