Line data    Source code

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