Skip to content

Commit 4819eff

Browse files
LittleBugLittleBug
committed
improve SYK.jl demo
1 parent 4f9208b commit 4819eff

File tree

2 files changed

+99
-89
lines changed

2 files changed

+99
-89
lines changed

‎example/SYK.jl

Lines changed: 69 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -1,88 +1,98 @@
11
"""
2-
Demonstration of solution of SYK equation by Matsubara frequency/imaginary time
3-
alternation method using the discrete Lehmann representation
4-
5-
The SYK equation is the nonlinear Dyson equation corresponding to self-energy
6-
7-
Σ(τ) = J^2 * G(τ) * G(τ) * G(β-τ).
8-
9-
We solve the Dyson equation self-consistently by a weighted
10-
fixed point iteration, with weight w assigned to the new iterate
11-
and weight 1-w assigned to the previous iterate. The self-energy
12-
is evaluated in the imaginary time domain, and each linear Dyson
13-
equation, corresponding to fixed self-energy, is solved in the
14-
Matsubara frequency domain, where it is diagonal.
15-
16-
To solve the equation with a desired chemical potential μ,
17-
we pick a number nmu>1, and solve a sequence intermediate
18-
problems to obtain a good initial guess. First we solve the
19-
equation with chemical potential μ_0 = 0, then use this
20-
solution as an initial guess for the fixed point iteration with
21-
μ_1 = μ/nmu, then use this the solution as an initial guess
22-
for the fixed point iteration with μ = 2*μ/nmu, and so on,
23-
until we reach μ_{nmu} = μ.
2+
A SYK model solver based on a forward fixed-point iteration method.
3+
4+
The self-energy of the SYK model is given by,
5+
6+
Σ(τ) = J^2 * G(τ) * G(τ) * G(β-τ),
7+
8+
where Green's function of the SYK model is given by the Dyson equation,
9+
10+
G(iωₙ) = -1/(iωₙ -μ + Σ(iωₙ))
11+
12+
We solve the Dyson equation self-consistently by a weighted fixed point iteration,
13+
with weight `mix` assigned to the new iterate and weight `1-mix` assigned to the previous iterate.
14+
15+
The self-energy is evaluated in the imaginary time domain,
16+
and the Dyson equation is solved in the Matsubara frequency domain.
17+
18+
The SYK Green's function has particle-hole symmetry when μ=0.
19+
You may enforce such symmetry by setting `symmetry = :ph` when initialize the DLR grids.
20+
A symmetrized solver tends to be more robust than the unsymmetrized counterpart.
2421
"""
2522

2623
using Lehmann
24+
using Printf
2725

2826
diff(a, b) = maximum(abs.(a - b)) # return the maximum deviation between a and b
2927

30-
conformal_tau(τ, β) = -π^(1 / 4) / sqrt(2β) * 1 / sqrt(sin(π * τ / β))
28+
conformal_tau(τ, β) = π^(1 / 4) / sqrt(2β) * 1 / sqrt(sin(π * τ / β))
3129

3230
function syk_sigma_dlr(d, G_x, J = 1.0)
3331

3432
tau_k = d.τ # DLR imaginary time nodes
3533
tau_k_rev = d.β .- tau_k # Reversed imaginary time nodes
3634

37-
G_k = dlr2tau(d, G_x) # G at tau_k
38-
G_k_rev = dlr2tau(d, G_x, tau_k_rev) # G at beta - tau_k
35+
G_x_rev = tau2tau(d, G_x, tau_k_rev) # G at beta - tau_k
3936

40-
# for i in 1:length(G_k)
41-
# println("$(d.τ[i]) $(real(G_k[i])) $(imag(G_k[i]))")
42-
# end
37+
Sigma_x = J .^ 2 .* G_x .^ 2 .* G_x_rev # SYK self-energy in imaginary time
4338

44-
# for i = 1:length(tau_k)
45-
# G_k[i] = conformal_tau(tau_k[i], d.β)
46-
# G_k_rev[i] = conformal_tau(tau_k_rev[i], d.β)
47-
# end
48-
49-
Sigma_k = J .^ 2 .* G_k .^ 2 .* G_k_rev # SYK self-energy in imaginary time
50-
Sigma_x = tau2dlr(d, Sigma_k) # DLR coeffs of self-energy
51-
52-
# println("sigma diff: ", diff(Sigma_k, dlr2tau(d, Sigma_x)))
53-
# for i in 1:length(Sigma_x)
54-
# println("$(d.τ[i]) $(real(Sigma_x[i])) $(imag(Sigma_x[i]))")
55-
# end
56-
# exit(0)
5739
return Sigma_x
5840
end
5941

60-
function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol, mix = 0.3, maxiter = 100)
42+
function dyson(d, sigma_q, mu)
43+
if d.symmetry == :ph #symmetrized G
44+
@assert mu ≈ 0.0 "Only μ=0 case enjoys the particle-hole symmetry."
45+
return 1im * imag.(1 ./ (d.ωn * 1im .+ mu .- sigma_q))
46+
elseif d.symmetry == :none
47+
return -1 ./ (d.ωn * 1im .- mu .+ sigma_q)
48+
else
49+
error("Not implemented!")
50+
end
51+
end
52+
53+
function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol * 10; mix = 0.1, maxiter = 1000, G_x = zeros(ComplexF64, length(d)))
6154

62-
Sigma_q = zeros(length(d)) # Initial guess
63-
G_q = zeros(ComplexF64, length(d))
6455
for iter in 1:maxiter
6556

66-
G_q .= -1 ./ (d.ωn * 1im .- mu .+ Sigma_q) # Solve Dyson
67-
G_x = matfreq2dlr(d, G_q) # Get DLR coeffs
68-
Sigma_x_new = syk_sigma_dlr(d, G_x)
69-
Sigma_q_new = dlr2matfreq(d, Sigma_x_new)
57+
Sigma_x = syk_sigma_dlr(d, G_x)
7058

71-
# println(diff(Sigma_q_new, Sigma_q))
72-
if maximum(abs.(Sigma_q_new .- Sigma_q)) < tol
59+
G_q_new = dyson(d, tau2matfreq(d, Sigma_x), mu)
60+
61+
G_x_new = matfreq2tau(d, G_q_new)
62+
63+
64+
if iter % (maxiter / 10) == 0
65+
println("round $iter: change $(diff(G_x_new, G_x))")
66+
end
67+
if maximum(abs.(G_x_new .- G_x)) < tol && iter > 10
7368
break
7469
end
75-
Sigma_q = mix * Sigma_q_new + (1 - mix) * Sigma_q # Linear mixing
70+
71+
G_x = mix * G_x_new + (1 - mix) * G_x # Linear mixing
7672
end
77-
return G_q
73+
return G_x
7874
end
7975

76+
function printG(d, G_x)
77+
@printf("%15s%40s%40s%40s\n", "τ", "DLR imag", "DLR real", "asymtotically exact")
78+
for i in 1:d.size
79+
if d.τ[i] <= d.β / 2
80+
@printf("%15.8f%40.15f%40.15f%40.15f\n", d.τ[i], imag(G_x[i]), real(G_x[i]), conformal_tau(d.τ[i], d.β))
81+
end
82+
end
83+
println()
84+
end
85+
86+
printstyled("===== Symmetrized DLR solver for SYK model =======\n", color = :yellow)
87+
mix = 0.1
88+
d = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-6, symmetry = :ph) # Initialize DLR object
89+
G_x = solve_syk_with_fixpoint_iter(d, 0.00, mix = mix)
90+
printG(d, G_x)
91+
92+
printstyled("===== Unsymmetrized DLR solver for SYK model =======\n", color = :yellow)
93+
mix = 0.01
94+
d = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-6, symmetry = :none) # Initialize DLR object
95+
G_x = solve_syk_with_fixpoint_iter(d, 0.00, mix = mix)
96+
printG(d, G_x)
8097

81-
d = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-14) # Initialize DLR object
82-
G_q = solve_syk_with_fixpoint_iter(d, 0.0)
8398

84-
G_q = matfreq2tau(d, G_q)
85-
println(G_q)
86-
for i in 1:length(G_q)
87-
println("$(d.τ[i]) $(real(G_q[i])) $(imag(G_q[i]))")
88-
end

‎example/demo3.jl renamed to ‎example/SYK_backup.jl

Lines changed: 30 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -34,55 +34,55 @@ function syk_sigma_dlr(d, G_x, J = 1.0)
3434
tau_k = d.τ # DLR imaginary time nodes
3535
tau_k_rev = d.β .- tau_k # Reversed imaginary time nodes
3636

37-
G_x_rev = tau2tau(d, G_x, tau_k_rev) # G at beta - tau_k
37+
G_k = dlr2tau(d, G_x) # G at tau_k
38+
G_k_rev = dlr2tau(d, G_x, tau_k_rev) # G at beta - tau_k
3839

39-
Sigma_x = J .^ 2 .* G_x .^ 2 .* G_x_rev # SYK self-energy in imaginary time
40+
# for i in 1:length(G_k)
41+
# println("$(d.τ[i]) $(real(G_k[i])) $(imag(G_k[i]))")
42+
# end
43+
44+
# for i = 1:length(tau_k)
45+
# G_k[i] = conformal_tau(tau_k[i], d.β)
46+
# G_k_rev[i] = conformal_tau(tau_k_rev[i], d.β)
47+
# end
48+
49+
Sigma_k = J .^ 2 .* G_k .^ 2 .* G_k_rev # SYK self-energy in imaginary time
50+
Sigma_x = tau2dlr(d, Sigma_k) # DLR coeffs of self-energy
4051

4152
# println("sigma diff: ", diff(Sigma_k, dlr2tau(d, Sigma_x)))
4253
# for i in 1:length(Sigma_x)
4354
# println("$(d.τ[i]) $(real(Sigma_x[i])) $(imag(Sigma_x[i]))")
4455
# end
4556
# exit(0)
46-
# return real.(Sigma_k)
4757
return Sigma_x
4858
end
4959

50-
function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol; mix = 0.2, maxiter = 1000, G_x = zeros(ComplexF64, length(d)))
60+
function solve_syk_with_fixpoint_iter(d, mu, tol = d.rtol, mix = 0.3, maxiter = 100)
5161

62+
Sigma_q = zeros(length(d)) # Initial guess
63+
G_q = zeros(ComplexF64, length(d))
5264
for iter in 1:maxiter
5365

54-
# println("G diff: ", diff(G_q, dlr2matfreq(d, G_x)))
55-
Sigma_x = syk_sigma_dlr(d, G_x)
56-
57-
# for i in 1:d.size
58-
# println("$(d.τ[i]) $(real(Sigma_x_new[i])) $(imag(Sigma_x_new[i]))")
59-
# end
60-
# println(typeof(Sigma_x_new))
61-
62-
# G_q_new = -1 ./ (d.ωn * 1im .- mu .+ tau2matfreq(d, Sigma_x))
63-
# +1 ./ (-d.ωn * 1im .- mu .+ tau2matfreq(d, Sigma_x)) # Solve Dyson
64-
G_q_new = -1 ./ (d.ωn * 1im .- mu .+ tau2matfreq(d, Sigma_x)) # Solve Dyson
66+
G_q .= -1 ./ (d.ωn * 1im .- mu .+ Sigma_q) # Solve Dyson
67+
G_x = matfreq2dlr(d, G_q) # Get DLR coeffs
68+
Sigma_x_new = syk_sigma_dlr(d, G_x)
69+
Sigma_q_new = dlr2matfreq(d, Sigma_x_new)
6570

66-
G_x_new = matfreq2tau(d, G_q_new)
67-
68-
println("imag", maximum(abs.(imag.(G_x_new))))
69-
70-
println(diff(G_x_new, G_x))
71-
if maximum(abs.(G_x_new .- G_x)) < tol && iter > 5
71+
# println(diff(Sigma_q_new, Sigma_q))
72+
if maximum(abs.(Sigma_q_new .- Sigma_q)) < tol
7273
break
7374
end
74-
75-
G_x = mix * G_x_new + (1 - mix) * G_x # Linear mixing
75+
Sigma_q = mix * Sigma_q_new + (1 - mix) * Sigma_q # Linear mixing
7676
end
77-
return G_x
77+
return G_q
7878
end
7979

80-
d = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-7, symmetry = :ph) # Initialize DLR object
81-
G_q = solve_syk_with_fixpoint_iter(d, 0.00)
82-
# G_q = solve_syk_with_fixpoint_iter(d, 0.15, G_x = G_q)
83-
# G_q = solve_syk_with_fixpoint_iter(d, 0.12, G_x = G_q)
84-
# G_q = solve_syk_with_fixpoint_iter(d, 0.11, G_x = G_q)
85-
# println(G_q)
80+
81+
d = DLRGrid(Euv = 5.0, β = 1000.0, isFermi = true, rtol = 1e-14) # Initialize DLR object
82+
G_q = solve_syk_with_fixpoint_iter(d, 0.0)
83+
84+
G_q = matfreq2tau(d, G_q)
85+
println(G_q)
8686
for i in 1:length(G_q)
8787
println("$(d.τ[i]) $(real(G_q[i])) $(imag(G_q[i]))")
8888
end

0 commit comments

Comments
 (0)