diff --git a/.github/workflows/generate_database.yml b/.github/workflows/generate_database.yml index d8b19d7..648591e 100644 --- a/.github/workflows/generate_database.yml +++ b/.github/workflows/generate_database.yml @@ -14,16 +14,15 @@ jobs: matrix: species: [Sr87, Sr88, Yb171, Yb173, Yb174] include: - # use n_max=220 for tags and master, but only n_max=100 for PRs, etc. - - n_max: ${{ (startsWith(github.ref, 'refs/tags/') || contains(github.ref, 'master')) && 220 || 100 }} - # FIXME: the following species do not yet support high \ell states + # use n_max=220 for tags and main, but only n_max=110 for PRs, etc. + - n_max: ${{ (startsWith(github.ref, 'refs/tags/') || contains(github.ref, 'main')) && 220 || 110 }} + # for isotopes with hyperfine structure, we only include high-l states up to n=110 (otherwise the database size gets too large) - species: Sr87 - extra_args: "--skip-high-l" + extra_args: "--n-max-high-l 110" - species: Yb171 - extra_args: "--skip-high-l" + extra_args: "--n-max-high-l 110" - species: Yb173 - extra_args: "--skip-high-l" # FIXME: Yb173 fails for n > 115 - n_max: ${{ (startsWith(github.ref, 'refs/tags/') || contains(github.ref, 'master')) && 115 || 100 }} + extra_args: "--n-max-high-l 110" name: ${{ matrix.species }} runs-on: ubuntu-latest timeout-minutes: 300 @@ -37,7 +36,7 @@ jobs: - uses: julia-actions/setup-julia@v2 with: - version: 1.11 + version: 1.12 arch: x64 - uses: julia-actions/cache@v2 diff --git a/MQDT.jl b/MQDT.jl index df4e47d..9cb11fd 160000 --- a/MQDT.jl +++ b/MQDT.jl @@ -1 +1 @@ -Subproject commit df4e47dc2f4df56f0a42260e921ce56389b36657 +Subproject commit 9cb11fd5f2a930c8ee13ae388d9a7acbc228c43a diff --git a/generate_database.jl b/generate_database.jl index a497187..d38fd1e 100644 --- a/generate_database.jl +++ b/generate_database.jl @@ -10,8 +10,9 @@ using LRUCache include("utils.jl") -include("tables.jl") -version = "v1.1" +HIGH_L = 10 + +version = "v1.2" function parse_commandline() s = ArgParseSettings( @@ -24,13 +25,17 @@ function parse_commandline() help = "The species to generate the database for" required = true arg_type = String - "--n-min-high-l" - help = "The minimal principal quantum number n for the high angular momentum (l) states to be included in the database." + "--n-min-sqdt" + help = "The minimal principal quantum number n for the SQDT models to be included in the database." default = 25 arg_type = Int + "--n-max-high-l" + help = "The maximal principal quantum number n for high angular momentum states (l > $HIGH_L) to be included in the database." + default = Inf + arg_type = Float64 "--n-max" help = "The maximum principal quantum number n for the states to be included in the database." - default = 90 + default = 110 arg_type = Int "--directory" help = "The directory where the database will be saved" @@ -39,9 +44,6 @@ function parse_commandline() "--overwrite" help = "Delete the species folder if it exists and create a new one" action = :store_true - "--skip-high-l" - help = "Include high angular momentum (l) states in the calculation" - action = :store_true end return parse_args(s) @@ -51,8 +53,9 @@ end function main() # Parse command line arguments args = parse_commandline() - n_min_high_l = args["n-min-high-l"] n_max = args["n-max"] + n_min_sqdt = args["n-min-sqdt"] + n_max_high_l = args["n-max-high-l"] directory = args["directory"] overwrite = args["overwrite"] species = Symbol(args["species"]) @@ -73,49 +76,72 @@ function main() global_logger(combined_logger) @info "Starting database generation for $species with version $version" - @info "Parameters: n_min_high_l=$n_min_high_l, n_max=$n_max" + @info "Parameters: $args" start_time = time() # initialize Wigner symbol calculation - if args["skip-high-l"] - CGcoefficient.wigner_init_float(max(FMODEL_MAX_L[species], 5), "Jmax", 9) - else - CGcoefficient.wigner_init_float(n_max - 1, "Jmax", 9) - end - parameters = PARA_TABLE[species] - models = MODELS_TABLE[species] - - if args["skip-high-l"] - @info "Skipping high ℓ SQDT models." - high_l_models = [] - else - @info "Generating high ℓ SQDT models..." - l_max = n_max - 1 - l_start = FMODEL_MAX_L[species] + 1 - high_l_models = MQDT.single_channel_models(species, l_start:l_max) - models = vcat(models, high_l_models) + CGcoefficient.wigner_init_float(n_max - 1, "Jmax", 9) + + parameters = MQDT.get_species_parameters(species) + all_models = Vector{MQDT.fModel}() + + s_r = 1 / 2 + j_c = 1 / 2 + i_c = parameters.spin + seen_models = Set{String}() + for l_r = 0:(n_max-1) + for j_r = abs(l_r-s_r):1:(l_r+s_r) + for f_c = abs(j_c-i_c):1:(j_c+i_c) + for f_tot = abs(f_c-j_r):1:(f_c+j_r) + models = MQDT.get_fmodels(species, l_r, j_r, f_c, f_tot, parameters) + for model in models + if !(model.name in seen_models) + push!(seen_models, model.name) + push!(all_models, model) + end + end + end + end + end end @info "Calculating MQDT states..." - states = Vector{MQDT.EigenStates}(undef, length(models)) - for (i, M) in enumerate(models) - n_min = M in high_l_models ? n_min_high_l : NaN + states = Vector{MQDT.EigenStates}(undef, length(all_models)) + for (i, M) in enumerate(all_models) + _n_min = NaN + _n_max = n_max + if startswith(M.name, "SQDT") + _n_min = n_min_sqdt + if parameters.spin > 0 + l_ryd = MQDT.get_lr(M)[1] + if l_ryd >= n_max_high_l + _n_max = 0 + elseif l_ryd > HIGH_L + _n_max = min(n_max_high_l, n_max) + end + end + end @info "$(M.name)" - states[i] = MQDT.eigenstates(n_min, n_max, M, parameters) - @info " found nu_min=$(minimum(states[i].n)), nu_max=$(maximum(states[i].n)), total states=$(length(states[i].n))" + states[i] = MQDT.eigenstates(_n_min, _n_max, M, parameters) + if length(states[i].n) > 0 + @info " found nu_min=$(minimum(states[i].n)), nu_max=$(maximum(states[i].n)), total states=$(length(states[i].n))" + else + @info " found no states" + end end - basis = MQDT.basisarray(states, models) + basis = MQDT.basisarray(states, all_models) @info "Generated state table with $(length(basis.states)) states" @info "Converting states to database table..." - states_df = basis_to_df(basis, parameters) + @timelog states_df = basis_to_df(basis, parameters) @info "Calculating matrix elements..." @timelog row_col_value_dict = all_matrix_element(basis, parameters) @info "Converting matrix elements to database table..." - matrix_elements_df_dict = Dict(k => rcv_to_df(v) for (k, v) in row_col_value_dict) + @timelog matrix_elements_df_dict = + Dict(k => rcv_to_df(v) for (k, v) in row_col_value_dict) @info "Storing database tables as parquet files..." tables = merge(Dict("states" => states_df), matrix_elements_df_dict) diff --git a/rydstate b/rydstate index b2d5597..2095748 160000 --- a/rydstate +++ b/rydstate @@ -1 +1 @@ -Subproject commit b2d559760779cfb2ff0158b9a05a00b9d3b53937 +Subproject commit 2095748c6faa61a9221f375a24d52a14901c9630 diff --git a/tables.jl b/tables.jl deleted file mode 100644 index a1c8eec..0000000 --- a/tables.jl +++ /dev/null @@ -1,92 +0,0 @@ -using MQDT - -const MODELS_TABLE = Dict( - :Sr87 => [ - MQDT.Sr87.FMODEL_LOWN_P45, - MQDT.Sr87.FMODEL_HIGHN_S35, - MQDT.Sr87.FMODEL_HIGHN_S45, - MQDT.Sr87.FMODEL_HIGHN_S55, - MQDT.Sr87.FMODEL_HIGHN_P25, - MQDT.Sr87.FMODEL_HIGHN_P35, - MQDT.Sr87.FMODEL_HIGHN_P45, - MQDT.Sr87.FMODEL_HIGHN_P55, - MQDT.Sr87.FMODEL_HIGHN_P65, - MQDT.Sr87.FMODEL_HIGHN_D15, - MQDT.Sr87.FMODEL_HIGHN_D25, - MQDT.Sr87.FMODEL_HIGHN_D35, - MQDT.Sr87.FMODEL_HIGHN_D45, - MQDT.Sr87.FMODEL_HIGHN_D55, - MQDT.Sr87.FMODEL_HIGHN_D65, - MQDT.Sr87.FMODEL_HIGHN_D75, - MQDT.Sr87.FMODEL_HIGHN_F45, - ], - :Sr88 => [ - MQDT.Sr88.FMODEL_LOWN_P1, - MQDT.Sr88.FMODEL_HIGHN_S0, - MQDT.Sr88.FMODEL_HIGHN_S1, - MQDT.Sr88.FMODEL_HIGHN_P0, - MQDT.Sr88.FMODEL_HIGHN_P1, - MQDT.Sr88.FMODEL_HIGHN_P2, - MQDT.Sr88.FMODEL_HIGHN_D1, - MQDT.Sr88.FMODEL_HIGHN_D2, - MQDT.Sr88.FMODEL_HIGHN_D3, - MQDT.Sr88.FMODEL_HIGHN_F2, - MQDT.Sr88.FMODEL_HIGHN_F3, - MQDT.Sr88.FMODEL_HIGHN_F4, - ], - :Yb171 => [ - MQDT.Yb171.FMODEL_HIGHN_S05, - MQDT.Yb171.FMODEL_HIGHN_S15, - MQDT.Yb171.FMODEL_HIGHN_P05, - MQDT.Yb171.FMODEL_HIGHN_P15, - MQDT.Yb171.FMODEL_HIGHN_D05, - MQDT.Yb171.FMODEL_HIGHN_D15, - MQDT.Yb171.FMODEL_HIGHN_D25, - MQDT.Yb171.FMODEL_HIGHN_D35, - MQDT.Yb171.FMODEL_HIGHN_F25, - MQDT.Yb171.FMODEL_HIGHN_F35, - MQDT.Yb171.FMODEL_HIGHN_F45, - MQDT.Yb171.FMODEL_HIGHN_G25, - MQDT.Yb171.FMODEL_HIGHN_G35, - MQDT.Yb171.FMODEL_HIGHN_G45, - MQDT.Yb171.FMODEL_HIGHN_G55, - ], - :Yb173 => [ - 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, - ], - :Yb174 => [ - MQDT.Yb174.FMODEL_LOWN_P1, - MQDT.Yb174.FMODEL_HIGHN_S0, - MQDT.Yb174.FMODEL_HIGHN_S1, - MQDT.Yb174.FMODEL_HIGHN_P0, - MQDT.Yb174.FMODEL_HIGHN_P1, - MQDT.Yb174.FMODEL_HIGHN_P2, - MQDT.Yb174.FMODEL_HIGHN_D1, - MQDT.Yb174.FMODEL_HIGHN_D2, - MQDT.Yb174.FMODEL_HIGHN_D3, - MQDT.Yb174.FMODEL_HIGHN_F2, - MQDT.Yb174.FMODEL_HIGHN_F3, - MQDT.Yb174.FMODEL_HIGHN_F4, - MQDT.Yb174.FMODEL_HIGHN_G3, - MQDT.Yb174.FMODEL_HIGHN_G4, - MQDT.Yb174.FMODEL_HIGHN_G5, - ], -) - -const FMODEL_MAX_L = Dict(:Sr87 => 3, :Sr88 => 3, :Yb171 => 4, :Yb173 => 1, :Yb174 => 4) - - -const PARA_TABLE = Dict( - :Sr87 => MQDT.Sr87.PARA, - :Sr88 => MQDT.Sr88.PARA, - :Yb171 => MQDT.Yb171.PARA, - :Yb173 => MQDT.Yb173.PARA, - :Yb174 => MQDT.Yb174.PARA, -) diff --git a/utils.jl b/utils.jl index 50145f8..d98848d 100644 --- a/utils.jl +++ b/utils.jl @@ -39,7 +39,10 @@ function all_matrix_element(B::MQDT.BasisArray, parameters::MQDT.Parameters) "matrix_elements_q0" => Tuple{Int64,Int64,Float64}[], ) - states_indexed = [(ids - 1 + START_ID, state) for (ids, state) in enumerate(B.states)] + states_indexed = [(ids, state) for (ids, state) in enumerate(B.states)] + states_lr = [get_relevant_lr(s) for (_, s) in states_indexed] + states_nu = [get_relevant_nu(s) for (_, s) in states_indexed] + states_sorted = sort( states_indexed, by = x -> @@ -47,17 +50,16 @@ function all_matrix_element(B::MQDT.BasisArray, parameters::MQDT.Parameters) ) for (i1, (id1, s1)) in enumerate(states_sorted) - lr1 = get_relevant_lr(s1) - nus1 = get_relevant_nu(s1) + lr1 = states_lr[id1] + nus1 = states_nu[id1] for (id2, s2) in states_sorted[i1:end] - lr2 = get_relevant_lr(s2) - nus2 = get_relevant_nu(s2) - + lr2 = states_lr[id2] # Skip if all contributions of the two states are far apart in angular momentum if minimum(lr2) - maximum(lr1) > k_angular_max continue end + nus2 = states_nu[id2] # Skip if all contributions of the two states are far apart in n and None of them is low-n if all(abs(nu1-nu2) >= 11 for nu1 in nus1 for nu2 in nus2) && all(nu1 > 25 for nu1 in nus1) && @@ -78,9 +80,12 @@ function all_matrix_element(B::MQDT.BasisArray, parameters::MQDT.Parameters) for (i, key) in enumerate(table_keys) if m[i] != 0 - push!(row_col_value[key], (id1, id2, m[i])) + # start IDs from 0 for consistency with python + _id1 = id1 - 1 + START_ID + _id2 = id2 - 1 + START_ID + push!(row_col_value[key], (_id1, _id2, m[i])) if id1 != id2 - push!(row_col_value[key], (id2, id1, m[i] * prefactor_transposed)) + push!(row_col_value[key], (_id2, _id1, m[i] * prefactor_transposed)) end end end