Line data    Source code

       1             : # SPDX-FileCopyrightText: 2024 Pairinteraction Developers
       2             : # SPDX-License-Identifier: LGPL-3.0-or-later
       3             : 
       4             : """Test the conversion from BaseUnits (=atomic units) used in the backend and the units input and output to the user."""
       5             : 
       6           1 : import numpy as np
       7             : 
       8           1 : import pairinteraction.real as pi
       9           1 : from pairinteraction.units import BaseUnits, QuantityScalar, ureg
      10             : 
      11             : 
      12           1 : def test_magnetic() -> None:
      13             :     """Test magnetic units."""
      14           1 :     ket = pi.KetAtom("Rb", n=60, l=0, m=0.5)
      15             : 
      16           1 :     mu = ket.get_matrix_element(ket, "magnetic_dipole", q=0)
      17           1 :     mu = mu.to("bohr_magneton")
      18           1 :     lande_factor = 2.002319304363
      19           1 :     assert np.isclose(mu.magnitude, -1 / 2 * lande_factor)
      20             : 
      21             :     # check magnetic field conversion is correct
      22           1 :     B_z = QuantityScalar.from_unit(1, "gauss", "magnetic_field")
      23           1 :     B_z_pint = ureg.Quantity(1, "gauss").to("T", "Gaussian")
      24           1 :     assert np.isclose(B_z.to_base_unit(), B_z_pint.to_base_units().magnitude)
      25             : 
      26             :     # such that mu * B_z is of dimension energy
      27           1 :     zeeman_energy = -mu * B_z_pint
      28           1 :     assert zeeman_energy.dimensionality == BaseUnits["energy"].dimensionality
      29             : 
      30             :     # check against constructed Hamiltonian
      31           1 :     basis = pi.BasisAtom("Rb", n=(1, 1), additional_kets=[ket])
      32           1 :     system = pi.SystemAtom(basis)
      33           1 :     system.set_diamagnetism_enabled(False).set_magnetic_field([0, 0, B_z_pint])
      34           1 :     zeeman_energy_from_hamiltonian = system.get_hamiltonian("MHz")[0, 0] - ket.get_energy("MHz")
      35           1 :     assert np.isclose(zeeman_energy_from_hamiltonian, zeeman_energy.to("MHz", "spectroscopy").magnitude)
      36             : 
      37             : 
      38           1 : def test_electric_dipole() -> None:
      39             :     """Test electric dipole units."""
      40           1 :     ket_a = pi.KetAtom("Rb", n=60, l=0, m=0.5)
      41           1 :     ket_b = pi.KetAtom("Rb", n=61, l=0, m=0.5)
      42           1 :     ket_c = pi.KetAtom("Rb", n=60, l=1, j=3 / 2, m=0.5)
      43             : 
      44           1 :     dipole_a_c = ket_a.get_matrix_element(ket_c, "electric_dipole", q=0)
      45           1 :     dipole_b_c = ket_b.get_matrix_element(ket_c, "electric_dipole", q=0)
      46             : 
      47           1 :     kappa = ureg.Quantity(1 / (4 * np.pi), "1 / epsilon_0")
      48           1 :     c3 = kappa * dipole_a_c * dipole_b_c
      49             : 
      50           1 :     GHz = ureg.Quantity(1, "GHz")
      51           1 :     c3 = c3 * GHz / GHz.to("J", "spectroscopy")
      52           1 :     c3 = c3.to("GHz micrometer^3")
      53             : 
      54           1 :     distance = ureg.Quantity(10, "micrometer")
      55           1 :     basis = pi.BasisAtom("Rb", additional_kets=[ket_a, ket_b, ket_c])
      56           1 :     system = pi.SystemAtom(basis)
      57           1 :     basis_pair = pi.BasisPair([system, system])
      58           1 :     system_pair = pi.SystemPair(basis_pair)
      59           1 :     system_pair.set_interaction_order(3)
      60           1 :     system_pair.set_distance(distance)
      61           1 :     system.get_hamiltonian()
      62             : 
      63           1 :     ket_ab_idx = np.argmax(basis_pair.get_overlaps([ket_a, ket_b]))
      64           1 :     ket_cc_idx = np.argmax(basis_pair.get_overlaps([ket_c, ket_c]))
      65           1 :     H = system_pair.get_hamiltonian("GHz") * distance.to("micrometer").magnitude ** 3  # GHz * micrometer^3
      66             : 
      67           1 :     assert np.isclose(-2 * c3.magnitude, H[ket_ab_idx, ket_cc_idx])
      68           1 :     assert np.isclose(-2 * c3.magnitude, 5.73507543166919)