Skip to content

Commit

Permalink
Cleanup GW scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
@dcerkoney
dcerkoney committed Feb 26, 2024
1 parent e978c7f commit f6cb077
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 60 deletions.
64 changes: 15 additions & 49 deletions example/sigma/get_meff_gw.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
32 changes: 21 additions & 11 deletions example/sigma/plot_sigma_gw.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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"])
Expand All @@ -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
Expand Down Expand Up @@ -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")
Expand All @@ -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()

0 comments on commit f6cb077

Please sign in to comment.