Line data    Source code

       1             : # SPDX-FileCopyrightText: 2025 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 scipy.constants as const
      10           1 : from typing_extensions import override
      11             : 
      12           1 : from pairinteraction.green_tensor.dynamic_green_tensor import (
      13             :     dynamic_green_tensor_homogeneous,
      14             :     dynamic_green_tensor_scattered,
      15             : )
      16           1 : from pairinteraction.green_tensor.green_tensor_base import GreenTensorBase
      17           1 : from pairinteraction.green_tensor.utils import (
      18             :     evaluate_relative_permittivity,
      19             :     get_lab_to_local_rotation_matrix,
      20             :     rotate_tensor_to_lab,
      21             :     rotate_vector_to_local,
      22             : )
      23           1 : from pairinteraction.units import QuantityScalar, ureg
      24             : 
      25             : if TYPE_CHECKING:
      26             :     from typing_extensions import Self
      27             : 
      28             :     from pairinteraction.green_tensor.utils import PermittivityLike
      29             :     from pairinteraction.units import ArrayLike, NDArray, PintArrayLike
      30             : 
      31             : 
      32           1 : class GreenTensorSurface(GreenTensorBase):
      33             :     """Green tensor for two atoms near a single infinite surface.
      34             : 
      35             :     Examples:
      36             :         >>> from pairinteraction.green_tensor import GreenTensorSurface
      37             :         >>> gt = GreenTensorSurface(
      38             :         ...     [0, 0, 0], [10, 0, 0], point_on_plane=[0, 0, -5], surface_normal=[0, 0, 1], unit="micrometer"
      39             :         ... )
      40             :         >>> transition_energy = 2  # h * GHz
      41             :         >>> gt_dipole_dipole = gt.get(1, 1, transition_energy, "planck_constant * GHz")
      42             :         >>> print(f"{gt_dipole_dipole[0, 0]:.2f}")
      43             :         -4.37 / bohr
      44             : 
      45             :     """
      46             : 
      47           1 :     def __init__(
      48             :         self,
      49             :         pos1: ArrayLike | PintArrayLike,
      50             :         pos2: ArrayLike | PintArrayLike,
      51             :         point_on_plane: ArrayLike | PintArrayLike,
      52             :         surface_normal: ArrayLike,
      53             :         unit: str | None = None,
      54             :         static_limit: bool = True,
      55             :         interaction_order: int = 3,
      56             :         *,
      57             :         without_vacuum_contribution: bool = False,
      58             :     ) -> None:
      59             :         """Create a Green tensor for two atoms near a single infinite surface.
      60             : 
      61             :         The surface is an infinite plane specified by a point on the plane and its normal vector.
      62             :         If not specified otherwise (see `set_relative_permittivities`), the surface is treated as a perfect mirror.
      63             : 
      64             : 
      65             :         Args:
      66             :             pos1: Position of the first atom in the given unit.
      67             :             pos2: Position of the second atom in the given unit.
      68             :             point_on_plane: A point on the surface plane in the given unit.
      69             :             surface_normal: The surface normal vector.
      70             :             unit: The unit of the distance, e.g. "micrometer".
      71             :                 Default None expects a `pint.Quantity`.
      72             :             static_limit: If True, the static limit is used.
      73             :                 Default True.
      74             :             interaction_order: The order of interaction, e.g., 3 for dipole-dipole.
      75             :                 Defaults to 3.
      76             :             without_vacuum_contribution: If True, return only the scattered contribution and
      77             :                 omit the homogeneous vacuum term. Defaults to False.
      78             : 
      79             :         """
      80           1 :         super().__init__(
      81             :             pos1, pos2, unit, static_limit, interaction_order, without_vacuum_contribution=without_vacuum_contribution
      82             :         )
      83           1 :         self.point_on_plane_au = np.array(
      84             :             [QuantityScalar.convert_user_to_au(v, unit, "distance") for v in point_on_plane]
      85             :         )
      86           1 :         if np.isclose(np.linalg.norm(surface_normal), 0):
      87           1 :             raise ValueError("Normal vector cannot be zero.")
      88           1 :         self.surface_normal = np.array(surface_normal, dtype=float)
      89             :         # Almost perfect mirror # TODO make utils be able to handle inf
      90           1 :         self.surface_epsilon: PermittivityLike = 1e9
      91             : 
      92           1 :     def set_relative_permittivities(self, epsilon: PermittivityLike, surface_epsilon: PermittivityLike) -> Self:
      93             :         """Set the relative permittivities of the system.
      94             : 
      95             :         Args:
      96             :             epsilon: The relative permittivity (dimensionless) of the medium inside the cavity.
      97             :             surface_epsilon: The relative permittivity (dimensionless) of the surface.
      98             : 
      99             : 
     100             :         """
     101           0 :         if self.without_vacuum_contribution and epsilon != 1.0:  # NOSONAR
     102           0 :             raise ValueError(
     103             :                 "If the Green tensor is provided without the vacuum contribution, "
     104             :                 "the relative permittivity of the medium needs to be set to 1."
     105             :             )
     106             : 
     107           0 :         self.epsilon = epsilon
     108           0 :         self.surface_epsilon = surface_epsilon
     109           0 :         return self
     110             : 
     111           1 :     @override
     112           1 :     def _get_scaled_au(self, kappa1: int, kappa2: int, transition_energy_au: float) -> NDArray:
     113           1 :         if kappa1 == 1 and kappa2 == 1:
     114           1 :             return self._get_scaled_dipole_dipole_au(transition_energy_au)
     115           0 :         raise NotImplementedError("Only dipole-dipole Green tensors are currently implemented.")
     116             : 
     117           1 :     def _get_scaled_dipole_dipole_au(self, transition_energy_au: float) -> NDArray:
     118             :         """Calculate the dipole dipole Green tensor in cartesian coordinates for a single surface in atomic units.
     119             : 
     120             :         Args:
     121             :             transition_energy_au: The transition energy in atomic units at which to evaluate the Green tensor.
     122             : 
     123             :         Returns:
     124             :             The dipole dipole Green tensor in cartesian coordinates as a 3x3 array in atomic units (i.e. 1/bohr).
     125             : 
     126             :         """
     127           1 :         au_to_meter: float = ureg.Quantity(1, "atomic_unit_of_length").to("meter").magnitude
     128           1 :         lab_to_local_rotation = get_lab_to_local_rotation_matrix(self.surface_normal)
     129             : 
     130           1 :         pos1_local_m = rotate_vector_to_local(self.pos1_au * au_to_meter, lab_to_local_rotation)
     131           1 :         pos2_local_m = rotate_vector_to_local(self.pos2_au * au_to_meter, lab_to_local_rotation)
     132           1 :         point_on_plane_local_m = rotate_vector_to_local(self.point_on_plane_au * au_to_meter, lab_to_local_rotation)
     133             : 
     134           1 :         z1_m = point_on_plane_local_m[2]
     135             :         # Assume two surfaces, where the further apart atom is located in the center
     136             :         # but the second surface has the same permittivity as the inbetween medium
     137           1 :         height = 2 * max(abs(pos1_local_m[2] - z1_m), abs(pos2_local_m[2] - z1_m))
     138           1 :         if pos1_local_m[2] < z1_m and pos2_local_m[2] < z1_m:
     139           1 :             z2_m = z1_m - height
     140           1 :         elif pos1_local_m[2] > z1_m and pos2_local_m[2] > z1_m:
     141           1 :             z2_m = z1_m + height
     142             :         else:
     143           0 :             raise ValueError("Both atoms must be located either above or below the surface.")
     144             : 
     145           1 :         omega_hz = ureg.Quantity(transition_energy_au, "hartree").to("hbar Hz").magnitude
     146           1 :         epsilon = evaluate_relative_permittivity(self.epsilon, transition_energy_au, "hartree")
     147           1 :         epsilon1 = evaluate_relative_permittivity(self.surface_epsilon, transition_energy_au, "hartree")
     148           1 :         epsilon2 = epsilon
     149             : 
     150             :         # unit: # m^(-3) [hbar]^(-1) [epsilon_0]^(-1)
     151           1 :         gt = dynamic_green_tensor_scattered(
     152             :             pos1_local_m, pos2_local_m, z1_m, z2_m, omega_hz, epsilon, epsilon1, epsilon2, only_real_part=True
     153             :         )
     154           1 :         if not self.without_vacuum_contribution:
     155           1 :             gt += dynamic_green_tensor_homogeneous(pos1_local_m, pos2_local_m, omega_hz, epsilon, only_real_part=True)
     156           1 :         gt = rotate_tensor_to_lab(gt, lab_to_local_rotation)
     157           1 :         to_au = au_to_meter ** (-3) * ((4 * np.pi) ** (-1)) / (const.epsilon_0 * const.hbar)
     158             :         # hbar * epsilon_0 = (4*np.pi)**(-1) in atomic units
     159           1 :         return np.real(gt) / to_au