diff --git a/example/sigma/get_meff_gw.jl b/example/sigma/get_meff_gw.jl index 748ca5b..b2c2dc1 100644 --- a/example/sigma/get_meff_gw.jl +++ b/example/sigma/get_meff_gw.jl @@ -1,37 +1,11 @@ +using CompositeGrids +using ElectronGas +using ElectronLiquid using JLD2 using PyCall -using PyPlot -using ElectronGas, ElectronLiquid, Parameters -using Lehmann, GreenFunc, CompositeGrids @pyimport numpy as np # for saving/loading numpy data -# Vibrant qualitative colour scheme from https://personal.sron.nl/~pault/ -const cdict = Dict([ - "orange" => "#EE7733", - "blue" => "#0077BB", - "cyan" => "#33BBEE", - "magenta" => "#EE3377", - "red" => "#CC3311", - "teal" => "#009988", - "grey" => "#BBBBBB", -]); -style = PyPlot.matplotlib."style" -style.use(["science", "std-colors"]) -const color = [ - "black", - cdict["orange"], - cdict["blue"], - cdict["cyan"], - cdict["magenta"], - cdict["red"], - cdict["teal"], -] -rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams") -rcParams["font.size"] = 16 -rcParams["mathtext.fontset"] = "cm" -rcParams["font.family"] = "Times New Roman" - @inline function get_Fs(rs) return get_Fs_PW(rs) # if rs < 1.0 || rs > 5.0 @@ -74,7 +48,7 @@ susceptibility ratio data (c.f. Kukkonen & Chen, 2021) end function GW0(G_prev, param, Euv, rtol, Nk, maxK, minK, order, int_type, kgrid::Union{AbstractGrid,AbstractVector,Nothing}=nothing; kwargs...) - @unpack dim = param + dim = param.dim if dim == 2 if isnothing(kgrid) @@ -122,15 +96,15 @@ function get_meff_from_Σ_GW(param::Parameter.Para; int_type=:rpa, δK=5e-6, max end # Helper function for linear interpolation with mixing parameter α -function LERP(M_start, M_end) +function LERP(M_start, M_end, alpha) return (1 - alpha) * M_start + alpha * M_end end function get_Σ_GW(param::Parameter.Para, int_type, max_steps, atol, alpha, δK, save_sigma) # Make sigma output directory if needed if save_sigma - mkpath("sigma_GW_$(param.dim)d") - dir = joinpath(@__DIR__, "sigma_GW_$(param.dim)d") + mkpath("results/sigma_GW_$(param.dim)d") + dir = joinpath(@__DIR__, "results/sigma_GW_$(param.dim)d") end # Make sure we are using parameters for the bare UEG theory @@ -187,9 +161,9 @@ function get_Σ_GW(param::Parameter.Para, int_type, max_steps, atol, alpha, δK, zfactor_prev = SelfEnergy.zfactor(param, Σ; ngrid=[-1, 0])[1] @assert kamp ≈ param.kF - # Write initial (G0W0) self-energy to JLD2 file + # Write initial (G0W0) self-energy to JLD2 file, overwriting if it already exists if save_sigma - jldopen(joinpath(dir, "sigma_$(int_type)_rs$(rs).jl"), "a+") do file + jldopen(joinpath(dir, "sigma_$(int_type)_rs$(rs).jl"), "w") do file file["Σ_0"] = Σ file["Σ_ins_0"] = Σ_ins end @@ -216,12 +190,12 @@ function get_Σ_GW(param::Parameter.Para, int_type, max_steps, atol, alpha, δK, # Get Σ[G_mix, W], where the input Green's function is a linear interpolation # of the previous and current G: G_mix = (1 - α) * G_prev + α * G - G_mix = LERP(G_prev, G) + G_mix = LERP(G_prev, G, alpha) Σ, Σ_ins = Σ_GW0(G_mix) - # Write self-energy at this step to JLD2 file + # Append self-energy at this step to JLD2 file if save_sigma - jldopen(joinpath(dir, "sigma_$(int_type).jl"), "a+") do file + jldopen(joinpath(dir, "sigma_$(int_type)_rs$(rs).jl"), "a+") do file file["Σ_$(i_step + 1)"] = Σ file["Σ_ins_$(i_step + 1)"] = Σ_ins end @@ -316,23 +290,15 @@ function main() save_sigma = true # Output directory - mkpath("finalized_meff_results/$(dim)d") - dir = joinpath(@__DIR__, "finalized_meff_results/$(dim)d") - - # rslist = [2.0] - # rslist = [1.0, 5.0, 10.0] + mkpath("results/finalized_meff_results/$(dim)d") + dir = joinpath(@__DIR__, "results/finalized_meff_results/$(dim)d") # rslist = [0.001; collect(LinRange(0.0, 1.1, 111))[2:end]] # for accurate 2D HDL # rslist = [0.005; collect(LinRange(0.0, 5.0, 101))[2:end]] # for 2D # rslist = [0.01; collect(LinRange(0.0, 10.0, 101))[2:end]] # for 3D # rslist = [0.001; collect(range(0.0, 1.1, step=0.05))[2:end]] # for accurate 2D HDL - # rslist = [0.01; collect(range(0.0, 10.0, step=0.5))[2:end]] # for 3D - - # rslist = [0.001; collect(range(0.0, 1.1, step=0.1))[2:end]] # for accurate 2D HDL - # rslist = [0.01; collect(range(0.0, 10.0, step=1.0))[2:end]] # for 3D - - rslist = [10.0] + rslist = [0.01; collect(range(0.0, 10.0, step=0.5))[2:end]] # for 3D # NOTE: int_type ∈ [:ko_const, :ko_takada_plus, :ko_takada, :ko_moroni, :ko_simion_giuliani] # NOTE: KO interaction using G+ and/or G- is currently only available in 3D diff --git a/example/sigma/plot_sigma_gw.jl b/example/sigma/plot_sigma_gw.jl index 0597afd..5863e97 100644 --- a/example/sigma/plot_sigma_gw.jl +++ b/example/sigma/plot_sigma_gw.jl @@ -73,7 +73,9 @@ end function main() # System parameters dim = 3 - rs = 2.0 + # rs = 0.01 + # rs = 2.0 + rs = 10.0 δK = 5e-6 beta = 1000.0 param = Parameter.rydbergUnit(1.0 / beta, rs, dim) @@ -92,7 +94,7 @@ function main() s = [] si = [] local n_max - jldopen("sigma_GW_$(dim)d/sigma_$(int_type).jl", "r") do f + jldopen("results/sigma_GW_$(dim)d/sigma_$(int_type)_rs$(rs).jl", "r") do f for i in 0:max_steps try push!(s, f["Σ_$i"]) @@ -116,7 +118,8 @@ function main() push!(s_iw0s, s_iw[idx_iw0, :]) push!(si_iw0s, si_iw[idx_iw0, :]) end - + println(length(s)) + # DLR cutoffs kF = param.kF Euv, rtol = 1000 * param.EF, 1e-11 @@ -165,43 +168,50 @@ function main() end ax[1].annotate( + # "\$r_s = $(rs)\$", "\$r_s = $(Int(rs))\$", - xy=(0.775, 0.225), + # xy=(0.775, 0.225), + # xy=(0.8, 0.225), + xy=(0.786, 0.225), xycoords="axes fraction", ha="center", va="center", ) ax[2].annotate( "\$\\left(\\frac{m^*}{m}\\right)_{G_0 $wstring} = $(round(meff_g0w; sigdigits=4))\$", - xy=(0.75, 0.75), + xy=(0.737, 0.7), + # xy=(0.75, 0.75), xycoords="axes fraction", ha="center", va="center", ) ax[2].annotate( "\$\\left(\\frac{m^*}{m}\\right)_{G $wstring} = $(round(meff_gw; sigdigits=4))\$", - xy=(0.757, 0.6), + xy=(0.757, 0.55), + # xy=(0.757, 0.6), xycoords="axes fraction", ha="center", va="center", ) ax[2].annotate( "\$Z_{G_0 $wstring} = $(round(zfactor_g0w; sigdigits=4))\$", - xy=(0.783, 0.45), + xy=(0.783, 0.4), + # xy=(0.783, 0.45), xycoords="axes fraction", ha="center", va="center", ) ax[2].annotate( "\$Z_{G $wstring} = $(round(zfactor_gw; sigdigits=4))\$", - xy=(0.79, 0.3), + xy=(0.79, 0.25), + # xy=(0.79, 0.3), xycoords="axes fraction", ha="center", va="center", ) - ax[1].set_ylim(-0.35, 0.25) - ax[1].set_yticks(-0.3:0.1:0.2) + # ax[1].set_ylim(-0.35, 0.25) + # ax[1].set_yticks(-0.3:0.1:0.2) ax[1].set_xlim(0, 5) ax[2].set_xlim(0, 5) ax[1].legend(loc="best") @@ -210,7 +220,7 @@ function main() ax[2].set_ylabel("\$\\text{Im}\\Sigma_c(k, i\\omega_0)\$") ax[2].set_xlabel("\$k / k_F\$") plt.tight_layout() - fig.savefig("static_gw_self_energy_$(int_type).pdf") + fig.savefig("figures/static_gw_self_energy_$(int_type)_rs$(rs).pdf") end main()