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