Line data Source code
1 : # SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : # SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 1 : from __future__ import annotations
5 :
6 1 : from typing import TYPE_CHECKING
7 :
8 1 : import numpy as np
9 1 : import pytest
10 1 : from pairinteraction import ureg
11 1 : from pairinteraction.green_tensor import GreenTensorFreeSpace
12 1 : from pairinteraction.green_tensor.green_tensor_cavity import GreenTensorCavity
13 1 : from pairinteraction.green_tensor.green_tensor_surface import GreenTensorSurface
14 :
15 : if TYPE_CHECKING:
16 : from pairinteraction.green_tensor.green_tensor_base import GreenTensorBase
17 : from pairinteraction.units import NDArray
18 :
19 : from .utils import PairinteractionModule
20 :
21 1 : DISTANCE_VECTOR_MUM_LIST = [[1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 0, 3.5]]
22 :
23 :
24 1 : @pytest.mark.parametrize("distance_vector_mum", DISTANCE_VECTOR_MUM_LIST)
25 1 : def test_static_green_tensor(distance_vector_mum: list[float]) -> None:
26 1 : gt_reference = reference_green_tensor_vacuum(distance_vector_mum)
27 :
28 : gt: GreenTensorBase
29 1 : gt = GreenTensorFreeSpace([0, 0, 0], distance_vector_mum, unit="micrometer", static_limit=True)
30 1 : gt_dipole_dipole = gt.get(1, 1, transition_energy=0, scaled=True)
31 1 : np.testing.assert_allclose(gt_dipole_dipole.m, gt_reference)
32 :
33 1 : gt = GreenTensorSurface(
34 : [0, 0, 0],
35 : distance_vector_mum,
36 : point_on_plane=[0, 0, 1000],
37 : surface_normal=[0, 0, 1],
38 : unit="micrometer",
39 : static_limit=True,
40 : )
41 1 : gt_dipole_dipole = gt.get(1, 1, transition_energy=0, scaled=True)
42 1 : np.testing.assert_allclose(gt_dipole_dipole.m, gt_reference, rtol=1e-6, atol=1e-16)
43 :
44 1 : gt = GreenTensorCavity(
45 : [0, 0, 0],
46 : distance_vector_mum,
47 : point_on_plane1=[0, 0, -1000],
48 : point_on_plane2=[0, 0, 1000],
49 : surface_normal=[0, 0, 1],
50 : unit="micrometer",
51 : static_limit=True,
52 : )
53 1 : gt_dipole_dipole = gt.get(1, 1, transition_energy=0, scaled=True)
54 1 : np.testing.assert_allclose(gt_dipole_dipole.m, gt_reference, rtol=1e-6, atol=1e-16)
55 :
56 :
57 1 : @pytest.mark.parametrize("distance_vector_mum", DISTANCE_VECTOR_MUM_LIST)
58 1 : def test_static_green_tensor_pair_potential(
59 : pi_module: PairinteractionModule, use_real: bool, distance_vector_mum: list[float]
60 : ) -> None:
61 1 : if use_real and distance_vector_mum[1] != 0:
62 1 : return # run tests with y-component only for complex Green tensors
63 :
64 1 : ket = pi_module.KetAtom("Rb", n=60, l=0, m=0.5)
65 1 : basis = pi_module.BasisAtom("Rb", n=(ket.n - 2, ket.n + 2), l=(0, 2))
66 1 : system = pi_module.SystemAtom(basis)
67 1 : system.set_magnetic_field([0, 0, 1], "gauss")
68 1 : system.set_electric_field([0, 0, 1], "V/cm")
69 1 : system.diagonalize()
70 :
71 1 : pair_energy_ghz = 2 * system.get_corresponding_energy(ket, unit="GHz")
72 1 : energy_range_ghz = (pair_energy_ghz - 3, pair_energy_ghz + 3)
73 1 : basis_pair = pi_module.BasisPair([system, system], energy=energy_range_ghz, energy_unit="GHz")
74 :
75 1 : system_pair_vacuum = pi_module.SystemPair(basis_pair)
76 1 : system_pair_vacuum.set_distance_vector(distance_vector_mum, unit="micrometer")
77 1 : system_pair_vacuum.set_interaction_order(3)
78 :
79 1 : gt = GreenTensorFreeSpace([0, 0, 0], distance_vector_mum, unit="micrometer", interaction_order=3, static_limit=True)
80 1 : gt.set_relative_permittivity(1.0)
81 1 : system_pair_free_space = pi_module.SystemPair(basis_pair).set_green_tensor(gt)
82 :
83 1 : pi_module.diagonalize([system_pair_free_space, system_pair_vacuum])
84 1 : np.testing.assert_allclose(
85 : system_pair_vacuum.get_eigenenergies("GHz") - pair_energy_ghz,
86 : system_pair_free_space.get_eigenenergies("GHz") - pair_energy_ghz,
87 : )
88 :
89 :
90 1 : @pytest.mark.parametrize("distance_vector_mum", DISTANCE_VECTOR_MUM_LIST)
91 1 : def test_vacuum_green_tensor(
92 : pi_module: PairinteractionModule, use_real: bool, distance_vector_mum: list[float]
93 : ) -> None:
94 1 : if use_real and distance_vector_mum[1] != 0:
95 1 : return # run tests with y-component only for complex Green tensors
96 :
97 1 : ket1 = pi_module.KetAtom("Rb", n=60, l=0, j=0.5, m=0.5)
98 1 : ket2 = pi_module.KetAtom("Rb", n=60, l=1, j=0.5, m=0.5)
99 :
100 1 : basis = pi_module.BasisAtom("Rb", n=(0, 0), additional_kets=[ket1, ket2])
101 1 : system = pi_module.SystemAtom(basis)
102 1 : pair_energy = ket1.get_energy("GHz") + ket2.get_energy("GHz")
103 1 : basis_pair = pi_module.BasisPair(
104 : (system, system), energy=(pair_energy - 0.1, pair_energy + 0.1), energy_unit="GHz", symmetrize=False
105 : )
106 :
107 1 : dd = ket1.get_matrix_element(ket2, "electric_dipole", q=0)
108 1 : gt_reference = reference_green_tensor_vacuum(distance_vector_mum)
109 1 : reference = dd * gt_reference[2, 2] * dd
110 :
111 : # test internal vacuum green tensor
112 1 : system_pair = pi_module.SystemPair(basis_pair)
113 1 : system_pair.set_distance_vector(distance_vector_mum, "micrometer")
114 1 : hamiltonian = system_pair.get_hamiltonian("hartree").toarray()
115 1 : np.testing.assert_allclose(hamiltonian[1, 0], reference)
116 :
117 : # test custom free space vacuum green tensor
118 1 : system_pair = pi_module.SystemPair(basis_pair)
119 1 : gt = GreenTensorFreeSpace([0, 0, 0], distance_vector_mum, unit="micrometer", static_limit=True)
120 1 : system_pair.set_green_tensor(gt)
121 1 : hamiltonian = system_pair.get_hamiltonian("hartree").toarray()
122 1 : np.testing.assert_allclose(hamiltonian[1, 0], reference)
123 :
124 :
125 1 : def reference_green_tensor_vacuum(distance_vector_mum: list[float]) -> NDArray:
126 1 : distance_mum = np.linalg.norm(distance_vector_mum)
127 1 : distance_au = ureg.Quantity(distance_mum, "micrometer").to_base_units().m
128 1 : gt: NDArray = (
129 : np.eye(3) - 3 * np.outer(distance_vector_mum, distance_vector_mum) / distance_mum**2
130 : ) / distance_au**3
131 1 : return gt
132 :
133 :
134 1 : def rotation_matrix_from_axis_angle(axis: list[float], angle: float) -> NDArray:
135 1 : axis_array = np.array(axis, dtype=float)
136 1 : axis_array /= np.linalg.norm(axis_array)
137 1 : cross_matrix = np.array(
138 : [
139 : [0, -axis_array[2], axis_array[1]],
140 : [axis_array[2], 0, -axis_array[0]],
141 : [-axis_array[1], axis_array[0], 0],
142 : ]
143 : )
144 1 : return np.eye(3) + np.sin(angle) * cross_matrix + (1 - np.cos(angle)) * cross_matrix @ cross_matrix # type: ignore [no-any-return]
145 :
146 :
147 1 : def test_surface_green_tensor_rotation_invariance() -> None:
148 1 : rotation = rotation_matrix_from_axis_angle([1, 2, 0.5], np.deg2rad(37))
149 1 : pos1 = np.array([0.4, -0.2, 0.8])
150 1 : pos2 = np.array([1.1, 0.3, 1.4])
151 1 : point_on_plane = np.array([0.0, 0.0, 0.0])
152 :
153 1 : gt_reference = GreenTensorSurface(
154 : pos1,
155 : pos2,
156 : point_on_plane=point_on_plane,
157 : surface_normal=[0, 0, 1],
158 : unit="micrometer",
159 : static_limit=True,
160 : ).get(1, 1, transition_energy=0, scaled=True)
161 1 : gt_rotated = GreenTensorSurface(
162 : rotation.T @ pos1,
163 : rotation.T @ pos2,
164 : point_on_plane=rotation.T @ point_on_plane,
165 : surface_normal=rotation.T @ np.array([0, 0, 1]),
166 : unit="micrometer",
167 : static_limit=True,
168 : ).get(1, 1, transition_energy=0, scaled=True)
169 :
170 1 : np.testing.assert_allclose(
171 : gt_rotated.m,
172 : rotation.T @ gt_reference.m @ rotation,
173 : rtol=1e-10,
174 : atol=1e-12,
175 : )
176 :
177 :
178 1 : def test_cavity_green_tensor_rotation_invariance() -> None:
179 1 : rotation = rotation_matrix_from_axis_angle([0.3, 1.0, -0.2], np.deg2rad(51))
180 1 : pos1 = np.array([0.2, -0.1, -0.2])
181 1 : pos2 = np.array([0.8, 0.5, 0.4])
182 1 : point_on_plane1 = np.array([0.0, 0.0, -1.0])
183 1 : point_on_plane2 = np.array([0.0, 0.0, 1.5])
184 :
185 1 : gt_reference = GreenTensorCavity(
186 : pos1,
187 : pos2,
188 : point_on_plane1=point_on_plane1,
189 : point_on_plane2=point_on_plane2,
190 : surface_normal=[0, 0, 1],
191 : unit="micrometer",
192 : static_limit=True,
193 : ).get(1, 1, transition_energy=0, scaled=True)
194 1 : gt_rotated = GreenTensorCavity(
195 : rotation.T @ pos1,
196 : rotation.T @ pos2,
197 : point_on_plane1=rotation.T @ point_on_plane1,
198 : point_on_plane2=rotation.T @ point_on_plane2,
199 : surface_normal=rotation.T @ np.array([0, 0, 1]),
200 : unit="micrometer",
201 : static_limit=True,
202 : ).get(1, 1, transition_energy=0, scaled=True)
203 :
204 1 : np.testing.assert_allclose(
205 : gt_rotated.m,
206 : rotation.T @ gt_reference.m @ rotation,
207 : rtol=1e-10,
208 : atol=1e-12,
209 : )
210 :
211 :
212 1 : def test_surface_green_tensor_rejects_zero_normal() -> None:
213 1 : with pytest.raises(ValueError, match="cannot be zero"):
214 1 : GreenTensorSurface(
215 : [0, 0, 0],
216 : [0, 0, 1],
217 : point_on_plane=[0, 0, 0],
218 : surface_normal=[0, 0, 0],
219 : unit="micrometer",
220 : static_limit=True,
221 : )
222 :
223 :
224 1 : def test_cavity_green_tensor_rejects_coincident_planes() -> None:
225 1 : gt = GreenTensorCavity(
226 : [0, 0, 0],
227 : [0, 0, 1],
228 : point_on_plane1=[0, 0, 0],
229 : point_on_plane2=[1, 0, 0],
230 : surface_normal=[0, 0, 1],
231 : unit="micrometer",
232 : static_limit=True,
233 : )
234 :
235 1 : with pytest.raises(ValueError, match="must be distinct"):
236 1 : gt.get(1, 1, transition_energy=0, scaled=True)
|
