Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

HalfInt usage broken #7

Open
Lala5th opened this issue Oct 30, 2024 · 0 comments
Open

HalfInt usage broken #7

Lala5th opened this issue Oct 30, 2024 · 0 comments

Comments

@Lala5th
Copy link

Lala5th commented Oct 30, 2024

The examples as well as code containing half-integer moments fail. This can be reproduced with a fresh install. An example trace is pasted below. I tracked this to a change in sympy and could circumvent by reverting the sympy version to 1.12. One solution might be to use Rational from sympy for the half integers moving forward.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 10
      7 mol.Nmax = 2
      9 # Generate Hamiltonians
---> 10 H0 = diatomic.operators.hyperfine_ham(mol)
     11 Hz = diatomic.operators.zeeman_ham(mol)
     13 # Parameter Space

File ~\.conda\envs\py310\lib\site-packages\diatomic\__init__.py:35, in log_time.<locals>.log_time_wrapper(*args, **kwargs)
     33 logger.info(f"Starting Function {func.__name__}...")
     34 start_time = time.perf_counter()
---> 35 result = func(*args, **kwargs)
     36 end_time = time.perf_counter()
     37 total_time = end_time - start_time

File ~\.conda\envs\py310\lib\site-packages\diatomic\operators.py:942, in hyperfine_ham(mol)
    936 scalar_part = (
    937     scalar_nuclear(mol.Ci[0], N, I1)
    938     + scalar_nuclear(mol.Ci[1], N, I2)
    939     + scalar_nuclear(mol.C4, I1, I2)
    940 )
    941 tensor_nuclear_part = tensor_nuclear(mol.C3, I1, I2, mol.Ii[0], mol.Ii[1], mol.Nmax)
--> 942 quadrupole_part = quadrupole(mol)
    943 H = rotational_part + scalar_part + tensor_nuclear_part + quadrupole_part
    944 return H

File ~\.conda\envs\py310\lib\site-packages\diatomic\operators.py:755, in quadrupole(mol)
    739 """Calculates the nuclear electric quadrupole interaction energy.
    740
    741 Computes the quadrupole interaction terms for the full hyperfine Hamiltonian using
   (...)
    751     np.ndarray: The quadrupole interaction term of the hyperfine Hamiltonian.
    752 """
    754 TdE = expanded_electric_gradient(mol)
--> 755 Tq0 = expanded_quad_moment(mol, 0)
    756 Tq1 = expanded_quad_moment(mol, 1)
    758 if mol.Ii[0] < 1:

File ~\.conda\envs\py310\lib\site-packages\diatomic\operators.py:715, in expanded_quad_moment(mol, nucleus)
    704 def expanded_quad_moment(mol, nucleus):
    705     """Expands the quadrupole moment tensor into the full hyperfine basis.
    706
    707     Args:
   (...)
    713         np.ndarray: The expanded nuclear quadrupole moment tensor.
    714     """
--> 715     T_nucleus = quad_moment(mol.Ii[nucleus])
    717     # Expand into full hyperfine basis
    718     num_N_proj = num_proj_with_below(mol.Nmax)

File ~\.conda\envs\py310\lib\site-packages\diatomic\operators.py:697, in quad_moment(I_nuc)
    693     for j, M2 in enumerate(proj_iter(I_nuc)):
    694         for n, q in enumerate(range(-2, 2 + 1)):
    695             T[n][i, j] = (
    696                 (-1) ** int(I_nuc - M1)
--> 697                 * wigner_3j(I_nuc, 2, I_nuc, -M1, q, M2)
    698                 / wigner_3j(I_nuc, 2, I_nuc, -I_nuc, 0, I_nuc)
    699             )
    701 return T

File ~\.conda\envs\py310\lib\site-packages\sympy\physics\wigner.py:219, in wigner_3j(j_1, j_2, j_3, m_1, m_2, m_3)
    130 def wigner_3j(j_1, j_2, j_3, m_1, m_2, m_3):
    131     r"""
    132     Calculate the Wigner 3j symbol `\operatorname{Wigner3j}(j_1,j_2,j_3,m_1,m_2,m_3)`.
    133
   (...)
    216     - Jens Rasch (2009-03-24): initial version
    217     """
--> 219     j_1, j_2, j_3, m_1, m_2, m_3 = map(_int_or_halfint,
    220                                        [j_1, j_2, j_3, m_1, m_2, m_3])
    222     if m_1 + m_2 + m_3 != 0:
    223         return S.Zero

File ~\.conda\envs\py310\lib\site-packages\sympy\physics\wigner.py:127, in _int_or_halfint(value)
    125 elif isinstance(value, Float):
    126     return _int_or_halfint(float(value))
--> 127 raise ValueError("expecting integer or half-integer, got %s" % value)

ValueError: expecting integer or half-integer, got (5/2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant