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