More changes to fit with database-sqdt

docs/examples/angular/angular_state.ipynb

@@ -207,7 +207,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": null,
    "id": "51756a1f",
    "metadata": {},
    "outputs": [
@@ -224,7 +224,7 @@
     "ket1 = AngularKetFJ(f_c=2, l_r=1, j_r=1.5, f_tot=2.5, species=\"Yb173\")\n",
     "ket2 = AngularKetFJ(f_c=2, l_r=2, j_r=1.5, f_tot=2.5, species=\"Yb173\")\n",
     "\n",
-    "print(ket1.calc_reduced_matrix_element(ket2, operator=\"SPHERICAL\", kappa=1))\n",
+    "print(ket1.calc_reduced_matrix_element(ket2, operator=\"spherical\", kappa=1))\n",
     "print(ket1.calc_reduced_matrix_element(ket1, operator=\"s_tot\", kappa=1))"
    ]
   }

docs/examples/comparisons/compare_dipole_matrix_element.ipynb

@@ -76,7 +76,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -94,7 +94,7 @@
     "    q = round(qn2[-1] - qn1[-1])\n",
     "    state_i = RydbergStateAlkali(\"Rb\", n=qn1[0], l=qn1[1], j=qn1[2], m=qn1[3])\n",
     "    state_f = RydbergStateAlkali(\"Rb\", n=qn2[0], l=qn2[1], j=qn2[2], m=qn2[3])\n",
-    "    dipole_me = state_i.calc_matrix_element(state_f, \"ELECTRIC_DIPOLE\", q, unit=\"a.u.\")\n",
+    "    dipole_me = state_i.calc_matrix_element(state_f, \"electric_dipole\", q, unit=\"a.u.\")\n",
     "    matrixelements.append(dipole_me)\n",
     "\n",
     "results[\"ryd-numerov\"] = np.array(matrixelements)"

docs/examples/dipole_matrix_elements.ipynb

@@ -39,12 +39,12 @@
     "\n",
     "kappa = 1\n",
     "radial = state_i.radial.calc_matrix_element(state_f.radial, k_radial=1)\n",
-    "angular = state_i.angular.calc_matrix_element(state_f.angular, \"SPHERICAL\", kappa=kappa, q=0)\n",
+    "angular = state_i.angular.calc_matrix_element(state_f.angular, \"spherical\", kappa=kappa, q=0)\n",
     "prefactor = np.sqrt(4 * np.pi / (2 * kappa + 1))\n",
     "print(f\"Numerov radial matrix element: {radial}\")\n",
     "print(f\"Numerov angular matrix element: {angular}\")\n",
     "\n",
-    "dipole = state_i.calc_matrix_element(state_f, \"ELECTRIC_DIPOLE\", q=0)\n",
+    "dipole = state_i.calc_matrix_element(state_f, \"electric_dipole\", q=0)\n",
     "print(f\"Numerov dipole matrix element: {dipole}\")\n",
     "\n",
     "assert np.isclose(dipole.magnitude, (prefactor * radial * angular).magnitude), (\n",

src/ryd_numerov/__init__.py

@@ -1,10 +1,9 @@
 from ryd_numerov import angular, radial, species
-from ryd_numerov.rydberg_state import RydbergStateAlkali, RydbergStateAlkaliHyperfine, RydbergStateAlkalineLS
+from ryd_numerov.rydberg_state import RydbergStateAlkali, RydbergStateAlkalineLS
 from ryd_numerov.units import ureg
 
 __all__ = [
     "RydbergStateAlkali",
-    "RydbergStateAlkaliHyperfine",
     "RydbergStateAlkalineLS",
     "angular",
     "radial",

src/ryd_numerov/angular/angular_ket.py

@@ -17,7 +17,7 @@
     check_spin_addition_rule,
     clebsch_gordan_6j,
     clebsch_gordan_9j,
-    get_possible_quantum_number_list,
+    get_possible_quantum_number_values,
     minus_one_pow,
     try_trivial_spin_addition,
 )
@@ -97,9 +97,11 @@ def __init__(
         if species is not None:
             if isinstance(species, str):
                 species = SpeciesObject.from_name(species)
-            if i_c is not None and i_c != species.i_c:
+            # use i_c = 0 for species without defined nuclear spin (-> ignore hyperfine)
+            species_i_c = species.i_c if species.i_c is not None else 0
+            if i_c is not None and i_c != species_i_c:
                 raise ValueError(f"Nuclear spin i_c={i_c} does not match the species {species} with i_c={species.i_c}.")
-            i_c = species.i_c
+            i_c = species_i_c
             s_c = 0.5 * (species.number_valence_electrons - 1)
         if i_c is None:
             raise ValueError("Nuclear spin i_c must be set or a species must be given.")
@@ -222,11 +224,11 @@ def _to_state_ls(self) -> AngularState[AngularKetLS]:
         kets: list[AngularKetLS] = []
         coefficients: list[float] = []
 
-        s_tot_list = get_possible_quantum_number_list(self.s_c, self.s_r, getattr(self, "s_tot", None))
-        l_tot_list = get_possible_quantum_number_list(self.l_c, self.l_r, getattr(self, "l_tot", None))
+        s_tot_list = get_possible_quantum_number_values(self.s_c, self.s_r, getattr(self, "s_tot", None))
+        l_tot_list = get_possible_quantum_number_values(self.l_c, self.l_r, getattr(self, "l_tot", None))
         for s_tot in s_tot_list:
             for l_tot in l_tot_list:
-                j_tot_list = get_possible_quantum_number_list(s_tot, l_tot, getattr(self, "j_tot", None))
+                j_tot_list = get_possible_quantum_number_values(s_tot, l_tot, getattr(self, "j_tot", None))
                 for j_tot in j_tot_list:
                     try:
                         ls_ket = AngularKetLS(
@@ -257,11 +259,11 @@ def _to_state_jj(self) -> AngularState[AngularKetJJ]:
         kets: list[AngularKetJJ] = []
         coefficients: list[float] = []
 
-        j_c_list = get_possible_quantum_number_list(self.s_c, self.l_c, getattr(self, "j_c", None))
-        j_r_list = get_possible_quantum_number_list(self.s_r, self.l_r, getattr(self, "j_r", None))
+        j_c_list = get_possible_quantum_number_values(self.s_c, self.l_c, getattr(self, "j_c", None))
+        j_r_list = get_possible_quantum_number_values(self.s_r, self.l_r, getattr(self, "j_r", None))
         for j_c in j_c_list:
             for j_r in j_r_list:
-                j_tot_list = get_possible_quantum_number_list(j_c, j_r, getattr(self, "j_tot", None))
+                j_tot_list = get_possible_quantum_number_values(j_c, j_r, getattr(self, "j_tot", None))
                 for j_tot in j_tot_list:
                     try:
                         jj_ket = AngularKetJJ(
@@ -292,10 +294,10 @@ def _to_state_fj(self) -> AngularState[AngularKetFJ]:
         kets: list[AngularKetFJ] = []
         coefficients: list[float] = []
 
-        j_c_list = get_possible_quantum_number_list(self.s_c, self.l_c, getattr(self, "j_c", None))
-        j_r_list = get_possible_quantum_number_list(self.s_r, self.l_r, getattr(self, "j_r", None))
+        j_c_list = get_possible_quantum_number_values(self.s_c, self.l_c, getattr(self, "j_c", None))
+        j_r_list = get_possible_quantum_number_values(self.s_r, self.l_r, getattr(self, "j_r", None))
         for j_c in j_c_list:
-            f_c_list = get_possible_quantum_number_list(j_c, self.i_c, getattr(self, "f_c", None))
+            f_c_list = get_possible_quantum_number_values(j_c, self.i_c, getattr(self, "f_c", None))
             for f_c in f_c_list:
                 for j_r in j_r_list:
                     try:
@@ -387,7 +389,7 @@ def calc_reduced_matrix_element(  # noqa: C901
             return self.to_state().calc_reduced_matrix_element(other.to_state(), operator, kappa)
 
         qn_name: AngularMomentumQuantumNumbers
-        if operator == "SPHERICAL":
+        if operator == "spherical":
             qn_name = "l_r"
             complete_reduced_matrix_element = calc_reduced_spherical_matrix_element(self.l_r, other.l_r, kappa)
         elif operator in self.quantum_number_names:
@@ -434,9 +436,9 @@ def calc_matrix_element(self, other: AngularKetBase, operator: AngularOperatorTy
         Args:
             other: The other AngularKet :math:`|other>`.
             operator: The operator type :math:`\hat{O}_{kq}` for which to calculate the matrix element.
-                Can be one of "MAGNETIC", "ELECTRIC", "SPHERICAL".
-            kappa: The quantum number :math:`\kappa` of the angular momentum operator.
-            q: The quantum number :math:`q` of the angular momentum operator.
+                E.g. 'spherical', 's_tot', 'l_r', etc.
+            kappa: The rank :math:`\kappa` of the angular momentum operator.
+            q: The component :math:`q` of the angular momentum operator.
 
         Returns:
             The dimensionless angular matrix element.

src/ryd_numerov/angular/angular_matrix_element.py

@@ -35,7 +35,7 @@ def lru_cache(maxsize: int) -> Callable[[Callable[P, R]], Callable[P, R]]: ...
     "identity_f_tot",
 ]
 AngularOperatorType = Literal[
-    "SPHERICAL",
+    "spherical",
     AngularMomentumQuantumNumbers,
     IdentityOperators,
 ]

src/ryd_numerov/angular/utils.py

@@ -229,7 +229,7 @@ def try_trivial_spin_addition(s_1: float, s_2: float, s_tot: float | None, name:
     """
     if s_tot is None:
         if s_1 != 0 and s_2 != 0:
-            msg = f"{name} must be set if both parts ({s_1=} and {s_2=}) are non-zero."
+            msg = f"{name} must be set if both parts ({s_1} and {s_2}) are non-zero."
             raise ValueError(msg)
         s_tot = s_1 + s_2
     return float(s_tot)
@@ -245,8 +245,8 @@ def check_spin_addition_rule(s_1: float, s_2: float, s_tot: float) -> bool:
     return abs(s_1 - s_2) <= s_tot <= s_1 + s_2 and (s_1 + s_2 + s_tot) % 1 == 0
 
 
-def get_possible_quantum_number_list(s_1: float, s_2: float, s_tot: float | None) -> list[float]:
-    """Determine a list of possible s_tot from s_1 and s_2 if s_tot is not given, else return [s_tot]."""
+def get_possible_quantum_number_values(s_1: float, s_2: float, s_tot: float | None) -> list[float]:
+    """Determine a list of possible s_tot values from s_1 and s_2 if s_tot is not given, else return [s_tot]."""
     if s_tot is not None:
         return [s_tot]
     return [float(s) for s in np.arange(abs(s_1 - s_2), s_1 + s_2 + 1, 1)]

src/ryd_numerov/radial/radial_state.py

@@ -4,8 +4,6 @@
 import math
 from typing import TYPE_CHECKING, Literal, overload
 
-import numpy as np
-
 from ryd_numerov.radial.grid import Grid
 from ryd_numerov.radial.model import Model
 from ryd_numerov.radial.radial_matrix_element import calc_radial_matrix_element_from_w_z
@@ -46,12 +44,12 @@ def __init__(
         self.species = species
 
         self.n: int | None = None
-        self.nu = nu
-        self.l_r = l_r
-
-        # sanity checks
         if not nu > 0:
             raise ValueError(f"nu must be larger than 0, but is {nu=}")
+        self.nu = nu
+        if not ((isinstance(l_r, int) or l_r.is_integer()) and l_r >= 0):
+            raise ValueError(f"l_r must be an integer, and larger or equal 0, but {l_r=}")
+        self.l_r = int(l_r)
 
     def set_n_for_sanity_check(self, n: int) -> None:
         """Provide n for additional sanity checks of the radial wavefunction.
@@ -62,19 +60,17 @@ def set_n_for_sanity_check(self, n: int) -> None:
             n: Principal quantum number of the rydberg electron.
 
         """
+        if not ((isinstance(n, int) or n.is_integer()) and n >= 1):
+            raise ValueError(f"n must be an integer, and larger or equal 1, but {n=}")
+
+        if n > 10 and n < (self.nu - 1e-5):
+            # if n <= 10, we use NIST energy data for low n, which sometimes results in nu > n
+            # -1e-5: avoid issues due to numerical precision and due to NIST data
+            raise ValueError(f"n must be larger or equal to nu, but {n=}, nu={self.nu} for {self}")
+        if n <= self.l_r:
+            raise ValueError(f"n must be larger than l_r, but {n=}, l_r={self.l_r} for {self}")
         self.n = n
 
-        if self.nu > n and abs(self.nu - n) < 1e-10:
-            self.nu = n  # avoid numerical issues
-
-        if n is not None and not (isinstance(n, (int, np.integer)) and n >= 1 and n >= self.nu):
-            raise ValueError(
-                f"n must be an integer larger than 0 and larger (or equal) than nu, but is {n=}, {self.nu=}"
-            )
-
-        if not (isinstance(self.l_r, (int, np.integer)) and self.l_r >= 0 and (n is None or self.l_r <= n - 1)):
-            raise ValueError(f"l_r must be an integer, and between 0 and n - 1, but is {self.l_r=}, {n=}")
-
     def __repr__(self) -> str:
         species, nu, l_r, n = self.species, self.nu, self.l_r, self.n
         n_str = "" if n is None else f", ({n=})"