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