Skip to content

Add high-l for fj channels #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
johannes-moegerle wants to merge 2 commits into from
Closed

Add high-l for fj channels #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2596851
add single_channel_fj_models
johannes-moegerle Dec 15, 2025
3656d74
add Yb173 example
johannes-moegerle Dec 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions examples/.gitignore
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
*.csv
*.parquet
./data/
70 changes: 70 additions & 0 deletions examples/generate_Yb173_database.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
using MQDT

# load Yb173 data
parameters = MQDT.Yb173.PARA
low_l_models = [
MQDT.Yb173.FMODEL_HIGHN_S15,
MQDT.Yb173.FMODEL_HIGHN_S25,
MQDT.Yb173.FMODEL_HIGHN_S35,
MQDT.Yb173.FMODEL_HIGHN_P05,
MQDT.Yb173.FMODEL_HIGHN_P15,
MQDT.Yb173.FMODEL_HIGHN_P25,
MQDT.Yb173.FMODEL_HIGHN_P35,
MQDT.Yb173.FMODEL_HIGHN_P45,
]

# bounds
n_min = 20
n_max = 30

# calculate high \nu, low \ell MQDT states
low_l_states = [eigenstates(n_min, n_max, low_l_models[i], parameters) for i in eachindex(low_l_models)]

# calculate high \ell SQDT states
l_max = n_max - 1
MQDT.wigner_init_float(n_max, "Jmax", 9) # initialize Wigner symbol caluclation
high_l_models = single_channel_fj_models(:Yb173, 5:l_max, parameters)
high_l_states = [eigenstates(25, n_max, M, parameters) for M in high_l_models]

# generate basis and calculate matrix elements
basis = basisarray(vcat(low_l_states, high_l_states), vcat(low_l_models, high_l_models))
@time me = matrix_elements(basis, parameters)

# prepare tables
using DataFrames
col_names = [:id_initial, :id_final, :value]

e1 = DataFrame(me["dipole"], col_names)
e2 = DataFrame(me["quadrupole"], col_names)
m1 = DataFrame(me["paramagnetic"], col_names)
m2 = DataFrame(me["diamagnetic"], col_names)

states_table = DataFrame(;
id=collect(1:size(basis)),
energy=MQDT.get_e(basis, parameters),
parity=MQDT.get_p(basis),
f=MQDT.get_f(basis),
nu=MQDT.get_nu(basis),
term=MQDT.get_term(basis),
lead=MQDT.get_lead(basis),
)
sort!(states_table, [:nu])

# store tables as csv files
output_dir = "data/Yb173_mqdt/"
mkpath(output_dir)

using CSV
CSV.write("$(output_dir)states.csv", states_table)
CSV.write("$(output_dir)matrix_elements_d.csv", e1)
CSV.write("$(output_dir)matrix_elements_q.csv", e2)
CSV.write("$(output_dir)matrix_elements_mu.csv", m1)
CSV.write("$(output_dir)matrix_elements_q0.csv", m2)

# store tables as parquet files
using Parquet2
Parquet2.writefile("$(output_dir)states.parquet", states_table)
Parquet2.writefile("$(output_dir)matrix_elements_d.parquet", e1)
Parquet2.writefile("$(output_dir)matrix_elements_q.parquet", e2)
Parquet2.writefile("$(output_dir)matrix_elements_mu.parquet", m1)
Parquet2.writefile("$(output_dir)matrix_elements_q0.parquet", m2)
25 changes: 14 additions & 11 deletions examples/generate_Yb174_database.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ low_l_states = [eigenstates(n_min[i], n_max, low_l_models[i], parameters) for i
# calculate high \ell SQDT states
l_max = n_max - 1
MQDT.wigner_init_float(n_max, "Jmax", 9) # initialize Wigner symbol calculation
high_l_models = single_channel_models(:Yb174, 5:l_max)
high_l_models = single_channel_jj_models(:Yb174, 5:l_max)
high_l_states = [eigenstates(25, n_max, M, parameters) for M in high_l_models]

# generate basis and calculate matrix elements
Expand Down Expand Up @@ -61,17 +61,20 @@ states_table = DataFrame(;
sort!(states_table, [:nu])

# store tables as csv files
output_dir = "data/Yb174_mqdt/"
mkpath(output_dir)

using CSV
CSV.write("Yb174_mqdt_states.csv", states_table)
CSV.write("Yb174_mqdt_matrix_elements_d.csv", e1)
CSV.write("Yb174_mqdt_matrix_elements_q.csv", e2)
CSV.write("Yb174_mqdt_matrix_elements_mu.csv", m1)
CSV.write("Yb174_mqdt_matrix_elements_q0.csv", m2)
CSV.write("$(output_dir)states.csv", states_table)
CSV.write("$(output_dir)matrix_elements_d.csv", e1)
CSV.write("$(output_dir)matrix_elements_q.csv", e2)
CSV.write("$(output_dir)matrix_elements_mu.csv", m1)
CSV.write("$(output_dir)matrix_elements_q0.csv", m2)

# store tables as parquet files
using Parquet2
Parquet2.writefile("Yb174_mqdt_states.parquet", states_table)
Parquet2.writefile("Yb174_mqdt_matrix_elements_d.parquet", e1)
Parquet2.writefile("Yb174_mqdt_matrix_elements_q.parquet", e2)
Parquet2.writefile("Yb174_mqdt_matrix_elements_mu.parquet", m1)
Parquet2.writefile("Yb174_mqdt_matrix_elements_q0.parquet", m2)
Parquet2.writefile("$(output_dir)states.parquet", states_table)
Parquet2.writefile("$(output_dir)matrix_elements_d.parquet", e1)
Parquet2.writefile("$(output_dir)matrix_elements_q.parquet", e2)
Parquet2.writefile("$(output_dir)matrix_elements_mu.parquet", m1)
Parquet2.writefile("$(output_dir)matrix_elements_q0.parquet", m2)
3 changes: 2 additions & 1 deletion src/MQDT.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,8 @@ export lsQuantumNumbers,
eigenstates,
basisarray,
matrix_elements,
single_channel_models
single_channel_jj_models,
single_channel_fj_models

include("general.jl")
include("boundstates.jl")
Expand Down
41 changes: 38 additions & 3 deletions src/general.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,7 @@ function get_thresholds(M::kModel, P::Parameters)
return thresholds
end

function single_channel_models(species::Symbol, l::Integer)
function single_channel_jj_models(species::Symbol, l::Integer)
@assert l > 0 "l must be positive and nonzero for this function"
jr = [l-1/2, l-1/2, l+1/2, l+1/2]
jt = [l-1, l, l, l+1]
Expand All @@ -667,10 +667,45 @@ function single_channel_models(species::Symbol, l::Integer)
return m
end

function single_channel_models(species::Symbol, l_list::UnitRange{Int64})
function single_channel_jj_models(species::Symbol, l_list::UnitRange{Int64})
m = Vector{fModel}()
for l in l_list
append!(m, single_channel_models(species, l))
append!(m, single_channel_jj_models(species, l))
end
return m
end

function single_channel_fj_models(species::Symbol, l_ryd::Integer, p::Parameters)
j_ryd_list = collect(abs(l_ryd - 1 / 2):1:(l_ryd + 1 / 2))
f_core_list = collect(abs(p.spin - 1 / 2):1:(p.spin + 1 / 2))
model_list = Vector{fModel}()
for j_ryd in j_ryd_list
for f_core in f_core_list
for f_tot in abs(f_core - j_ryd):1:(f_core + j_ryd)
model = fModel(
species,
"Lr=$l_ryd, Jr=$j_ryd, Fc=$f_core, F=$f_tot",
1,
[""],
Bool[1],
[0;;],
[""],
[0;;],
fjChannels([fjQuantumNumbers(0.5, 0, 0.5, f_core, l_ryd, j_ryd, f_tot)]),
fjChannels([fjQuantumNumbers(0.5, 0, 0.5, f_core, l_ryd, j_ryd, f_tot)]),
[1;;],
)
push!(model_list, model)
end
end
end
return model_list
end

function single_channel_fj_models(species::Symbol, l_ryd_list::UnitRange{Int64}, p::Parameters)
m = Vector{fModel}()
for l_ryd in l_ryd_list
append!(m, single_channel_fj_models(species, l_ryd, p))
end
return m
end
Expand Down