Skip to content

Commit 691e6ea

Browse files
LouhoksesonLouhokseson
committed
wigner distribution illustrate
1 parent 09036c4 commit 691e6ea

File tree

4 files changed

+475
-0
lines changed

4 files changed

+475
-0
lines changed

‎scripts/plots/Wigner_distribution/3D.jl‎

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
1+
using DrWatson
2+
@quickactivate "HokseonModelSimulation"
3+
using Unitful, UnitfulAtomic
4+
using CairoMakie
5+
using HokseonPlots
6+
using ColorSchemes
7+
using Colors
8+
using CSV
9+
using DataFrames
10+
using DelimitedFiles
11+
using KernelDensity, Distributions
12+
using LaTeXStrings
13+
using QuadGK, Interpolations
14+
colorscheme = ColorScheme(parse.(Colorant, ["#045275", "#089099", "#7CCBA2", "#FCDE9C", "#F0746E", "#DC3977", "#7C1D6F"]));
15+
colormap = HokseonPlots.NICECOLORS;
16+
17+
using Glob
18+
## Load the scripts in folder HokseonModelSimulation/src
19+
for file in readdir(srcdir(), join=true) # join=true returns the full path
20+
if occursin(r"\.jl$", file) # Check if the file ends with .jl
21+
include(file)
22+
end
23+
end
24+
25+
26+
is_Wigner = [true] # can be a vector of Bool
27+
28+
all_params = Dict{String, Any}(
29+
"trajectories" => [500],
30+
"nstates" => [150],
31+
"dt" => [0.05],
32+
"width" => [50],
33+
"mass" => [1.00784],
34+
"temperature" => [300.0],
35+
"tmax" => [1001],
36+
"discretisation" => [:GapGaussLegendre],
37+
"impuritymodel" => :Hokseon,
38+
"method" => [:AdiabaticIESH],
39+
"incident_energy" => [6.17],
40+
"couplings_rescale" => [2.5],
41+
"centre" => [0],
42+
"gap" => [0.49],
43+
"decoherence" => [:EDC],
44+
"is_Wigner" => is_Wigner,
45+
# 💡 Elementwise sigma assignment:
46+
"sigma" => [w ? 0.3 : nothing for w in is_Wigner],
47+
)
48+
49+
params_list = dict_list(all_params)
50+
# just make sure that params_list is a list with Dicts
51+
if typeof(params_list) != Vector{Dict{String, Any}}
52+
params_list = [params_list]
53+
end
54+
55+
56+
@unpack mass,incident_energy = params_list[1]
57+
###Initial Conditions
58+
m = ustrip(auconvert(mass*u"u"))
59+
position = austrip(5u"Å")
60+
ke = austrip(incident_energy * u"eV")
61+
velocity = - sqrt(2ke / m)
62+
@unpack is_Wigner = params_list[1]
63+
if is_Wigner
64+
@unpack sigma = params_list[1]
65+
lambda = 1/(sigma*2*m) #in a.u.
66+
# Define Gaussian distributions for position and velocity
67+
x_distribution = Normal(position, sigma) #position
68+
v_distribution = Normal(velocity, lambda) #velocity
69+
@info "Nuclei initialization: Wigner distributed"
70+
end
71+
72+
73+
function threeD_wigner_distribution(x_distribution, v_distribution)
74+
75+
fig = Figure(
76+
size = (800, 600),
77+
fonts = (; regular = "MinionPro-Capt.otf"),
78+
fontsize = 12
79+
)
80+
81+
ax = Axis3(
82+
fig[1, 1],
83+
xlabel = "Position / Å",
84+
ylabel = "Kinetic Energy / eV",
85+
zlabel = "Probability Density / (Å·eV)⁻¹",
86+
xgridvisible = false,
87+
ygridvisible = false,
88+
zgridvisible = false
89+
)
90+
91+
## Position Distribution
92+
μ_position = x_distribution.μ
93+
σ_position = x_distribution.σ
94+
x_range_au = collect(range(μ_position - 4σ_position, μ_position + 4σ_position, length=100))
95+
f_x_au = pdf.(x_distribution, x_range_au)
96+
Hartree_to_SI = 0.529177210903
97+
x_range_SI = ustrip.(auconvert.(u"Å", x_range_au))
98+
f_x_SI = f_x_au / Hartree_to_SI
99+
100+
## Velocity → Kinetic Energy Distribution
101+
μ_velocity = v_distribution.μ
102+
σ_velocity = v_distribution.σ
103+
v_range_au = collect(range(μ_velocity - 4σ_velocity, μ_velocity + 4σ_velocity, length=100))
104+
f_v_au = pdf.(v_distribution, v_range_au)
105+
KE_range_au = 0.5 .* m .* v_range_au .^ 2
106+
KE_range_SI = ustrip.(auconvert.(u"eV", KE_range_au))
107+
f_KE_au = f_v_au ./ (m .* sqrt.(2 .* KE_range_au ./ m))
108+
Hartree_to_SI_KE = 27.211386245988
109+
f_KE_SI = f_KE_au / Hartree_to_SI_KE
110+
111+
## Create 2D surface (outer product)
112+
probability_density = f_x_SI .* f_KE_SI' # 100×100 matrix
113+
114+
surface!(ax, x_range_SI, KE_range_SI, probability_density)
115+
116+
fig
117+
end
118+
119+
120+
threeD_wigner_distribution(x_distribution, v_distribution)

‎scripts/plots/Wigner_distribution/heatmap.jl‎

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
using DrWatson
2+
@quickactivate "HokseonModelSimulation"
3+
using Unitful, UnitfulAtomic
4+
using CairoMakie
5+
using HokseonPlots
6+
using ColorSchemes
7+
using Colors
8+
using CSV
9+
using DataFrames
10+
using DelimitedFiles
11+
using KernelDensity, Distributions
12+
using LaTeXStrings
13+
using QuadGK, Interpolations
14+
colorscheme = ColorScheme(parse.(Colorant, ["#045275", "#089099", "#7CCBA2", "#FCDE9C", "#F0746E", "#DC3977", "#7C1D6F"]));
15+
colormap = HokseonPlots.NICECOLORS;
16+
17+
using Glob
18+
## Load the scripts in folder HokseonModelSimulation/src
19+
for file in readdir(srcdir(), join=true) # join=true returns the full path
20+
if occursin(r"\.jl$", file) # Check if the file ends with .jl
21+
include(file)
22+
end
23+
end
24+
25+
26+
is_Wigner = [true] # can be a vector of Bool
27+
28+
all_params = Dict{String, Any}(
29+
"trajectories" => [500],
30+
"nstates" => [150],
31+
"dt" => [0.05],
32+
"width" => [50],
33+
"mass" => [1.00784],
34+
"temperature" => [300.0],
35+
"tmax" => [1001],
36+
"discretisation" => [:GapGaussLegendre],
37+
"impuritymodel" => :Hokseon,
38+
"method" => [:AdiabaticIESH],
39+
"incident_energy" => [6.17],
40+
"couplings_rescale" => [2.5],
41+
"centre" => [0],
42+
"gap" => [0.49],
43+
"decoherence" => [:EDC],
44+
"is_Wigner" => is_Wigner,
45+
# 💡 Elementwise sigma assignment:
46+
"sigma" => [w ? 0.4 : nothing for w in is_Wigner],
47+
)
48+
49+
params_list = dict_list(all_params)
50+
# just make sure that params_list is a list with Dicts
51+
if typeof(params_list) != Vector{Dict{String, Any}}
52+
params_list = [params_list]
53+
end
54+
55+
56+
@unpack mass,incident_energy = params_list[1]
57+
###Initial Conditions
58+
m = ustrip(auconvert(mass*u"u"))
59+
position = austrip(5u"Å")
60+
ke = austrip(incident_energy * u"eV")
61+
velocity = - sqrt(2ke / m)
62+
@unpack is_Wigner = params_list[1]
63+
if is_Wigner
64+
@unpack sigma = params_list[1]
65+
lambda = 1/(sigma*2*m) #in a.u.
66+
# Define Gaussian distributions for position and velocity
67+
x_distribution = Normal(position, sigma) #position
68+
v_distribution = Normal(velocity, lambda) #velocity
69+
@info "Nuclei initialization: Wigner distributed"
70+
end
71+
72+
73+
74+
function heatmap_wigner_distribution(x_distribution, v_distribution)
75+
76+
fig = Figure(
77+
size = (800, 600),
78+
fonts = (; regular = "MinionPro-Capt.otf"),
79+
fontsize = 17
80+
)
81+
82+
ax = Axis(
83+
fig[1, 1],
84+
xlabel = "Position / Å",
85+
ylabel = "Kinetic Energy / eV",
86+
#zlabel = "Probability Density / (Å·eV)⁻¹",
87+
xgridvisible = false,
88+
ygridvisible = false,
89+
#zgridvisible = false
90+
)
91+
92+
## Position Distribution
93+
μ_position = x_distribution.μ
94+
σ_position = x_distribution.σ
95+
x_range_au = collect(range(μ_position - 4σ_position, μ_position + 4σ_position, length=200))
96+
f_x_au = pdf.(x_distribution, x_range_au)
97+
Hartree_to_SI = 0.529177210903
98+
x_range_SI = ustrip.(auconvert.(u"Å", x_range_au))
99+
f_x_SI = f_x_au / Hartree_to_SI
100+
101+
## Velocity → Kinetic Energy Distribution
102+
μ_velocity = v_distribution.μ
103+
σ_velocity = v_distribution.σ
104+
v_range_au = collect(range(μ_velocity - 4σ_velocity, μ_velocity + 4σ_velocity, length=200))
105+
f_v_au = pdf.(v_distribution, v_range_au)
106+
KE_range_au = 0.5 .* m .* v_range_au .^ 2
107+
KE_range_SI = ustrip.(auconvert.(u"eV", KE_range_au))
108+
f_KE_au = f_v_au ./ (m .* sqrt.(2 .* KE_range_au ./ m))
109+
Hartree_to_SI_KE = 27.211386245988
110+
f_KE_SI = f_KE_au / Hartree_to_SI_KE
111+
112+
## Create 2D surface (outer product)
113+
probability_density = f_x_SI .* f_KE_SI' # 100×100 matrix
114+
115+
#surface!(ax, x_range_SI, KE_range_SI, probability_density)
116+
117+
hm = heatmap!(ax, x_range_SI, KE_range_SI, probability_density)
118+
119+
Colorbar(fig[1, 2], hm, label = "Probability Density / (Å·eV)⁻¹", labelrotation = -π/2)
120+
121+
fig
122+
end
123+
124+
heatmap_wigner_distribution(x_distribution, v_distribution)

‎scripts/plots/Wigner_distribution/kinetic.jl‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
using DrWatson
2+
@quickactivate "HokseonModelSimulation"
3+
using Unitful, UnitfulAtomic
4+
using CairoMakie
5+
using HokseonPlots
6+
using ColorSchemes
7+
using Colors
8+
using CSV
9+
using DataFrames
10+
using DelimitedFiles
11+
using KernelDensity, Distributions
12+
using LaTeXStrings
13+
using QuadGK, Interpolations
14+
colorscheme = ColorScheme(parse.(Colorant, ["#045275", "#089099", "#7CCBA2", "#FCDE9C", "#F0746E", "#DC3977", "#7C1D6F"]));
15+
colormap = HokseonPlots.NICECOLORS;
16+
17+
using Glob
18+
## Load the scripts in folder HokseonModelSimulation/src
19+
for file in readdir(srcdir(), join=true) # join=true returns the full path
20+
if occursin(r"\.jl$", file) # Check if the file ends with .jl
21+
include(file)
22+
end
23+
end
24+
25+
26+
is_Wigner = [true] # can be a vector of Bool
27+
28+
all_params = Dict{String, Any}(
29+
"trajectories" => [500],
30+
"nstates" => [150],
31+
"dt" => [0.05],
32+
"width" => [50],
33+
"mass" => [1.00784],
34+
"temperature" => [300.0],
35+
"tmax" => [1001],
36+
"discretisation" => [:GapGaussLegendre],
37+
"impuritymodel" => :Hokseon,
38+
"method" => [:AdiabaticIESH],
39+
"incident_energy" => [6.17],
40+
"couplings_rescale" => [2.5],
41+
"centre" => [0],
42+
"gap" => [0.49],
43+
"decoherence" => [:EDC],
44+
"is_Wigner" => is_Wigner,
45+
# 💡 Elementwise sigma assignment:
46+
"sigma" => [w ? 1.0 : nothing for w in is_Wigner],
47+
)
48+
49+
params_list = dict_list(all_params)
50+
# just make sure that params_list is a list with Dicts
51+
if typeof(params_list) != Vector{Dict{String, Any}}
52+
params_list = [params_list]
53+
end
54+
55+
56+
@unpack mass,incident_energy = params_list[1]
57+
###Initial Conditions
58+
m = ustrip(auconvert(mass*u"u"))
59+
position = austrip(5u"Å")
60+
ke = austrip(incident_energy * u"eV")
61+
velocity = - sqrt(2ke / m)
62+
@unpack is_Wigner = params_list[1]
63+
if is_Wigner
64+
@unpack sigma = params_list[1]
65+
lambda = 1/(sigma*2*m) #in a.u.
66+
# Define Gaussian distributions for position and velocity
67+
x_distribution = Normal(position, sigma) #position
68+
v_distribution = Normal(velocity, lambda) #velocity
69+
@info "Nuclei initialization: Wigner distributed"
70+
end
71+
72+
function plot_wigner_kinetic(v_distribution)
73+
# Create your figure with Minion Pro as the default font
74+
fig = Figure(
75+
size = (HokseonPlots.RESOLUTION[1] * 3, 2.5 * HokseonPlots.RESOLUTION[2]),
76+
figure_padding = (1, 20, 2, 2),
77+
fonts = (; regular = projectdir("fonts", "MinionPro-Capt.otf")),
78+
fontsize = 23
79+
)
80+
81+
82+
@unpack incident_energy = params_list[1]
83+
# Create the main axis
84+
ax = MyAxis(
85+
fig[1, 1],
86+
xlabel = "Kinetic Energy / eV", # Wrap text parts in \text{}
87+
ylabel = "Probability Density / eV⁻¹",
88+
limits = (nothing, nothing, -0.1, nothing),
89+
xgridvisible = false,
90+
ygridvisible = false
91+
)
92+
93+
94+
μ = v_distribution.μ
95+
σ = v_distribution.σ
96+
v_range_au = collect(range(μ-4σ, μ+4σ, length=1000))
97+
f_v_au = pdf.(v_distribution, v_range_au)
98+
# Convert velocity to kinetic energy
99+
KE_range_au = 0.5 * m .* v_range_au .^ 2
100+
# Transform PDF using change-of-variable: f_E(E) = f_v(v(E)) * |dv/dE|
101+
# dv/dE = 1 / (sqrt(2 * m * E))
102+
f_KE_au = f_v_au ./ (m .* sqrt.(2 .* KE_range_au ./ m))
103+
104+
# Replace original ranges with KE
105+
x_range_au = KE_range_au
106+
y_values_au = f_KE_au
107+
Hartree_to_SI = 27.211386245988
108+
x_range_SI = ustrip.(auconvert.(u"eV", collect(x_range_au)))
109+
y_values_SI = y_values_au / Hartree_to_SI
110+
lines!(ax, x_range_SI, y_values_SI, color=colormap[3], linewidth=2, label="Wigner Distribution \nsigma = $(sigma)")
111+
112+
113+
Legend(fig[1,1], ax, tellwidth=false, tellheight=false, valign=:top, halign=:right, margin=(5, 0, 5, 0), orientation=:vertical)
114+
115+
return fig
116+
end
117+
118+
119+
plot_wigner_kinetic(v_distribution)

0 commit comments

Comments
 (0)