Get models dynamically and add high l for I!=0

.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
 

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)

utils.jl

@@ -39,25 +39,27 @@ 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 ->
             (minimum(get_relevant_lr(x[2])), minimum(get_relevant_nu(x[2])), x[1]),
     )
 
     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