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)
      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)
      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)