Refactor input.jl and plot scripts

.gitignore

@@ -23,6 +23,7 @@ deps/src/
 *.aux
 *.DS_Store
 *.xdv
+*.xcf
 *.dat
 *.err
 *.out

example/freeE/plot_freeE.jl

@@ -12,7 +12,7 @@ Fs = [-0.0,]
 # beta = [25.0]
 beta = [80.0]
 order = [5,]
-const fileName = "freeE_$(dim)d.dat"
+const filename = "freeE_$(dim)d.dat"
 const val_benchmark =
     cdict = Dict(["blue" => "#0077BB", "cyan" => "#33BBEE", "teal" => "#009988", "orange" => "#EE7733", "red" => "#CC3311", "magenta" => "#EE3377", "grey" => "#BBBBBB"]);
 
@@ -53,7 +53,7 @@ end
 
 if abspath(PROGRAM_FILE) == @__FILE__
 
-    freeE_data = readdlm(fileName)
+    freeE_data = readdlm(filename)
     num_data = size(freeE_data)[1]
 
     freeE_total, error_total, mass2_total = [], [], []

example/input.jl

@@ -1,24 +1,90 @@
-dim = 2. # dimension of the problem
-rs = [2.0,]
-# mass2 = [1.0, 2.0, 3.0]
-mass2 = [0.1,]  # screening parameter
-Fs = [-0.0,]    # Fermi liquid parameter with zero angular momentum
-beta = [25.0]   # inverse temperature beta = β*E_F 
-order = [3,]    # order of diagrams
-neval = 1e6    # number of Monte Carlo samples
-# neval = 1e8
+dim = 3 # dimension of the problem
+
+### rs = 0.5 ###
+# rs = [0.5]
+# mass2 = [3.5]  # screening parameter
+
+### rs = 1 ###
+rs = [1.0]
+mass2 = [1.75]  # screening parameter
+
+### rs = 2 ###
+# rs = [2.0]
+# mass2 = [2.0]  # screening parameter
+
+### rs = 3 ###
+# rs = [3.0]
+# mass2 = [1.5]  # screening parameter
+
+### rs = 4 ###
+# rs = [4.0]
+# mass2 = [1.125]  # screening parameter
+
+### rs = 5 ###
+# rs = [5.0]
+# mass2 = [1.125]  # screening parameter
+
+### rs = 6 ###
+# rs = [6.0]
+# mass2 = [0.75]  # screening parameter
+
+Fs = [-0.0]    # Fermi liquid parameter with zero angular momentum
+beta = [40.0]   # inverse temperature beta = β*E_F 
+order = [5]    # order of diagrams
+neval = 1e11    # number of Monte Carlo samples
 isDynamic = false # whether to use effective field theory with dynamic screening or not 
 isFock = false # whether to use Fock renormalization or not
 
 diagGenerate = :GV # :GV or :Parquet, algorithm to generate diagrams
-# diagGenerate = :Parquet
-isLayered2D = true # whether to use layered 2D system or not
+isLayered2D = false # whether to use layered 2D system or not
 
 spin = 2    # 2 for unpolarized, 1 for polarized
 spinPolarPara = 2 / spin - 1 # spin-polarization parameter (n_up - n_down) / (n_up + n_down) ∈ [0,1]
+ispolarized = spinPolarPara != 0.0
+
+# Build file base names
+basenames = [
+    "para_wn_1minus0",
+    "data$(dim)d_Z",
+    "data$(dim)d_K",
+    "meff_$(dim)d",
+    "zfactor_$(dim)d",
+    "inverse_zfactor_$(dim)d",
+    "dispersion_ratio_$(dim)d",
+    "inverse_dispersion_ratio_$(dim)d",
+]
+for i in eachindex(basenames)
+    if spin != 2
+        basenames[i] *= "_spin$(spin)"
+    end
+    if ispolarized
+        basenames[i] *= "_polarized"
+    end
+    if isLayered2D
+        @assert dim == 2 "Layered 2D mode is only available for dim = 2!"
+        basenames[i] *= "_layered2d"
+    end
+end
+para_basename,
+sigma_z_basename,
+sigma_k_basename,
+meff_basename,
+zfactor_basename,
+inverse_zfactor_basename,
+dispersion_ratio_basename,
+inverse_dispersion_ratio_basename = basenames
+
+# Directory paths
+data_directory = joinpath(@__DIR__, "sigma/data$(dim)d")
+res_directory = joinpath(@__DIR__, "sigma")
+para_directory = ""  # src directory
 
-if isLayered2D
-    const parafilename = "para_wn_1minus0_layered2d.csv"
-else
-    const parafilename = "para_wn_1minus0.csv"
-end
+# File paths
+const parafilename = joinpath(para_directory, para_basename * ".csv")
+const simga_z_filename = joinpath(data_directory, sigma_z_basename * ".jld2")
+const sigma_k_filename = joinpath(data_directory, sigma_k_basename * ".jld2")
+const meff_filename = joinpath(res_directory, meff_basename * ".dat")
+const zfactor_filename = joinpath(res_directory, zfactor_basename * ".dat")
+const zinv_filename = joinpath(res_directory, inverse_zfactor_basename * ".dat")
+const dispersion_ratio_filename = joinpath(res_directory, dispersion_ratio_basename * ".dat")
+const inverse_dispersion_ratio_filename = joinpath(res_directory, inverse_dispersion_ratio_basename * ".dat")

example/sigma/derive_dispersion_ratio.jl

@@ -0,0 +1,36 @@
+using UniElectronGas, ElectronLiquid
+using JLD2, DelimitedFiles
+
+include("../input.jl")
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    isSave = false
+    if length(ARGS) >= 1 && (ARGS[1] == "s" || ARGS[1] == "-s" || ARGS[1] == "--save" || ARGS[1] == " save")
+        # the second parameter may be set to save the derived parameters
+        isSave = true
+    end
+
+    f = jldopen(sigma_k_filename, "r")
+    results = Any[]
+    for (_rs, _mass2, _F, _beta, _order) in Iterators.product(rs, mass2, Fs, beta, order)
+        para = ParaMC(rs=_rs, beta=_beta, Fs=_F, order=_order, mass2=_mass2, isDynamic=isDynamic, dim=dim, spin=spin)
+        kF = para.kF
+        for key in keys(f)
+            loadpara = ParaMC(key)
+            if UEG.paraid(loadpara) == UEG.paraid(para)
+                println(UEG.paraid(para))
+
+                ngrid, kgrid, rSw_k, iSw_k = UniElectronGas.getSigma(para, sigma_k_filename; parafile=parafilename)
+                dispersion_ratio, fit_p = UniElectronGas.getDispersionRatio(para, rSw_k, kgrid; parafile=parafilename)
+                println("ϵ_qp / ϵ_0 = ", dispersion_ratio)
+                push!(results, Any[_rs, _beta, _mass2, _order, dispersion_ratio...])
+                writedlm("fit_dispersion.dat", fit_p)
+            end
+        end
+    end
+    if isSave
+        open(dispersion_ratio_filename, "a+") do io
+            writedlm(io, results)
+        end
+    end
+end

example/sigma/derive_inverse_dispersion_ratio.jl

@@ -0,0 +1,36 @@
+using UniElectronGas, ElectronLiquid
+using JLD2, DelimitedFiles
+
+include("../input.jl")
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    isSave = false
+    if length(ARGS) >= 1 && (ARGS[1] == "s" || ARGS[1] == "-s" || ARGS[1] == "--save" || ARGS[1] == " save")
+        # the second parameter may be set to save the derived parameters
+        isSave = true
+    end
+
+    f = jldopen(sigma_k_filename, "r")
+    results = Any[]
+    for (_rs, _mass2, _F, _beta, _order) in Iterators.product(rs, mass2, Fs, beta, order)
+        para = ParaMC(rs=_rs, beta=_beta, Fs=_F, order=_order, mass2=_mass2, isDynamic=isDynamic, dim=dim, spin=spin)
+        kF = para.kF
+        for key in keys(f)
+            loadpara = ParaMC(key)
+            if UEG.paraid(loadpara) == UEG.paraid(para)
+                println(UEG.paraid(para))
+
+                ngrid, kgrid, rSw_k, iSw_k = UniElectronGas.getSigma(para, sigma_k_filename; parafile=parafilename)
+                dispersion_ratio, fit_p = UniElectronGas.getDispersionRatioInv(para, rSw_k, kgrid; parafile=parafilename)
+                println("ϵ_qp / ϵ_0 = ", dispersion_ratio)
+                push!(results, Any[_rs, _beta, _mass2, _order, dispersion_ratio...])
+                writedlm("fit_inverse_dispersion.dat", fit_p)
+            end
+        end
+    end
+    if isSave
+        open(inverse_dispersion_ratio_filename, "a+") do io
+            writedlm(io, results)
+        end
+    end
+end

example/sigma/derive_meff.jl

@@ -3,37 +3,14 @@ using JLD2, DelimitedFiles
 
 include("../input.jl")
 
-# # dim = 3
-# dim = 2
-# spin = 1
-# # spin = 2
-# # rs = [0.5, 1.0, 2.0]
-# rs = [1.0]
-# # rs = [0.1, 0.3]
-# # mass2 = [0.5, 1.0, 1.5, 2.0, 2.2, 2.4, 2.5, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0]
-# # mass2 = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5]
-# mass2 = [1.75]
-# Fs = [-0.0,]
-# # beta = [25.0,]
-# beta = [40.0,]
-# order = [4,]
-# # order = [5,]
-# isDynamic = false
-### dSigma/dk = (Sigma[kF_label+idx_dk] - Sigma[kF_label-idx_dk]) / (kgrid[kF_label+idx_dk] - kgrid[kF_label-idx_dk])
-
-# const parafilename = "para_wn_1minus0.csv"
-const filename = "./data$(dim)d/data$(dim)d_K.jld2"
-# const filename = "./data$(dim)d/data$(dim)d_K_o5.jld2"
-const savefilename = spin == 2 ? "meff_$(dim)d.dat" : "meff_$(dim)d_spin$spin.dat"
-
 if abspath(PROGRAM_FILE) == @__FILE__
     isSave = false
     if length(ARGS) >= 1 && (ARGS[1] == "s" || ARGS[1] == "-s" || ARGS[1] == "--save" || ARGS[1] == " save")
         # the second parameter may be set to save the derived parameters
         isSave = true
     end
 
-    f = jldopen(filename, "r")
+    f = jldopen(sigma_k_filename, "r")
     results = Any[]
     for (_rs, _mass2, _F, _beta, _order) in Iterators.product(rs, mass2, Fs, beta, order)
         para = ParaMC(rs=_rs, beta=_beta, Fs=_F, order=_order, mass2=_mass2, isDynamic=isDynamic, dim=dim, spin=spin)
@@ -43,8 +20,8 @@ if abspath(PROGRAM_FILE) == @__FILE__
             if UEG.paraid(loadpara) == UEG.paraid(para)
                 println(UEG.paraid(para))
 
-                ngrid, kgrid, rSw_k, iSw_k = UniElectronGas.getSigma(para, filename)
-                meff, fit_p = UniElectronGas.getMeff(para, rSw_k, kgrid)
+                ngrid, kgrid, rSw_k, iSw_k = UniElectronGas.getSigma(para, sigma_k_filename; parafile=parafilename)
+                meff, fit_p = UniElectronGas.getMeff(para, rSw_k, kgrid; parafile=parafilename)
                 # meff = UniElectronGas.getMeff(para, rSw_k, kgrid)
                 println("m* / m = ", meff)
                 push!(results, Any[_rs, _beta, _mass2, _order, meff...])
@@ -53,8 +30,8 @@ if abspath(PROGRAM_FILE) == @__FILE__
         end
     end
     if isSave
-        open(savefilename, "a+") do io
+        open(meff_filename, "a+") do io
             writedlm(io, results)
         end
     end
-end
+end

example/sigma/derive_zinv.jl

@@ -0,0 +1,35 @@
+using UniElectronGas, ElectronLiquid
+using JLD2, DelimitedFiles
+
+include("../input.jl")
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    isSave = false
+    if length(ARGS) >= 1 && (ARGS[1] == "s" || ARGS[1] == "-s" || ARGS[1] == "--save" || ARGS[1] == " save")
+        # the second parameter may be set to save the derived parameters
+        isSave = true
+    end
+
+    f = jldopen(sigma_k_filename, "r")
+    results = Any[]
+    for (_rs, _mass2, _F, _beta, _order) in Iterators.product(rs, mass2, Fs, beta, order)
+        para = ParaMC(rs=_rs, beta=_beta, Fs=_F, order=_order, mass2=_mass2, isDynamic=isDynamic, dim=dim, spin=spin)
+        kF = para.kF
+        for key in keys(f)
+            loadpara = ParaMC(key)
+            if UEG.paraid(loadpara) == UEG.paraid(para)
+                println(UEG.paraid(para))
+
+                ngrid, kgrid, rSw_k, iSw_k = UniElectronGas.getSigma(para, sigma_k_filename; parafile=parafilename)
+                zinv = UniElectronGas.getZinv(para; parafile=parafilename)
+                println("1/z = ", zinv)
+                push!(results, Any[_rs, _beta, _mass2, _order, zinv...])
+            end
+        end
+    end
+    if isSave
+        open(zinv_filename, "a+") do io
+            writedlm(io, results)
+        end
+    end
+end

example/sigma/derive_zmu.jl

@@ -3,16 +3,6 @@ using JLD2, DelimitedFiles
 
 include("../input.jl")
 
-if isLayered2D
-    const filename = "./data$(dim)d/data_Z_layered2d.jld2"
-else
-    const filename = "./data$(dim)d/data_Z.jld2"
-end
-# const filename = "./data2d_Z_v0.jld2"
-# const filename = "./data$(dim)d/data$(dim)d_Z_beta80_rs$(rs[1]).jld2"
-const savefilename = spin == 2 ? "zfactor_$(dim)d.dat" : "zfactor_$(dim)d_spin$spin.dat"
-
-
 function zfactor_renorm(dz, dzinv; isRenorm=true)
     if isRenorm
         sumzinv = accumulate(+, dzinv)
@@ -38,8 +28,8 @@ if abspath(PROGRAM_FILE) == @__FILE__
         isSave = true
     end
 
-    f = jldopen(filename, "r")
-    println(filename)
+    f = jldopen(sigma_z_filename, "r")
+    println(sigma_z_filename)
     results = Any[]
     for (_rs, _mass2, _F, _beta, _order) in Iterators.product(rs, mass2, Fs, beta, order)
         para = ParaMC(rs=_rs, beta=_beta, Fs=_F, order=_order, mass2=_mass2, isDynamic=isDynamic, dim=dim, spin=spin)
@@ -54,8 +44,8 @@ if abspath(PROGRAM_FILE) == @__FILE__
         end
     end
     if isSave
-        open(savefilename, "a+") do io
+        open(zfactor_filename, "a+") do io
             writedlm(io, results)
         end
     end
-end
+end