Merge pull request #19 from numericalEFT/dev

Dev

.gitignore

@@ -1,4 +1,5 @@
 *.jl.*.cov
 *.jl.cov
 *.jl.mem
+Manifest.toml
 /docs/build/

Project.toml

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

build.jl

@@ -3,13 +3,14 @@ using Lehmann
 using Printf
 
 # rtol = [-6, -8, -10, -12, -14]
-rtol = [-7, -9, -11, -13]
-Λ = [100, 1000, 10000, 100000, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14]
+# rtol = [-7, -9, -11, -13]
+# Λ = [100, 1000, 10000, 100000, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14]
 # Λ = [1e11, 1e12, 1e13, 1e14]
-# rtol = [-12]
-# Λ = [1000000,]
+rtol = [-12]
+Λ = [1000000,]
 algorithm = :functional
-folder = "./basis/"
+# folder = "./basis/"
+folder = "./"
 
 for lambda in Λ
     for err in rtol

example/SYK.jl

@@ -17,7 +17,7 @@ A SYK model solver based on a forward fixed-point iteration method.
 
  The SYK Green's function has particle-hole symmetry when μ=0. 
  You may enforce such symmetry by setting `symmetry = :ph` when initialize the DLR grids.
- A symmetrized solver tends to be more robust than the unsymmetrized counterpart.
+ A symmetrized solver tends to be more robust than a unsymmetrized one.
 """
 
 using Lehmann
@@ -50,7 +50,7 @@ function dyson(d, sigma_q, mu)
     end
 end
 
-function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 1000, G_x = zeros(ComplexF64, length(d)))
+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)
 
     for iter in 1:maxiter
 
@@ -61,8 +61,10 @@ function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxit
         G_x_new = matfreq2tau(d, G_q_new)
 
 
-        if iter % (maxiter / 10) == 0
-            println("round $iter: change $(diff(G_x_new, G_x))")
+        if verbose
+            if iter % (maxiter / 10) == 0
+                println("round $iter: change $(diff(G_x_new, G_x))")
+            end
         end
         if maximum(abs.(G_x_new .- G_x)) < tol && iter > 10
             break
@@ -83,26 +85,51 @@ function printG(d, G_x)
     println()
 end
 
-printstyled("=====     Symmetrized DLR solver for SYK model     =======\n", color = :yellow)
-mix = 0.1
-dsym = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-8, symmetry = :ph) # Initialize DLR object
-G_x_ph = solve_syk_with_fixpoint_iter(dsym, 0.00, mix = mix)
-printG(dsym, G_x_ph)
-
-printstyled("=====     Unsymmetrized DLR solver for SYK model     =======\n", color = :yellow)
-mix = 0.02
-dnone = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-8, symmetry = :none) # Initialize DLR object
-G_x_none = solve_syk_with_fixpoint_iter(dnone, 0.00, mix = mix)
-printG(dnone, G_x_none)
-
-printstyled("=====     Unsymmetrized versus Symmetrized DLR solver    =======\n", color = :yellow)
-@printf("%15s%40s%40s%40s\n", "τ", "sym DLR (interpolated)", "unsym DLR", "difference")
-G_x_interp = tau2tau(dsym, G_x_ph, dnone.τ)
-for i in 1:dnone.size
-    if dnone.τ[i] <= dnone.β / 2
-        @printf("%15.8f%40.15f%40.15f%40.15f\n", dnone.τ[i], real(G_x_interp[i]), real(G_x_none[i]), abs(real(G_x_interp[i] - G_x_none[i])))
-    end
+verbose = false
+
+printstyled("=====    Prepare the expected Green's function of the SYK model     =======\n", color = :yellow)
+dsym_correct = DLRGrid(Euv = 5.0, β = 10000.0, isFermi = true, rtol = 1e-14, symmetry = :ph) # Initialize DLR object
+G_x_correct = solve_syk_with_fixpoint_iter(dsym_correct, 0.00, mix = 0.1, verbose = false)
+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")
+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
+    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)
+    # printG(dnone, G_x_none)
+
+    # printstyled("=====     Unsymmetrized versus Symmetrized DLR solver    =======\n", color = :yellow)
+    # @printf("%15s%40s%40s%40s\n", "τ", "sym DLR (interpolated)", "unsym DLR", "difference")
+    # G_x_interp = tau2tau(dsym_correct, G_x_correct, dnone.τ)
+    # for i in 1:dnone.size
+    #     if dnone.τ[i] <= dnone.β / 2
+    #         @printf("%15.8f%40.15f%40.15f%40.15f\n", dnone.τ[i], real(G_x_interp[i]), real(G_x_none[i]), abs(real(G_x_interp[i] - G_x_none[i])))
+    #     end
+    # end
+
+    G_x_interp_ph = tau2tau(dsym_correct, G_x_correct, dsym.τ)
+    G_x_interp_none = 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"
+
+    @printf("%30.15f%30.15e%30.15e%20s\n", Euv, d_ph, d_none, 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))
+
 end
-println("maximumal difference: ", diff(G_x_interp, G_x_none))