Line data Source code
1 : # SPDX-FileCopyrightText: 2025 Pairinteraction Developers
2 : # SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 1 : import logging
5 :
6 1 : import numpy as np
7 1 : import pytest
8 1 : from scipy import sparse
9 :
10 1 : import pairinteraction.real as pi
11 1 : from pairinteraction import perturbative
12 1 : from pairinteraction.perturbative.perturbative import _calculate_perturbative_hamiltonian
13 :
14 :
15 : # ---------------------------------------------------------------------------------------
16 : # Helper functions
17 : # ---------------------------------------------------------------------------------------
18 1 : def _check_sparse_matrices_equal(matrix_a: sparse.csr_matrix, matrix_b: sparse.csr_matrix) -> bool:
19 : """Check for equality of sparse matrices efficiently.
20 :
21 : This functions compares two sparse matrices in compressed sparse row format on their equality.
22 :
23 : Args:
24 : matrix_a: A sparse matrix in csr format.
25 : matrix_b: A sparse matrix in csr format.
26 :
27 : Returns:
28 : bool: True if matrices are equal, False if not.
29 :
30 : """
31 1 : matrix_a.sort_indices()
32 1 : matrix_b.sort_indices()
33 1 : if (
34 : matrix_a.format == "csr"
35 : and matrix_b.format == "csr"
36 : and len(matrix_a.indices) == len(matrix_b.indices)
37 : and len(matrix_a.indptr) == len(matrix_b.indptr)
38 : and len(matrix_a.data) == len(matrix_b.data)
39 : ):
40 1 : return (
41 : np.all(matrix_a.indices == matrix_b.indices)
42 : and np.all(matrix_a.indptr == matrix_b.indptr)
43 : and np.allclose(matrix_a.data, matrix_b.data, rtol=0, atol=1e-14)
44 : )
45 0 : return False
46 :
47 :
48 1 : def _create_system_pair_sample():
49 1 : basis = pi.BasisAtom(
50 : species="Rb",
51 : n=(59, 63),
52 : l=(0, 1),
53 : m=(-1.5, 1.5),
54 : )
55 1 : system = pi.SystemAtom(basis=basis)
56 1 : system.set_diamagnetism_enabled(False)
57 1 : system.set_magnetic_field([0, 0, 1e-3], "gauss")
58 1 : if not system.is_diagonal:
59 1 : pi.diagonalize([system], diagonalizer="eigen", sort_by_energy=False)
60 1 : basis_pair = pi.BasisPair([system, system])
61 1 : system_pair = pi.SystemPair(basis_pair)
62 1 : theta = 0
63 1 : r = 12
64 1 : system_pair.set_distance_vector(r * np.array([np.sin(theta), 0, np.cos(theta)]), "micrometer")
65 1 : system_pair.set_interaction_order(3)
66 1 : return system_pair
67 :
68 :
69 : # ---------------------------------------------------------------------------------------
70 : # Tests
71 : # ---------------------------------------------------------------------------------------
72 1 : def test_perturbative_calculation(caplog):
73 : """Test of mathematic functionality."""
74 1 : H = sparse.csr_matrix([[0, 1, 1, 0, 2], [1, 1, 0, 1, 3], [1, 0, 10, 1, 0], [0, 1, 1, 11, 1], [2, 3, 0, 1, 12]])
75 1 : model_space_indices = [0, 1]
76 :
77 1 : hamiltonian_eff, eig_perturb = _calculate_perturbative_hamiltonian(H, model_space_indices, order=0)
78 :
79 1 : assert np.any(hamiltonian_eff == np.array([[0, 0], [0, 1]]))
80 1 : assert _check_sparse_matrices_equal(eig_perturb, sparse.eye(2, 5, k=0, format="csr"))
81 :
82 1 : hamiltonian_eff, eig_perturb = _calculate_perturbative_hamiltonian(H, model_space_indices, order=1)
83 :
84 1 : assert np.any(hamiltonian_eff == np.array([[0, 1], [1, 1]]))
85 1 : assert _check_sparse_matrices_equal(eig_perturb, sparse.csr_matrix([[1, 0, 0, 0, 0], [0, 1, 0, 0, 0]]))
86 :
87 1 : hamiltonian_eff, eig_perturb = _calculate_perturbative_hamiltonian(H, model_space_indices, order=2)
88 1 : a_00 = 0 + 1 * 1 / (0 - 10) + 2 * 2 / (0 - 12)
89 1 : a_01 = 1 + 2 * 3 / (0 - 12)
90 1 : a_10 = 1 + 3 * 2 / (1 - 12)
91 1 : a_11 = 1 + 1 * 1 / (1 - 11) + 3 * 3 / (1 - 12)
92 1 : hamiltonian_new = np.array([[a_00, (a_01 + a_10) / 2], [(a_01 + a_10) / 2, a_11]])
93 1 : assert np.any(hamiltonian_eff == hamiltonian_new)
94 :
95 1 : v0 = (
96 : sparse.eye(1, 5, k=0, format="csr")
97 : + 1 / (0 - 10) * sparse.eye(1, 5, k=2, format="csr")
98 : + 2 / (0 - 12) * sparse.eye(1, 5, k=4, format="csr")
99 : )
100 1 : v1 = (
101 : sparse.eye(1, 5, k=1, format="csr")
102 : + 1 / (1 - 11) * sparse.eye(1, 5, k=3, format="csr")
103 : + 3 / (1 - 12) * sparse.eye(1, 5, k=4, format="csr")
104 : )
105 1 : assert _check_sparse_matrices_equal(eig_perturb, sparse.vstack([v0, v1]))
106 :
107 1 : with caplog.at_level(logging.ERROR):
108 1 : hamiltonian_eff, eig_perturb = _calculate_perturbative_hamiltonian(H, model_space_indices, order=3)
109 1 : a_00 -= 2 * 3 * 1 / ((0 - 12) * (1 - 12))
110 1 : a_01 += (
111 : 1 * 1 * 1 / ((1 - 10) * (1 - 11))
112 : + 2 * 1 * 1 / ((1 - 11) * (1 - 12))
113 : - 1 * 1 * 1 / ((0 - 10) * (1 - 10))
114 : - 2 * 2 * 1 / ((1 - 12) * (0 - 12))
115 : )
116 1 : a_10 += (
117 : 1 * 1 * 1 / ((0 - 10) * (0 - 11))
118 : + 2 * 1 * 1 / ((0 - 11) * (0 - 12))
119 : - 1 * 1 * 1 / ((0 - 11) * (1 - 11))
120 : - 3 * 3 * 1 / ((0 - 12) * (1 - 12))
121 : )
122 1 : a_11 += 1 * 1 * 3 / ((1 - 11) * (1 - 12)) + 3 * 1 * 1 / ((1 - 11) * (1 - 12)) - 3 * 2 * 1 / ((1 - 12) * (0 - 12))
123 1 : hamiltonian_new = sparse.csr_matrix([[a_00, (a_01 + a_10) / 2], [(a_01 + a_10) / 2, a_11]])
124 1 : assert np.any(hamiltonian_eff == hamiltonian_new)
125 :
126 1 : v0 += (
127 : -0.5 * (1 * 1 / (0 - 10) ** 2 + 2 * 2 / (0 - 12) ** 2) * sparse.eye(1, 5, k=0, format="csr")
128 : + (1 * 1 / ((0 - 10) * (0 - 11)) + 2 * 1 / ((0 - 11) * (0 - 12))) * sparse.eye(1, 5, k=3, format="csr")
129 : + 1 * 1 / ((0 - 11) * (0 - 1)) * sparse.eye(1, 5, k=3, format="csr")
130 : + 3 * 1 / ((0 - 1) * (0 - 12)) * sparse.eye(1, 5, k=4, format="csr")
131 : + 3 * 2 / ((0 - 1) * (0 - 12)) * sparse.eye(1, 5, k=1, format="csr")
132 : )
133 1 : v1 += (
134 : -0.5 * (1 * 1 / (1 - 11) ** 2 + 3 * 3 / (1 - 12) ** 2) * sparse.eye(1, 5, k=1, format="csr")
135 : + 1 * 1 / ((1 - 10) * (1 - 11)) * sparse.eye(1, 5, k=2, format="csr")
136 : + 1 * 3 / ((1 - 11) * (1 - 12)) * sparse.eye(1, 5, k=3, format="csr")
137 : + 1 * 1 / ((1 - 11) * (1 - 12)) * sparse.eye(1, 5, k=4, format="csr")
138 : + 1 * 1 / ((1 - 0) * (1 - 10)) * sparse.eye(1, 5, k=2, format="csr")
139 : + 2 * 1 / ((1 - 0) * (1 - 12)) * sparse.eye(1, 5, k=4, format="csr")
140 : + 3 * 2 / ((1 - 0) * (1 - 12)) * sparse.eye(1, 5, k=0, format="csr")
141 : )
142 1 : assert _check_sparse_matrices_equal(eig_perturb, sparse.vstack([v0, v1]))
143 :
144 :
145 1 : def test_c3_coefficient():
146 : """Test whether dispersion coefficients are correctly calculated."""
147 1 : system_pair = _create_system_pair_sample()
148 :
149 1 : ket_tuple_list = [
150 : (pi.KetAtom("Rb", n=61, l=0, j=0.5, m=0.5), pi.KetAtom("Rb", n=61, l=1, j=1.5, m=0.5)),
151 : (pi.KetAtom("Rb", n=61, l=1, j=1.5, m=0.5), pi.KetAtom("Rb", n=61, l=0, j=0.5, m=0.5)),
152 : ]
153 1 : c3 = perturbative.get_c3_from_system(
154 : system_pair=system_pair, ket_tuple_list=ket_tuple_list, unit="planck_constant * gigahertz * micrometer^3"
155 : )
156 1 : assert np.isclose(-0.5 * c3, 3.1515)
157 :
158 :
159 1 : def test_c6_coefficient():
160 : """Test whether the c6 coefficient is correct."""
161 1 : system_pair = _create_system_pair_sample()
162 1 : ket_atom = pi.KetAtom(species="Rb", n=61, l=0, j=0.5, m=0.5)
163 1 : c6 = perturbative.get_c6_from_system(
164 : ket_tuple=(ket_atom, ket_atom), system_pair=system_pair, unit="planck_constant * gigahertz * micrometer^6"
165 : )
166 1 : assert np.isclose(c6, 167.92)
167 :
168 :
169 1 : def test_resonance_detection(caplog):
170 : """Test whether resonance is correctly detected."""
171 1 : system_pair = _create_system_pair_sample()
172 1 : ket_tuple_list = [
173 : (pi.KetAtom(species="Rb", n=61, l=0, j=0.5, m=0.5), pi.KetAtom(species="Rb", n=61, l=1, j=0.5, m=0.5))
174 : ]
175 1 : with (
176 : pytest.warns(RuntimeWarning, match="divide by zero encountered in divide"),
177 : pytest.raises(ValueError, match="Perturbative Calculation not possible due to resonances."),
178 : caplog.at_level(logging.CRITICAL),
179 : ):
180 1 : perturbative.get_effective_hamiltonian_from_system(
181 : ket_tuple_list=ket_tuple_list, system_pair=system_pair, order=2, required_overlap=0.8
182 : )
|
