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