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 GreenTensorCavity(GreenTensorBase):
      33             :     """Green tensor for two atoms in a cavity (between two infinite planar surfaces).
      34             : 
      35             :     Examples:
      36             :         >>> from pairinteraction.green_tensor import GreenTensorCavity
      37             :         >>> gt = GreenTensorCavity(
      38             :         ...     [0, 0, 0],
      39             :         ...     [10, 0, 0],
      40             :         ...     point_on_plane1=[0, 0, -5],
      41             :         ...     point_on_plane2=[0, 0, 5],
      42             :         ...     surface_normal=[0, 0, 1],
      43             :         ...     unit="micrometer",
      44             :         ... )
      45             :         >>> transition_energy = 2  # h * GHz
      46             :         >>> gt_dipole_dipole = gt.get(1, 1, transition_energy, "planck_constant * GHz")
      47             :         >>> print(f"{gt_dipole_dipole[0, 0]:.2f}")
      48             :         -3.84 / bohr
      49             : 
      50             :     """
      51             : 
      52           1 :     def __init__(
      53             :         self,
      54             :         pos1: ArrayLike | PintArrayLike,
      55             :         pos2: ArrayLike | PintArrayLike,
      56             :         point_on_plane1: ArrayLike | PintArrayLike,
      57             :         point_on_plane2: ArrayLike | PintArrayLike,
      58             :         surface_normal: ArrayLike,
      59             :         unit: str | None = None,
      60             :         static_limit: bool = True,
      61             :         interaction_order: int = 3,
      62             :         *,
      63             :         without_vacuum_contribution: bool = False,
      64             :     ) -> None:
      65             :         """Create a Green tensor for two atoms inside a planar cavity formed by two infinite surfaces.
      66             : 
      67             :         The two surfaces of the cavity are parallel infinite planes defined by two points
      68             :         and a shared surface normal vector.
      69             :         If not specified otherwise (see `set_relative_permittivities`), the surfaces are treated as perfect mirrors.
      70             : 
      71             :         Args:
      72             :             pos1: Position of the first atom in the given unit.
      73             :             pos2: Position of the second atom in the given unit.
      74             :             point_on_plane1: A point on the first cavity surface in the given unit.
      75             :             point_on_plane2: A point on the second cavity surface in the given unit.
      76             :             surface_normal: The shared surface normal vector of both cavity surfaces.
      77             :             unit: The unit of the distance, e.g. "micrometer".
      78             :                 Default None expects a `pint.Quantity`.
      79             :             static_limit: If True, the static limit is used.
      80             :                 Default True.
      81             :             interaction_order: The order of interaction, e.g., 3 for dipole-dipole.
      82             :                 Defaults to 3.
      83             :             without_vacuum_contribution: If True, return only the scattered contribution and
      84             :                 omit the homogeneous vacuum term. Defaults to False.
      85             : 
      86             :         """
      87           1 :         super().__init__(
      88             :             pos1, pos2, unit, static_limit, interaction_order, without_vacuum_contribution=without_vacuum_contribution
      89             :         )
      90           1 :         self.point_on_plane1_au = np.array(
      91             :             [QuantityScalar.convert_user_to_au(v, unit, "distance") for v in point_on_plane1]
      92             :         )
      93           1 :         self.point_on_plane2_au = np.array(
      94             :             [QuantityScalar.convert_user_to_au(v, unit, "distance") for v in point_on_plane2]
      95             :         )
      96           1 :         if np.isclose(np.linalg.norm(surface_normal), 0):
      97           0 :             raise ValueError("Normal vector cannot be zero.")
      98           1 :         self.surface_normal = np.array(surface_normal, dtype=float)
      99             :         # Almost perfect mirrors # TODO make utils be able to handle inf
     100           1 :         self.surface1_epsilon: PermittivityLike = 1e9
     101           1 :         self.surface2_epsilon: PermittivityLike = 1e9
     102             : 
     103           1 :     def set_relative_permittivities(
     104             :         self, epsilon: PermittivityLike, epsilon1: PermittivityLike, epsilon2: PermittivityLike
     105             :     ) -> Self:
     106             :         """Set the relative permittivities of the system.
     107             : 
     108             :         Args:
     109             :             epsilon: The relative permittivity (dimensionless) of the medium inside the cavity.
     110             :             epsilon1: The relative permittivity (dimensionless) of the first surface.
     111             :             epsilon2: The relative permittivity (dimensionless) of the second surface.
     112             : 
     113             : 
     114             :         """
     115           0 :         if self.without_vacuum_contribution and epsilon != 1.0:  # NOSONAR
     116           0 :             raise ValueError(
     117             :                 "If the Green tensor is provided without the vacuum contribution, "
     118             :                 "the relative permittivity of the medium needs to be set to 1."
     119             :             )
     120             : 
     121           0 :         self.epsilon = epsilon
     122           0 :         self.surface1_epsilon = epsilon1
     123           0 :         self.surface2_epsilon = epsilon2
     124           0 :         return self
     125             : 
     126           1 :     @override
     127           1 :     def _get_scaled_au(self, kappa1: int, kappa2: int, transition_energy_au: float) -> NDArray:
     128           1 :         if kappa1 == 1 and kappa2 == 1:
     129           1 :             return self._get_scaled_dipole_dipole_au(transition_energy_au)
     130           0 :         raise NotImplementedError("Only dipole-dipole Green tensors are currently implemented.")
     131             : 
     132           1 :     def _get_scaled_dipole_dipole_au(self, transition_energy_au: float) -> NDArray:
     133             :         """Calculate the dipole dipole Green tensor in cartesian coordinates for a cavity in atomic units.
     134             : 
     135             :         Args:
     136             :             transition_energy_au: The transition energy in atomic units at which to evaluate the Green tensor.
     137             : 
     138             :         Returns:
     139             :             The dipole dipole Green tensor in cartesian coordinates as a 3x3 array in atomic units (i.e. 1/bohr).
     140             : 
     141             :         """
     142           1 :         au_to_meter: float = ureg.Quantity(1, "atomic_unit_of_length").to("meter").magnitude
     143           1 :         lab_to_local_rotation = get_lab_to_local_rotation_matrix(self.surface_normal)
     144             : 
     145           1 :         pos1_local_m = rotate_vector_to_local(self.pos1_au * au_to_meter, lab_to_local_rotation)
     146           1 :         pos2_local_m = rotate_vector_to_local(self.pos2_au * au_to_meter, lab_to_local_rotation)
     147           1 :         point_on_plane1_local_m = rotate_vector_to_local(self.point_on_plane1_au * au_to_meter, lab_to_local_rotation)
     148           1 :         point_on_plane2_local_m = rotate_vector_to_local(self.point_on_plane2_au * au_to_meter, lab_to_local_rotation)
     149             : 
     150           1 :         z1_m = point_on_plane1_local_m[2]
     151           1 :         z2_m = point_on_plane2_local_m[2]
     152           1 :         if np.isclose(z1_m, z2_m, atol=1e-12):
     153           1 :             raise ValueError("The two cavity planes must be distinct.")
     154             : 
     155           1 :         omega_hz = ureg.Quantity(transition_energy_au, "hartree").to("hbar Hz").magnitude
     156           1 :         epsilon = evaluate_relative_permittivity(self.epsilon, transition_energy_au, "hartree")
     157           1 :         epsilon1 = evaluate_relative_permittivity(self.surface1_epsilon, transition_energy_au, "hartree")
     158           1 :         epsilon2 = evaluate_relative_permittivity(self.surface2_epsilon, transition_energy_au, "hartree")
     159             : 
     160             :         # unit: # m^(-3) [hbar]^(-1) [epsilon_0]^(-1)
     161           1 :         gt = dynamic_green_tensor_scattered(
     162             :             pos1_local_m, pos2_local_m, z1_m, z2_m, omega_hz, epsilon, epsilon1, epsilon2, only_real_part=True
     163             :         )
     164           1 :         if not self.without_vacuum_contribution:
     165           1 :             gt += dynamic_green_tensor_homogeneous(pos1_local_m, pos2_local_m, omega_hz, epsilon, only_real_part=True)
     166           1 :         gt = rotate_tensor_to_lab(gt, lab_to_local_rotation)
     167           1 :         to_au = au_to_meter ** (-3) * ((4 * np.pi) ** (-1)) / (const.epsilon_0 * const.hbar)
     168             :         # hbar * epsilon_0 = (4*np.pi)**(-1) in atomic units
     169           1 :         return np.real(gt) / to_au