Line data    Source code

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