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 abc import ABC, abstractmethod
7 1 : from typing import TYPE_CHECKING, overload
8 :
9 1 : import numpy as np
10 :
11 1 : from pairinteraction.units import QuantityArray, QuantityScalar, ureg
12 :
13 : if TYPE_CHECKING:
14 : from collections.abc import Callable, Collection
15 :
16 : from pairinteraction.green_tensor.green_tensor_interpolator import GreenTensorInterpolator
17 : from pairinteraction.units import (
18 : ArrayLike,
19 : Dimension,
20 : NDArray,
21 : PintArray, # needed for sphinx to recognize PintArrayLike
22 : PintArrayLike,
23 : PintFloat,
24 : )
25 :
26 :
27 1 : class GreenTensorBase(ABC):
28 1 : epsilon: complex | Callable[[PintFloat], complex]
29 :
30 1 : def __init__(
31 : self,
32 : pos1: ArrayLike | PintArrayLike,
33 : pos2: ArrayLike | PintArrayLike,
34 : unit: str | None = None,
35 : static_limit: bool = True,
36 : interaction_order: int = 3,
37 : *,
38 : without_vacuum_contribution: bool = False,
39 : ) -> None:
40 1 : self.pos1_au = np.array([QuantityScalar.convert_user_to_au(v, unit, "distance") for v in pos1])
41 1 : self.pos2_au = np.array([QuantityScalar.convert_user_to_au(v, unit, "distance") for v in pos2])
42 1 : self.static_limit = static_limit
43 1 : self.without_vacuum_contribution = without_vacuum_contribution
44 :
45 1 : self.interaction_order = interaction_order
46 1 : if interaction_order != 3:
47 0 : raise NotImplementedError(
48 : "Only interaction order 3 (dipole-dipole) is currently implemented for custom Green tensors."
49 : )
50 :
51 1 : self.epsilon = 1.0
52 :
53 : @overload
54 : def get(
55 : self,
56 : kappa1: int,
57 : kappa2: int,
58 : transition_energy: float | PintFloat,
59 : transition_energy_unit: str | None = None,
60 : unit: None = None,
61 : *,
62 : scaled: bool = False,
63 : ) -> PintArray: ...
64 :
65 : @overload
66 : def get(
67 : self,
68 : kappa1: int,
69 : kappa2: int,
70 : transition_energy: float,
71 : transition_energy_unit: str,
72 : unit: str,
73 : *,
74 : scaled: bool = False,
75 : ) -> NDArray: ...
76 :
77 : @overload
78 : def get(
79 : self,
80 : kappa1: int,
81 : kappa2: int,
82 : transition_energy: PintFloat,
83 : *,
84 : unit: str,
85 : scaled: bool = False,
86 : ) -> NDArray: ...
87 :
88 1 : def get(
89 : self,
90 : kappa1: int,
91 : kappa2: int,
92 : transition_energy: float | PintFloat,
93 : transition_energy_unit: str | None = None,
94 : unit: str | None = None,
95 : *,
96 : scaled: bool = False,
97 : ) -> PintArray | NDArray:
98 : """Calculate the Green tensor in cartesian coordinates for the given ranks kappa1, kappa2 and frequency omega.
99 :
100 : kappa = 1 corresponds to dipole operator with the cartesian basis: [x, y, z]
101 :
102 : Args:
103 : kappa1: The rank of the first multipole operator.
104 : kappa2: The rank of the second multipole operator.
105 : transition_energy: The transition energy at which to evaluate the Green tensor.
106 : Only needed if the Green tensor is frequency dependent.
107 : transition_energy_unit: The unit of the transition energy.
108 : Default None, which means that the transition energy must be given as pint object.
109 : unit: The unit to which to convert the result.
110 : Default None, which means that the result is returned as pint object.
111 : scaled: If True, the Green tensor is returned with the prefactor for the interaction
112 : already included (the unit has to be adopted accordingly).
113 : Default False, which means that the bare Green tensor is returned.
114 :
115 : Returns:
116 : The Green tensor as a 2D array in cartesian coordinates.
117 :
118 : """
119 1 : omega_au = QuantityScalar.convert_user_to_au(transition_energy, transition_energy_unit, "energy")
120 1 : omega_for_calculation = 0 if self.static_limit else omega_au
121 :
122 1 : scaled_gt_au = self._get_scaled_au(kappa1, kappa2, omega_for_calculation)
123 :
124 1 : prefactor = 1 if scaled else self._get_prefactor_au(kappa1, kappa2, omega_au)
125 1 : dimension = self._get_dimension(kappa1, kappa2, scaled)
126 1 : return QuantityArray.convert_au_to_user(scaled_gt_au / prefactor, dimension, unit)
127 :
128 : @abstractmethod
129 : def _get_scaled_au(self, kappa1: int, kappa2: int, transition_energy_au: float) -> NDArray: ...
130 :
131 1 : @staticmethod
132 1 : def _get_prefactor_au(kappa1: int, kappa2: int, transition_energy_au: float) -> float:
133 : r"""Get the prefactor to get the interaction strength from the Green tensor.
134 :
135 : The interaction between two dipole moments is given as (see e.g. https://arxiv.org/pdf/2303.13564)
136 : .. math::
137 : V_{\alpha\beta} = \frac{\omega^2}{\hbar \epsilon_0 c^2}
138 : d_\alpha^T \mathrm{Re}\{G(r_\alpha, r_\beta, \omega)\} d_\beta
139 :
140 : This functions returns the prefactor
141 : .. math::
142 : \frac{\omega^2}{\hbar \epsilon_0 c^2}
143 :
144 : In C++ we use the convention that the Green tensor already contains this prefactor.
145 :
146 : Args:
147 : kappa1: The rank of the first multipole operator.
148 : kappa2: The rank of the second multipole operator.
149 : transition_energy_au: The transition energy at which to evaluate the Green tensor in atomic units (hartree).
150 :
151 : """
152 1 : if kappa1 == 1 and kappa2 == 1:
153 1 : speed_of_light_au: float = ureg.Quantity(1, "speed_of_light").to_base_units().m
154 1 : prefactor = (transition_energy_au / speed_of_light_au) ** 2
155 1 : prefactor /= 1 / (4 * np.pi) # hbar = 1, epsilon_0 = 1 / (4*np.pi) in atomic units
156 1 : prefactor *= -1 # minus sign from H = - \hbar \sum V_{\alpha\beta} ...
157 1 : return -prefactor
158 0 : raise NotImplementedError("Only dipole-dipole Green tensor prefactor is currently implemented.")
159 :
160 1 : @staticmethod
161 1 : def _get_dimension(kappa1: int, kappa2: int, scaled: bool) -> list[Dimension]:
162 1 : dimension_00: Dimension = "scaled_green_tensor_00" if scaled else "green_tensor_00"
163 1 : dimension: list[Dimension] = [dimension_00] + ["inverse_distance"] * (kappa1 + kappa2)
164 1 : return dimension
165 :
166 : @overload
167 : def get_interpolator(self, *, use_real: bool) -> GreenTensorInterpolator: ...
168 :
169 : @overload
170 : def get_interpolator(
171 : self,
172 : transition_energies: Collection[PintFloat] | PintArray,
173 : transition_energies_unit: None = None,
174 : *,
175 : use_real: bool,
176 : ) -> GreenTensorInterpolator: ...
177 :
178 : @overload
179 : def get_interpolator(
180 : self, transition_energies: Collection[float] | NDArray, transition_energies_unit: str, *, use_real: bool
181 : ) -> GreenTensorInterpolator: ...
182 :
183 1 : def get_interpolator(
184 : self,
185 : transition_energies: Collection[PintFloat] | PintArray | Collection[float] | NDArray | None = None,
186 : transition_energies_unit: str | None = None,
187 : *,
188 : use_real: bool,
189 : ) -> GreenTensorInterpolator:
190 : """Get a GreenTensorInterpolator from this Green tensor.
191 :
192 : The GreenTensorInterpolator can be used for the interaction of a SystemPair.
193 :
194 : Returns:
195 : A GreenTensorInterpolator that interpolates the Green tensor at given frequency points.
196 :
197 : """
198 : GTIClass: type[GreenTensorInterpolator] # noqa: N806
199 1 : if use_real:
200 1 : from pairinteraction.green_tensor.green_tensor_interpolator import GreenTensorInterpolatorReal as GTIClass
201 : else:
202 1 : from pairinteraction.green_tensor.green_tensor_interpolator import GreenTensorInterpolator as GTIClass
203 :
204 1 : if self.static_limit and not (transition_energies is None and transition_energies_unit is None):
205 0 : raise ValueError("You must not specify a frequency range when static limit is set.")
206 :
207 1 : if self.static_limit:
208 1 : gti = GTIClass()
209 1 : scaled_gt_au = self.get(1, 1, 0, scaled=True)
210 1 : gti.set_constant(1, 1, scaled_gt_au)
211 1 : return gti
212 :
213 0 : if transition_energies is not None:
214 0 : omegas_pint = [
215 : QuantityScalar.convert_user_to_pint(omega, transition_energies_unit, "energy")
216 : for omega in transition_energies
217 : ]
218 0 : gti = GTIClass()
219 0 : scaled_gt_list = [self.get(1, 1, omega, scaled=True) for omega in omegas_pint]
220 0 : gti.set_list(1, 1, scaled_gt_list, omegas_pint, scaled=True)
221 0 : return gti
222 :
223 0 : raise ValueError(
224 : "You must either specify transition_energies or set static_limit to True to get an interpolator object."
225 : )
|
