Line data    Source code

       1             : # SPDX-FileCopyrightText: 2024 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             : 
      10             : if TYPE_CHECKING:
      11             :     from .utils import PairinteractionModule
      12             : 
      13             : 
      14           1 : def test_diamagnetism(pi_module: PairinteractionModule) -> None:
      15             :     """Test calculating diamagnetic interactions."""
      16             :     # Create a basis
      17           1 :     ket = pi_module.KetAtom("Rb", n=60, l=0, j=0.5, m=0.5)
      18           1 :     basis = pi_module.BasisAtom("Rb", n=(ket.n - 2, ket.n + 2), l=(0, ket.l + 2))
      19           1 :     print(f"Number of basis states: {basis.number_of_states}")
      20             : 
      21             :     # Create system for a magnetic field of 1000 G
      22           1 :     bfield = 1000
      23           1 :     system = pi_module.SystemAtom(basis).set_magnetic_field([0, 0, bfield], unit="G").set_diamagnetism_enabled(True)
      24             : 
      25             :     # Diagonalize the system
      26           1 :     system = system.diagonalize(diagonalizer="eigen", sort_by_energy=True)
      27             : 
      28             :     # Get eigenenergies and the overlap with |ket>
      29           1 :     overlaps = system.basis.get_overlaps(ket)
      30           1 :     eigenenergies = system.get_eigenenergies(unit="GHz") - ket.get_energy(unit="GHz")
      31             : 
      32             :     # Get the overlap and eigenenergy corresponding to |ket>
      33           1 :     idx = np.argmax(overlaps)
      34           1 :     overlap = overlaps[idx]
      35           1 :     eigenenergy = eigenenergies[idx]
      36           1 :     print(
      37             :         f"The state |{ket}> in a field of {bfield} G has an energy of {eigenenergy:.3f} GHz "
      38             :         f"and overlap of {overlap:.2%} with the unperturbed state."
      39             :     )
      40             : 
      41             :     # Compare to reference data
      42           1 :     assert np.isclose(overlap, system.basis.get_corresponding_state(ket).get_overlap(ket))
      43           1 :     assert np.isclose(overlap, 0.9823116102876408, atol=1e-10)
      44           1 :     assert np.isclose(eigenenergy, 4.113262772909366)
      45             : 
      46             : 
      47           1 : def test_diamagnetism_angle_dependence() -> None:
      48             :     """Test calculating diamagnetic interactions for differently aligned magnetic fields."""
      49           1 :     import pairinteraction as pi_module  # only test with complex due to y component
      50             : 
      51             :     # Create a basis
      52           1 :     basis = pi_module.BasisAtom("Rb", n=(58, 62), l=(0, 2))
      53             : 
      54             :     # Create systems for fields in different directions
      55           1 :     bfield = 1000
      56           1 :     system_x = pi_module.SystemAtom(basis).set_magnetic_field([bfield, 0, 0], unit="G").set_diamagnetism_enabled(True)
      57           1 :     system_y = pi_module.SystemAtom(basis).set_magnetic_field([0, bfield, 0], unit="G").set_diamagnetism_enabled(True)
      58           1 :     system_z = pi_module.SystemAtom(basis).set_magnetic_field([0, 0, bfield], unit="G").set_diamagnetism_enabled(True)
      59             : 
      60             :     # Diagonalize the systems in parallel
      61           1 :     pi_module.diagonalize([system_x, system_y, system_z], diagonalizer="eigen", sort_by_energy=True)
      62             : 
      63             :     # Ensure that all eigenenergies are the same
      64           1 :     np.testing.assert_allclose(system_x.get_eigenenergies(unit="GHz"), system_y.get_eigenenergies(unit="GHz"))
      65           1 :     np.testing.assert_allclose(system_x.get_eigenenergies(unit="GHz"), system_z.get_eigenenergies(unit="GHz"))