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 pytest
      10              : 
      11              : if TYPE_CHECKING:
      12              :     from .utils import PairinteractionModule
      13              : 
      14              : 
      15            1 : @pytest.mark.parametrize("species", ["Yb171_mqdt", "Rb"])
      16            1 : def test_pair_potential(pi_module: PairinteractionModule, species: str) -> None:
      17              :     """Test multipole interaction."""
      18            1 :     ket = pi_module.KetAtom(species, nu=55.7, l=0, m=0.5)
      19              : 
      20              :     # Create a single-atom system
      21            1 :     basis = pi_module.BasisAtom(ket.species, nu=(ket.nu - 2, ket.nu + 2), l=(0, 3))
      22            1 :     print(f"Number of single-atom basis states: {basis.number_of_states}")
      23              : 
      24            1 :     system = pi_module.SystemAtom(basis)
      25              : 
      26              :     # Create two-atom systems for different interatomic distances and multipole orders
      27            1 :     delta_energy = 3  # GHz
      28            1 :     min_energy = 2 * ket.get_energy(unit="GHz") - delta_energy
      29            1 :     max_energy = 2 * ket.get_energy(unit="GHz") + delta_energy
      30              : 
      31            1 :     basis_pair = pi_module.BasisPair([system, system], energy=(min_energy, max_energy), energy_unit="GHz", m=(1, 1))
      32            1 :     print(f"Number of two-atom basis states: {basis_pair.number_of_states}")
      33              : 
      34            1 :     distances = np.linspace(0.2, 5, 5)
      35            1 :     system_pairs_0 = [pi_module.SystemPair(basis_pair) for d in distances]
      36            1 :     system_pairs_3 = [
      37              :         pi_module.SystemPair(basis_pair).set_interaction_order(3).set_distance(d, unit="micrometer") for d in distances
      38              :     ]
      39            1 :     system_pairs_4 = [
      40              :         pi_module.SystemPair(basis_pair).set_interaction_order(4).set_distance(d, unit="micrometer") for d in distances
      41              :     ]
      42              : 
      43              :     # Separate the contributions of the different multipole orders
      44            1 :     order_3 = [
      45              :         a.get_hamiltonian(unit="GHz").todense() - b.get_hamiltonian(unit="GHz").todense()
      46              :         for a, b in zip(system_pairs_3, system_pairs_0, strict=True)
      47              :     ]
      48            1 :     order_4 = [
      49              :         a.get_hamiltonian(unit="GHz").todense() - b.get_hamiltonian(unit="GHz").todense()
      50              :         for a, b in zip(system_pairs_4, system_pairs_3, strict=True)
      51              :     ]
      52              : 
      53              :     # Check that each order of the multipole expansion of the interaction has a significant contribution
      54              :     # at short distance
      55            1 :     norm_3 = np.linalg.norm(order_3, axis=(1, 2))
      56            1 :     norm_4 = np.linalg.norm(order_4, axis=(1, 2))
      57            1 :     assert norm_3[0] * distances[0] ** 3 > 1
      58            1 :     assert norm_4[0] * distances[0] ** 4 > 1
      59              : 
      60              :     # Check that for large/small distances, dipole-dipole/dipole-quadrupole interaction dominates
      61            1 :     assert norm_3[0] < norm_4[0]
      62            1 :     assert norm_3[-1] > norm_4[-1]
      63              : 
      64              :     # Check that each order of the multipole expansion scales as expected
      65            1 :     assert np.allclose(norm_3 * distances**3, norm_3[0] * distances[0] ** 3)
      66            1 :     assert np.allclose(norm_4 * distances**4, norm_4[0] * distances[0] ** 4)