Line data Source code
1 : # SPDX-FileCopyrightText: 2025 Pairinteraction Developers
2 : # SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : """Test receiving matrix elements from the databases."""
5 :
6 1 : from collections.abc import Generator
7 1 : from pathlib import Path
8 :
9 1 : import duckdb
10 1 : import numpy as np
11 1 : import pytest
12 1 : from packaging.version import Version
13 1 : from sympy.physics.wigner import wigner_3j
14 :
15 1 : import pairinteraction.real as pi
16 1 : from tests.constants import GAUSS_IN_ATOMIC_UNITS, HARTREE_IN_GHZ, SPECIES_TO_NUCLEAR_SPIN, SUPPORTED_SPECIES
17 :
18 1 : expected_operator = np.array([[3.58588117, 1.66420213], [1.66420213, 4.16645123]])
19 :
20 :
21 1 : def fetch_id(n: int, l: float, f: float, s: float, connection: duckdb.DuckDBPyConnection, table: str) -> int:
22 1 : return connection.execute(
23 : f"SELECT id FROM '{table}' WHERE n = {n} AND f = {f} ORDER BY (exp_s - {s})^2+(exp_l - {l})^2 LIMIT 1"
24 : ).fetchone()[0]
25 :
26 :
27 1 : def fetch_wigner_element(
28 : f_initial: float,
29 : f_final: float,
30 : m_initial: float,
31 : m_final: float,
32 : kappa: int,
33 : q: int,
34 : connection: duckdb.DuckDBPyConnection,
35 : table: str,
36 : ) -> float:
37 1 : result = connection.execute(
38 : f"SELECT val FROM '{table}' WHERE f_initial = {f_initial} AND f_final = {f_final} AND "
39 : f"m_initial = {m_initial} AND m_final = {m_final} AND kappa = {kappa} AND q = {q}"
40 : ).fetchone()
41 1 : return result[0] if result else 0
42 :
43 :
44 1 : def fetch_reduced_matrix_element(
45 : id_initial: int, id_final: int, connection: duckdb.DuckDBPyConnection, table: str
46 : ) -> float:
47 1 : result = connection.execute(
48 : f"SELECT val FROM '{table}' WHERE id_initial = {id_initial} AND id_final = {id_final}"
49 : ).fetchone()
50 1 : return result[0] if result else 0
51 :
52 :
53 1 : @pytest.fixture(scope="module")
54 1 : def connection() -> Generator[duckdb.DuckDBPyConnection]:
55 1 : with duckdb.connect(":memory:") as connection:
56 1 : yield connection
57 :
58 :
59 1 : @pytest.mark.parametrize("swap_states", [False, True])
60 1 : def test_database(connection: duckdb.DuckDBPyConnection, swap_states: bool) -> None:
61 : """Test receiving matrix elements from the databases."""
62 1 : database = pi.Database.get_global_database()
63 1 : bfield_in_gauss = 1500
64 :
65 : # Define initial and final quantum states
66 1 : n_initial, n_final = 54, 54
67 1 : l_initial, l_final = 1, 1
68 1 : f_initial, f_final = 1, 0
69 1 : m_initial, m_final = 0, 0
70 1 : s_initial, s_final = 0.6, 1.0
71 :
72 : # Swap states if required by the test parameter
73 1 : if swap_states:
74 1 : n_initial, n_final = n_final, n_initial
75 1 : l_initial, l_final = l_final, l_initial
76 1 : f_initial, f_final = f_final, f_initial
77 1 : m_initial, m_final = m_final, m_initial
78 1 : s_initial, s_final = s_final, s_initial
79 :
80 : # Get the Zeeman interaction operator from the database using pairinteraction
81 1 : ket_initial = pi.KetAtom("Yb174_mqdt", n=n_initial, l=l_initial, f=f_initial, m=m_initial, s=s_initial)
82 1 : ket_final = pi.KetAtom("Yb174_mqdt", n=n_final, l=l_final, f=f_final, m=m_final, s=s_final)
83 1 : basis = pi.BasisAtom("Yb174_mqdt", additional_kets=[ket_initial, ket_final])
84 1 : operator = (
85 : pi.SystemAtom(basis)
86 : .set_magnetic_field([0, 0, bfield_in_gauss], unit="G")
87 : .set_diamagnetism_enabled(True)
88 : .get_hamiltonian(unit="GHz")
89 : )
90 1 : operator -= np.diag(np.sort([ket_initial.get_energy(unit="GHz"), ket_final.get_energy(unit="GHz")]))
91 1 : assert np.allclose(operator, expected_operator)
92 :
93 : # Get the latest parquet files from the database directory
94 1 : parquet_files = {}
95 1 : parquet_versions = {}
96 1 : for path in list(Path(database.database_dir).rglob("*.parquet")):
97 1 : species, version_str = path.parent.name.rsplit("_v", 1)
98 1 : table = path.stem
99 1 : name = f"{species}_{table}"
100 1 : version = Version(version_str)
101 1 : if name not in parquet_files or version > parquet_versions[name]:
102 1 : parquet_files[name] = path
103 1 : parquet_versions[name] = version
104 1 : assert "misc_wigner" in parquet_files
105 1 : assert "Yb174_mqdt_states" in parquet_files
106 1 : assert "Yb174_mqdt_matrix_elements_mu" in parquet_files
107 1 : assert "Yb174_mqdt_matrix_elements_q" in parquet_files
108 1 : assert "Yb174_mqdt_matrix_elements_q0" in parquet_files
109 :
110 : # Obtain the ids of the initial and final states
111 1 : id_initial = fetch_id(n_initial, l_initial, f_initial, s_initial, connection, parquet_files["Yb174_mqdt_states"])
112 1 : assert id_initial == 362 if swap_states else 363
113 :
114 1 : id_final = fetch_id(n_final, l_final, f_final, s_final, connection, parquet_files["Yb174_mqdt_states"])
115 1 : assert id_final == 363 if swap_states else 362
116 :
117 : # Obtain a matrix element of the magnetic dipole operator (for the chosen kets, it is non-zero iff initial != final)
118 1 : kappa, q = 1, 0
119 1 : wigner_element = fetch_wigner_element(
120 : f_initial, f_final, m_initial, m_final, kappa, q, connection, parquet_files["misc_wigner"]
121 : )
122 1 : assert np.isclose(
123 : wigner_element,
124 : float((-1) ** (f_final - m_final) * wigner_3j(f_final, kappa, f_initial, -m_final, q, m_initial)),
125 : )
126 :
127 1 : mu = fetch_reduced_matrix_element(id_initial, id_final, connection, parquet_files["Yb174_mqdt_matrix_elements_mu"])
128 1 : matrix_element = -wigner_element * mu * bfield_in_gauss * GAUSS_IN_ATOMIC_UNITS * HARTREE_IN_GHZ
129 1 : assert np.isclose(matrix_element, operator[0, 1])
130 :
131 : # Obtain a matrix element of the diamagnetic operator (for the chosen kets, it is non-zero iff initial == final)
132 1 : n_final = n_initial
133 1 : l_final = l_initial
134 1 : f_final = f_initial
135 1 : m_final = m_initial
136 1 : s_final = s_initial
137 1 : id_final = id_initial
138 :
139 1 : kappa, q = 0, 0
140 1 : wigner_element = fetch_wigner_element(
141 : f_initial, f_final, m_initial, m_final, kappa, q, connection, parquet_files["misc_wigner"]
142 : )
143 1 : assert np.isclose(
144 : wigner_element,
145 : float((-1) ** (f_final - m_final) * wigner_3j(f_final, kappa, f_initial, -m_final, q, m_initial)),
146 : )
147 :
148 1 : q0 = fetch_reduced_matrix_element(id_initial, id_final, connection, parquet_files["Yb174_mqdt_matrix_elements_q0"])
149 1 : matrix_element = 1 / 12 * wigner_element * q0 * (bfield_in_gauss * GAUSS_IN_ATOMIC_UNITS) ** 2 * HARTREE_IN_GHZ
150 :
151 1 : kappa, q = 2, 0
152 1 : wigner_element = fetch_wigner_element(
153 : f_initial, f_final, m_initial, m_final, kappa, q, connection, parquet_files["misc_wigner"]
154 : )
155 1 : assert np.isclose(
156 : wigner_element,
157 : float((-1) ** (f_final - m_final) * wigner_3j(f_final, kappa, f_initial, -m_final, q, m_initial)),
158 : )
159 :
160 1 : q = fetch_reduced_matrix_element(id_initial, id_final, connection, parquet_files["Yb174_mqdt_matrix_elements_q"])
161 1 : matrix_element -= 1 / 12 * wigner_element * q * (bfield_in_gauss * GAUSS_IN_ATOMIC_UNITS) ** 2 * HARTREE_IN_GHZ
162 1 : assert np.isclose(matrix_element, operator[0, 0] if swap_states else operator[1, 1])
163 :
164 :
165 1 : @pytest.mark.parametrize("species", SUPPORTED_SPECIES)
166 1 : def test_obtaining_kets(species: str) -> None:
167 : """Test obtaining kets from the database."""
168 1 : is_mqdt = species.endswith("_mqdt")
169 1 : is_single_valence_electron = species in ["Rb"]
170 1 : is_triplet = species in ["Sr88_triplet"]
171 :
172 1 : quantum_number_i = SPECIES_TO_NUCLEAR_SPIN[species] if is_mqdt else 0
173 1 : quantum_number_s = 0.5 if is_single_valence_electron else (1 if is_triplet else 0)
174 1 : quantum_number_f = quantum_number_i + quantum_number_s
175 1 : quantum_number_m = quantum_number_i + quantum_number_s
176 :
177 : # Obtain a ket from the database
178 1 : ket = pi.KetAtom(species, n=60, l=0, f=quantum_number_f, m=quantum_number_m, s=quantum_number_s)
179 :
180 1 : print(f"|{ket}>")
181 :
182 : # Check the result
183 1 : assert ket.species == species
184 1 : assert ket.n == 60 if not is_mqdt else abs(ket.n - 60) < 1
185 1 : assert ket.l == 0 if not is_mqdt else abs(ket.l - 0) < 1
186 1 : assert ket.f == quantum_number_f
187 1 : assert ket.m == quantum_number_m
188 1 : assert ket.s == quantum_number_s if not is_mqdt else abs(ket.s - quantum_number_s) < 1
189 1 : assert {
190 : "Rb": "Rb:60,S_1/2,1/2",
191 : "Sr88_singlet": "Sr88_singlet:60,S_0,0",
192 : "Sr88_triplet": "Sr88_triplet:60,S_1,1",
193 : "Sr88_mqdt": "Sr88:S=0.0,nu=56.7,L=0.0,J=0,0",
194 : "Sr87_mqdt": "Sr87:S=0.5,nu=56.6,L=0.0,F=9/2,9/2",
195 : "Yb174_mqdt": "Yb174:S=0.0,nu=55.7,L=0.0,J=0,0",
196 : "Yb171_mqdt": "Yb171:S=0.3,nu=55.7,L=0.0,F=1/2,1/2",
197 : "Yb173_mqdt": "Yb173:S=0.4,nu=56.4,L=0.0,F=5/2,5/2",
198 : }[species] == repr(ket)
|
