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], sort_by_energy=True)
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((system, system), energy=(pair_energy - 0.1, pair_energy + 0.1), energy_unit="GHz")
104 :
105 1 : dd = ket1.get_matrix_element(ket2, "electric_dipole", q=0)
106 1 : gt_reference = reference_green_tensor_vacuum(distance_vector_mum)
107 1 : reference = dd * gt_reference[2, 2] * dd
108 :
109 : # test internal vacuum green tensor
110 1 : system_pair = pi_module.SystemPair(basis_pair)
111 1 : system_pair.set_distance_vector(distance_vector_mum, "micrometer")
112 1 : hamiltonian = system_pair.get_hamiltonian("hartree").toarray()
113 1 : np.testing.assert_allclose(hamiltonian[1, 0], reference)
114 :
115 : # test custom free space vacuum green tensor
116 1 : system_pair = pi_module.SystemPair(basis_pair)
117 1 : gt = GreenTensorFreeSpace([0, 0, 0], distance_vector_mum, unit="micrometer", static_limit=True)
118 1 : system_pair.set_green_tensor(gt)
119 1 : hamiltonian = system_pair.get_hamiltonian("hartree").toarray()
120 1 : np.testing.assert_allclose(hamiltonian[1, 0], reference)
121 :
122 :
123 1 : def reference_green_tensor_vacuum(distance_vector_mum: list[float]) -> NDArray:
124 1 : distance_mum = np.linalg.norm(distance_vector_mum)
125 1 : distance_au = ureg.Quantity(distance_mum, "micrometer").to_base_units().m
126 1 : gt: NDArray = (
127 : np.eye(3) - 3 * np.outer(distance_vector_mum, distance_vector_mum) / distance_mum**2
128 : ) / distance_au**3
129 1 : return gt
130 :
131 :
132 1 : def rotation_matrix_from_axis_angle(axis: list[float], angle: float) -> NDArray:
133 1 : axis_array = np.array(axis, dtype=float)
134 1 : axis_array /= np.linalg.norm(axis_array)
135 1 : cross_matrix = np.array(
136 : [
137 : [0, -axis_array[2], axis_array[1]],
138 : [axis_array[2], 0, -axis_array[0]],
139 : [-axis_array[1], axis_array[0], 0],
140 : ]
141 : )
142 1 : return np.eye(3) + np.sin(angle) * cross_matrix + (1 - np.cos(angle)) * cross_matrix @ cross_matrix # type: ignore [no-any-return]
143 :
144 :
145 1 : def test_surface_green_tensor_rotation_invariance() -> None:
146 1 : rotation = rotation_matrix_from_axis_angle([1, 2, 0.5], np.deg2rad(37))
147 1 : pos1 = np.array([0.4, -0.2, 0.8])
148 1 : pos2 = np.array([1.1, 0.3, 1.4])
149 1 : point_on_plane = np.array([0.0, 0.0, 0.0])
150 :
151 1 : gt_reference = GreenTensorSurface(
152 : pos1,
153 : pos2,
154 : point_on_plane=point_on_plane,
155 : surface_normal=[0, 0, 1],
156 : unit="micrometer",
157 : static_limit=True,
158 : ).get(1, 1, transition_energy=0, scaled=True)
159 1 : gt_rotated = GreenTensorSurface(
160 : rotation.T @ pos1,
161 : rotation.T @ pos2,
162 : point_on_plane=rotation.T @ point_on_plane,
163 : surface_normal=rotation.T @ np.array([0, 0, 1]),
164 : unit="micrometer",
165 : static_limit=True,
166 : ).get(1, 1, transition_energy=0, scaled=True)
167 :
168 1 : np.testing.assert_allclose(
169 : gt_rotated.m,
170 : rotation.T @ gt_reference.m @ rotation,
171 : rtol=1e-10,
172 : atol=1e-12,
173 : )
174 :
175 :
176 1 : def test_cavity_green_tensor_rotation_invariance() -> None:
177 1 : rotation = rotation_matrix_from_axis_angle([0.3, 1.0, -0.2], np.deg2rad(51))
178 1 : pos1 = np.array([0.2, -0.1, -0.2])
179 1 : pos2 = np.array([0.8, 0.5, 0.4])
180 1 : point_on_plane1 = np.array([0.0, 0.0, -1.0])
181 1 : point_on_plane2 = np.array([0.0, 0.0, 1.5])
182 :
183 1 : gt_reference = GreenTensorCavity(
184 : pos1,
185 : pos2,
186 : point_on_plane1=point_on_plane1,
187 : point_on_plane2=point_on_plane2,
188 : surface_normal=[0, 0, 1],
189 : unit="micrometer",
190 : static_limit=True,
191 : ).get(1, 1, transition_energy=0, scaled=True)
192 1 : gt_rotated = GreenTensorCavity(
193 : rotation.T @ pos1,
194 : rotation.T @ pos2,
195 : point_on_plane1=rotation.T @ point_on_plane1,
196 : point_on_plane2=rotation.T @ point_on_plane2,
197 : surface_normal=rotation.T @ np.array([0, 0, 1]),
198 : unit="micrometer",
199 : static_limit=True,
200 : ).get(1, 1, transition_energy=0, scaled=True)
201 :
202 1 : np.testing.assert_allclose(
203 : gt_rotated.m,
204 : rotation.T @ gt_reference.m @ rotation,
205 : rtol=1e-10,
206 : atol=1e-12,
207 : )
208 :
209 :
210 1 : def test_surface_green_tensor_rejects_zero_normal() -> None:
211 1 : with pytest.raises(ValueError, match="cannot be zero"):
212 1 : GreenTensorSurface(
213 : [0, 0, 0],
214 : [0, 0, 1],
215 : point_on_plane=[0, 0, 0],
216 : surface_normal=[0, 0, 0],
217 : unit="micrometer",
218 : static_limit=True,
219 : )
220 :
221 :
222 1 : def test_cavity_green_tensor_rejects_coincident_planes() -> None:
223 1 : gt = GreenTensorCavity(
224 : [0, 0, 0],
225 : [0, 0, 1],
226 : point_on_plane1=[0, 0, 0],
227 : point_on_plane2=[1, 0, 0],
228 : surface_normal=[0, 0, 1],
229 : unit="micrometer",
230 : static_limit=True,
231 : )
232 :
233 1 : with pytest.raises(ValueError, match="must be distinct"):
234 1 : gt.get(1, 1, transition_energy=0, scaled=True)
|
