diff --git a/docs/examples/angular/angular_state.ipynb b/docs/examples/angular/angular_state.ipynb index 0125aba1..5ca6b5f2 100644 --- a/docs/examples/angular/angular_state.ipynb +++ b/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))" ] } diff --git a/docs/examples/comparisons/compare_dipole_matrix_element.ipynb b/docs/examples/comparisons/compare_dipole_matrix_element.ipynb index 14363e74..088a5e46 100644 --- a/docs/examples/comparisons/compare_dipole_matrix_element.ipynb +++ b/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)" diff --git a/docs/examples/dipole_matrix_elements.ipynb b/docs/examples/dipole_matrix_elements.ipynb index b9104059..6d0532cb 100644 --- a/docs/examples/dipole_matrix_elements.ipynb +++ b/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", diff --git a/src/ryd_numerov/__init__.py b/src/ryd_numerov/__init__.py index b26c1b1b..296f68d0 100644 --- a/src/ryd_numerov/__init__.py +++ b/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", diff --git a/src/ryd_numerov/angular/angular_ket.py b/src/ryd_numerov/angular/angular_ket.py index e068d9c7..1a749c8d 100644 --- a/src/ryd_numerov/angular/angular_ket.py +++ b/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. diff --git a/src/ryd_numerov/angular/angular_matrix_element.py b/src/ryd_numerov/angular/angular_matrix_element.py index 5b045558..b114d319 100644 --- a/src/ryd_numerov/angular/angular_matrix_element.py +++ b/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, ] diff --git a/src/ryd_numerov/angular/utils.py b/src/ryd_numerov/angular/utils.py index ea6a7950..30bca267 100644 --- a/src/ryd_numerov/angular/utils.py +++ b/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)] diff --git a/src/ryd_numerov/radial/radial_state.py b/src/ryd_numerov/radial/radial_state.py index 4b87efce..ec3191be 100644 --- a/src/ryd_numerov/radial/radial_state.py +++ b/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=})" diff --git a/src/ryd_numerov/rydberg_state.py b/src/ryd_numerov/rydberg_state.py index 456db954..4f5f07c2 100644 --- a/src/ryd_numerov/rydberg_state.py +++ b/src/ryd_numerov/rydberg_state.py @@ -4,7 +4,7 @@ import math from abc import ABC, abstractmethod from functools import cached_property -from typing import TYPE_CHECKING, get_args, overload +from typing import TYPE_CHECKING, overload import numpy as np @@ -13,23 +13,16 @@ from ryd_numerov.radial import RadialState from ryd_numerov.species.species_object import SpeciesObject from ryd_numerov.species.utils import calc_energy_from_nu -from ryd_numerov.units import BaseQuantities, MatrixElementType, ureg +from ryd_numerov.units import BaseQuantities, MatrixElementOperatorRanks, ureg if TYPE_CHECKING: from typing_extensions import Self from ryd_numerov.angular.angular_ket import AngularKetBase - from ryd_numerov.units import PintFloat + from ryd_numerov.units import MatrixElementOperator, PintFloat logger = logging.getLogger(__name__) -OPERATOR_TO_KS = { # operator: (k_radial, k_angular) - "MAGNETIC_DIPOLE": (0, 1), - "ELECTRIC_DIPOLE": (1, 1), - "ELECTRIC_QUADRUPOLE": (2, 2), - "ELECTRIC_OCTUPOLE": (3, 3), - "ELECTRIC_QUADRUPOLE_ZERO": (2, 0), -} class RydbergStateBase(ABC): @@ -70,19 +63,21 @@ def get_energy(self, unit: str | None = None) -> PintFloat | float: energy_au = calc_energy_from_nu(self.species.reduced_mass_au, nu) if unit == "a.u.": return energy_au - energy: PintFloat = energy_au * BaseQuantities["ENERGY"] + energy: PintFloat = energy_au * BaseQuantities["energy"] if unit is None: return energy return energy.to(unit, "spectroscopy").magnitude @overload - def calc_reduced_matrix_element(self, other: Self, operator: MatrixElementType, unit: None = None) -> PintFloat: ... + def calc_reduced_matrix_element( + self, other: Self, operator: MatrixElementOperator, unit: None = None + ) -> PintFloat: ... @overload - def calc_reduced_matrix_element(self, other: Self, operator: MatrixElementType, unit: str) -> float: ... + def calc_reduced_matrix_element(self, other: Self, operator: MatrixElementOperator, unit: str) -> float: ... def calc_reduced_matrix_element( - self, other: Self, operator: MatrixElementType, unit: str | None = None + self, other: Self, operator: MatrixElementOperator, unit: str | None = None ) -> PintFloat | float: r"""Calculate the reduced matrix element. @@ -105,14 +100,16 @@ def calc_reduced_matrix_element( The reduced matrix element for the given operator. """ - if operator not in get_args(MatrixElementType): - raise ValueError(f"Operator {operator} not supported, must be one of {get_args(MatrixElementType)}") + if operator not in MatrixElementOperatorRanks: + raise ValueError( + f"Operator {operator} not supported, must be one of {list(MatrixElementOperatorRanks.keys())}." + ) - k_radial, k_angular = OPERATOR_TO_KS[operator] + k_radial, k_angular = MatrixElementOperatorRanks[operator] radial_matrix_element = self.radial.calc_matrix_element(other.radial, k_radial) matrix_element: PintFloat - if operator == "MAGNETIC_DIPOLE": + if operator == "magnetic_dipole": # Magnetic dipole operator: mu = - mu_B (g_l + g_s ) g_s = 2.0023192 value_s_tot = self.angular.calc_reduced_matrix_element(other.angular, "s_tot", k_angular) @@ -125,9 +122,9 @@ def calc_reduced_matrix_element( # as the same dimensionality as the Bohr magneton (mu = - mu_B (g_l l + g_s s_tot)) # such that - mu * B (where the magnetic field B is given in dimension Tesla) is an energy - elif operator in ["ELECTRIC_DIPOLE", "ELECTRIC_QUADRUPOLE", "ELECTRIC_OCTUPOLE", "ELECTRIC_QUADRUPOLE_ZERO"]: + elif operator in ["electric_dipole", "electric_quadrupole", "electric_octupole", "electric_quadrupole_zero"]: # Electric multipole operator: p_{k,q} = e r^k_radial * sqrt(4pi / (2k+1)) * Y_{k_angular,q}(\theta, phi) - angular_matrix_element = self.angular.calc_reduced_matrix_element(other.angular, "SPHERICAL", k_angular) + angular_matrix_element = self.angular.calc_reduced_matrix_element(other.angular, "spherical", k_angular) matrix_element = ( ureg.Quantity(1, "e") * math.sqrt(4 * np.pi / (2 * k_angular + 1)) @@ -145,13 +142,13 @@ def calc_reduced_matrix_element( return matrix_element.to(unit).magnitude @overload - def calc_matrix_element(self, other: Self, operator: MatrixElementType, q: int) -> PintFloat: ... + def calc_matrix_element(self, other: Self, operator: MatrixElementOperator, q: int) -> PintFloat: ... @overload - def calc_matrix_element(self, other: Self, operator: MatrixElementType, q: int, unit: str) -> float: ... + def calc_matrix_element(self, other: Self, operator: MatrixElementOperator, q: int, unit: str) -> float: ... def calc_matrix_element( - self, other: Self, operator: MatrixElementType, q: int, unit: str | None = None + self, other: Self, operator: MatrixElementOperator, q: int, unit: str | None = None ) -> PintFloat | float: r"""Calculate the matrix element. @@ -176,7 +173,7 @@ def calc_matrix_element( The matrix element for the given operator. """ - _k_radial, k_angular = OPERATOR_TO_KS[operator] + _k_radial, k_angular = MatrixElementOperatorRanks[operator] prefactor = self.angular._calc_wigner_eckart_prefactor(other.angular, k_angular, q) # noqa: SLF001 reduced_matrix_element = self.calc_reduced_matrix_element(other, operator, unit) return prefactor * reduced_matrix_element @@ -185,61 +182,6 @@ def calc_matrix_element( class RydbergStateAlkali(RydbergStateBase): """Create an Alkali Rydberg state, including the radial and angular states.""" - def __init__( - self, - species: str | SpeciesObject, - n: int, - l: int, - j: float | None = None, - m: float | None = None, - ) -> None: - r"""Initialize the Rydberg state. - - Args: - species: Atomic species. - n: Principal quantum number of the rydberg electron. - l: Orbital angular momentum quantum number of the rydberg electron. - j: Angular momentum quantum number of the rydberg electron. - m: Total magnetic quantum number. - Optional, only needed for concrete angular matrix elements. - - """ - if isinstance(species, str): - species = SpeciesObject.from_name(species) - self.species = species - self.n = n - self.l = l - self.j = try_trivial_spin_addition(l, 0.5, j, "j") - self.m = m - - if species.number_valence_electrons != 1: - raise ValueError(f"The species {species.name} is not an alkali atom.") - if not species.is_allowed_shell(n, l): - raise ValueError(f"The shell ({n=}, {l=}) is not allowed for the species {self.species}.") - - @cached_property - def angular(self) -> AngularKetLS: - """The angular/spin state of the Rydberg electron.""" - return AngularKetLS(l_r=self.l, j_tot=self.j, m=self.m, species=self.species) - - @cached_property - def radial(self) -> RadialState: - """The radial state of the Rydberg electron.""" - radial_state = RadialState(self.species, nu=self.get_nu(), l_r=self.l) - radial_state.set_n_for_sanity_check(self.n) - return radial_state - - def __repr__(self) -> str: - species, n, l, j, m = self.species, self.n, self.l, self.j, self.m - return f"{self.__class__.__name__}({species.name}, {n=}, {l=}, {j=}, {m=})" - - def get_nu(self) -> float: - return self.species.calc_nu(self.n, self.l, self.j, s_tot=1 / 2) - - -class RydbergStateAlkaliHyperfine(RydbergStateBase): - """Create an Alkali Rydberg state, including the radial and angular states.""" - def __init__( self, species: str | SpeciesObject, @@ -257,6 +199,7 @@ def __init__( l: Orbital angular momentum quantum number of the rydberg electron. j: Angular momentum quantum number of the rydberg electron. f: Total angular momentum quantum number of the atom (rydberg electron + core) + Optional, only needed if the species supports hyperfine structure (i.e. species.i_c is not None or 0). m: Total magnetic quantum number. Optional, only needed for concrete angular matrix elements. @@ -264,10 +207,11 @@ def __init__( if isinstance(species, str): species = SpeciesObject.from_name(species) self.species = species + i_c = species.i_c if species.i_c is not None else 0 self.n = n self.l = l self.j = try_trivial_spin_addition(l, 0.5, j, "j") - self.f = try_trivial_spin_addition(self.j, species.i_c, f, "f") + self.f = try_trivial_spin_addition(self.j, i_c, f, "f") self.m = m if species.number_valence_electrons != 1: @@ -303,8 +247,9 @@ def __init__( species: str | SpeciesObject, n: int, l: int, - s_tot: float, - j_tot: float | None = None, + s_tot: int, + j_tot: int | None = None, + f_tot: float | None = None, m: float | None = None, ) -> None: r"""Initialize the Rydberg state. @@ -315,6 +260,8 @@ def __init__( l: Orbital angular momentum quantum number of the rydberg electron. s_tot: Total spin quantum number of all electrons. j_tot: Total angular momentum quantum number of all electrons. + f_tot: Total angular momentum quantum number of the atom (rydberg electron + core) + Optional, only needed if the species supports hyperfine structure (i.e. species.i_c is not None or 0). m: Total magnetic quantum number. Optional, only needed for concrete angular matrix elements. @@ -322,10 +269,12 @@ def __init__( if isinstance(species, str): species = SpeciesObject.from_name(species) self.species = species + i_c = species.i_c if species.i_c is not None else 0 self.n = n self.l = l self.s_tot = s_tot self.j_tot = try_trivial_spin_addition(l, s_tot, j_tot, "j_tot") + self.f_tot = try_trivial_spin_addition(self.j_tot, i_c, f_tot, "f_tot") self.m = m if species.number_valence_electrons != 2: @@ -336,7 +285,9 @@ def __init__( @cached_property def angular(self) -> AngularKetLS: """The angular/spin state of the Rydberg electron.""" - return AngularKetLS(l_r=self.l, s_tot=self.s_tot, j_tot=self.j_tot, m=self.m, species=self.species) + return AngularKetLS( + l_r=self.l, s_tot=self.s_tot, j_tot=self.j_tot, f_tot=self.f_tot, m=self.m, species=self.species + ) @cached_property def radial(self) -> RadialState: diff --git a/src/ryd_numerov/species/species_object.py b/src/ryd_numerov/species/species_object.py index e9a6ebda..2ae1c1ef 100644 --- a/src/ryd_numerov/species/species_object.py +++ b/src/ryd_numerov/species/species_object.py @@ -34,8 +34,8 @@ class SpeciesObject(ABC): """The name of the atomic species.""" Z: ClassVar[int] """Atomic number of the species.""" - i_c: ClassVar[float] = 0 - """Nuclear spin, (default 0 to ignore hyperfine structure).""" + i_c: ClassVar[float | None] = None + """Nuclear spin, (default None to ignore hyperfine structure, will be treated like i_c = 0).""" number_valence_electrons: ClassVar[int] """Number of valence electrons (i.e. 1 for alkali atoms and 2 for alkaline earth atoms).""" ground_state_shell: ClassVar[tuple[int, int]] diff --git a/src/ryd_numerov/units.py b/src/ryd_numerov/units.py index 7572dbcd..77488655 100644 --- a/src/ryd_numerov/units.py +++ b/src/ryd_numerov/units.py @@ -17,54 +17,58 @@ ureg: UnitRegistry[float] = UnitRegistry(system="atomic") -MatrixElementType = Literal[ - "MAGNETIC_DIPOLE", # MAGNETIC with k_radial = 0, k_angular = 1 - "ELECTRIC_DIPOLE", # ELECTRIC with k_radial = 1, k_angular = 1 - "ELECTRIC_QUADRUPOLE", # ELECTRIC with k_radial = 2, k_angular = 2 - "ELECTRIC_OCTUPOLE", # ELECTRIC with k_radial = 3, k_angular = 3 - "ELECTRIC_QUADRUPOLE_ZERO", # ELECTRIC with k_radial = 2, k_angular = 0 +MatrixElementOperator = Literal[ + "magnetic_dipole", + "electric_dipole", + "electric_quadrupole", + "electric_octupole", + "electric_quadrupole_zero", ] +MatrixElementOperatorRanks: dict[MatrixElementOperator, tuple[int, int]] = { + # "operator": (k_radial, k_angular) + "magnetic_dipole": (0, 1), + "electric_dipole": (1, 1), + "electric_quadrupole": (2, 2), + "electric_octupole": (3, 3), + "electric_quadrupole_zero": (2, 0), +} Dimension = Literal[ - "ELECTRIC_FIELD", - "MAGNETIC_FIELD", - "DISTANCE", - "ENERGY", - "CHARGE", - "VELOCITY", - "TEMPERATURE", - "TIME", - "RADIAL_MATRIX_ELEMENT", - "ANGULAR_MATRIX_ELEMENT", - "ELECTRIC_DIPOLE", - "ELECTRIC_QUADRUPOLE", - "ELECTRIC_QUADRUPOLE_ZERO", - "ELECTRIC_OCTUPOLE", - "MAGNETIC_DIPOLE", - "ARBITRARY", - "ZERO", + MatrixElementOperator, + "electric_field", + "magnetic_field", + "distance", + "energy", + "charge", + "velocity", + "temperature", + "time", + "radial_matrix_element", + "angular_matrix_element", + "arbitrary", + "zero", ] DimensionLike = Union[Dimension, tuple[Dimension, Dimension]] # some abbreviations: au_time: atomic_unit_of_time; au_current: atomic_unit_of_current; m_e: electron_mass _CommonUnits: dict[Dimension, str] = { - "ELECTRIC_FIELD": "V/cm", # 1 V/cm = 1.9446903811524456e-10 bohr * m_e / au_current / au_time ** 3 - "MAGNETIC_FIELD": "T", # 1 T = 4.254382157342044e-06 m_e / au_current / au_time ** 2 - "DISTANCE": "micrometer", # 1 mum = 18897.26124622279 bohr - "ENERGY": "hartree", # 1 hartree = 1 bohr ** 2 * m_e / au_time ** 2 - "CHARGE": "e", # 1 e = 1 au_current * au_time - "VELOCITY": "speed_of_light", # 1 c = 137.03599908356244 bohr / au_time - "TEMPERATURE": "K", # 1 K = 3.1668115634555572e-06 atomic_unit_of_temperature - "TIME": "s", # 1 s = 4.134137333518244e+16 au_time - "RADIAL_MATRIX_ELEMENT": "bohr", # 1 bohr - "ANGULAR_MATRIX_ELEMENT": "", # 1 dimensionless - "ELECTRIC_DIPOLE": "e * a0", # 1 e * a0 = 1 au_current * au_time * bohr - "ELECTRIC_QUADRUPOLE": "e * a0^2", # 1 e * a0^2 = 1 au_current * au_time * bohr ** 2 - "ELECTRIC_QUADRUPOLE_ZERO": "e * a0^2", # 1 e * a0^2 = 1 au_current * au_time * bohr ** 2 - "ELECTRIC_OCTUPOLE": "e * a0^3", # 1 e * a0^3 = 1 au_current * au_time * bohr ** 3 - "MAGNETIC_DIPOLE": "bohr_magneton", # 1 bohr_magneton = 0.5 au_current * bohr ** 2' - "ARBITRARY": "", # 1 dimensionless - "ZERO": "", # 1 dimensionless + "electric_field": "V/cm", # 1 V/cm = 1.9446903811524456e-10 bohr * m_e / au_current / au_time ** 3 + "magnetic_field": "T", # 1 T = 4.254382157342044e-06 m_e / au_current / au_time ** 2 + "distance": "micrometer", # 1 mum = 18897.26124622279 bohr + "energy": "hartree", # 1 hartree = 1 bohr ** 2 * m_e / au_time ** 2 + "charge": "e", # 1 e = 1 au_current * au_time + "velocity": "speed_of_light", # 1 c = 137.03599908356244 bohr / au_time + "temperature": "K", # 1 K = 3.1668115634555572e-06 atomic_unit_of_temperature + "time": "s", # 1 s = 4.134137333518244e+16 au_time + "radial_matrix_element": "bohr", # 1 bohr + "angular_matrix_element": "", # 1 dimensionless + "electric_dipole": "e * a0", # 1 e * a0 = 1 au_current * au_time * bohr + "electric_quadrupole": "e * a0^2", # 1 e * a0^2 = 1 au_current * au_time * bohr ** 2 + "electric_quadrupole_zero": "e * a0^2", # 1 e * a0^2 = 1 au_current * au_time * bohr ** 2 + "electric_octupole": "e * a0^3", # 1 e * a0^3 = 1 au_current * au_time * bohr ** 3 + "magnetic_dipole": "bohr_magneton", # 1 bohr_magneton = 0.5 au_current * bohr ** 2' + "arbitrary": "", # 1 dimensionless + "zero": "", # 1 dimensionless } BaseUnits: dict[Dimension, PlainUnit] = { k: ureg.Quantity(1, unit).to_base_units().units for k, unit in _CommonUnits.items() @@ -73,8 +77,8 @@ Context = Literal["spectroscopy", "Gaussian"] BaseContexts: dict[Dimension, Context] = { - "MAGNETIC_FIELD": "Gaussian", - "ENERGY": "spectroscopy", + "magnetic_field": "Gaussian", + "energy": "spectroscopy", } diff --git a/tests/test_all_elements.py b/tests/test_all_elements.py index 05472b0d..dedf25e5 100644 --- a/tests/test_all_elements.py +++ b/tests/test_all_elements.py @@ -15,7 +15,10 @@ def test_magnetic(species_name: str) -> None: state: RydbergStateBase if species.number_valence_electrons == 1: - state = RydbergStateAlkali(species, n=50, l=0) + if species.i_c is None: + state = RydbergStateAlkali(species, n=50, l=0) + else: + state = RydbergStateAlkali(species, n=50, l=0, f=species.i_c + 0.5) state.radial.create_wavefunction() with pytest.raises(ValueError, match="j must be set"): RydbergStateAlkali(species, n=50, l=1) diff --git a/tests/test_angular_matrix_elements.py b/tests/test_angular_matrix_elements.py index 5c04279a..1f03ca26 100644 --- a/tests/test_angular_matrix_elements.py +++ b/tests/test_angular_matrix_elements.py @@ -79,10 +79,10 @@ def test_reduced_identity(ket: AngularKetBase) -> None: @pytest.mark.parametrize(("ket1", "ket2"), TEST_KET_PAIRS) def test_matrix_elements_in_different_coupling_schemes(ket1: AngularKetBase, ket2: AngularKetBase) -> None: example_list: list[tuple[AngularOperatorType, int]] = [ - ("SPHERICAL", 0), - ("SPHERICAL", 1), - ("SPHERICAL", 2), - ("SPHERICAL", 3), + ("spherical", 0), + ("spherical", 1), + ("spherical", 2), + ("spherical", 3), ("s_tot", 1), ("l_r", 1), ("i_c", 1), diff --git a/tests/test_matrix_elements.py b/tests/test_matrix_elements.py index 88313f74..a923abcc 100644 --- a/tests/test_matrix_elements.py +++ b/tests/test_matrix_elements.py @@ -13,13 +13,13 @@ def test_magnetic(l: int) -> None: state = RydbergStateAlkali("Rb", n=max(l + 1, 10), l=l, j=l + 0.5, m=l + 0.5) # Check that for m = j_tot = l + s_tot the magnetic matrix element is - mu_B * (g_l * l + g_s * s_tot) - mu = state.calc_matrix_element(state, "MAGNETIC_DIPOLE", q=0) + mu = state.calc_matrix_element(state, "magnetic_dipole", q=0) mu = mu.to("bohr_magneton") assert np.isclose(mu.magnitude, -(g_l * l + g_s * 0.5)), f"{mu.magnitude} != {-(g_l * l + g_s * 0.5)}" # Check dimensionality magnetic_field = ureg.Quantity(1, "T") zeeman_energy = -mu * magnetic_field - assert zeeman_energy.dimensionality == BaseUnits["ENERGY"].dimensionality, ( - f"{zeeman_energy.dimensionality} != {BaseUnits['ENERGY'].dimensionality}" + assert zeeman_energy.dimensionality == BaseUnits["energy"].dimensionality, ( + f"{zeeman_energy.dimensionality} != {BaseUnits['energy'].dimensionality}" )