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