90 lines
2.9 KiB
Python
90 lines
2.9 KiB
Python
|
#!/usr/bin/env python3
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
from scipy.special import gamma
|
||
|
|
|
||
|
|
# from sympy.physics.quantum.cg import CG
|
||
|
|
# from sympy import S
|
||
|
|
# def clebsch_gordan(j1, m1, j2, m2, j, m):
|
||
|
|
# cg = CG(S(j1), S(m1), S(j2), S(m2), S(j), S(m))
|
||
|
|
# result = cg.doit()
|
||
|
|
# return np.complex128(result)
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
from math import sqrt
|
||
|
|
|
||
|
|
def quantum_factorial(n):
|
||
|
|
"""
|
||
|
|
Calculate factorial for integer or half-integer numbers using gamma function.
|
||
|
|
For integer n: n! = n * (n-1) * ... * 1
|
||
|
|
For half-integer n: n! = Γ(n + 1)
|
||
|
|
"""
|
||
|
|
if n < 0:
|
||
|
|
return 0.0
|
||
|
|
return gamma(n + 1)
|
||
|
|
|
||
|
|
def clebsch_gordan(j1, m1,j2, m2, j, m):
|
||
|
|
"""
|
||
|
|
Calculate Clebsch-Gordan coefficient <j1 m1 j2 m2 | j m>
|
||
|
|
|
||
|
|
Parameters:
|
||
|
|
j1, j2: angular momentum quantum numbers
|
||
|
|
m1, m2: magnetic quantum numbers
|
||
|
|
j: total angular momentum quantum number
|
||
|
|
m: total magnetic quantum number
|
||
|
|
|
||
|
|
Returns:
|
||
|
|
float: Clebsch-Gordan coefficient value
|
||
|
|
"""
|
||
|
|
# Check validity of inputs using triangular inequalities and conservation
|
||
|
|
if not np.isclose(m, m1 + m2, atol=1e-10):
|
||
|
|
return 0.0
|
||
|
|
if abs(m1) > j1 or abs(m2) > j2 or abs(m) > j:
|
||
|
|
return 0.0
|
||
|
|
if not (abs(j1 - j2) <= j <= j1 + j2):
|
||
|
|
return 0.0
|
||
|
|
if j1 < 0 or j2 < 0 or j < 0:
|
||
|
|
return 0.0
|
||
|
|
|
||
|
|
# Ensure all quantum numbers are either integer or half-integer
|
||
|
|
if not (np.mod(2*j1, 1) < 1e-10 or np.isclose(np.mod(2*j1, 1), 1, atol=1e-10)):
|
||
|
|
return 0.0
|
||
|
|
if not (np.mod(2*j2, 1) < 1e-10 or np.isclose(np.mod(2*j2, 1), 1, atol=1e-10)):
|
||
|
|
return 0.0
|
||
|
|
if not (np.mod(2*j, 1) < 1e-10 or np.isclose(np.mod(2*j, 1), 1, atol=1e-10)):
|
||
|
|
return 0.0
|
||
|
|
|
||
|
|
# Calculate the coefficient
|
||
|
|
prefactor = sqrt((2*j + 1) * quantum_factorial(j1 + j2 - j) *
|
||
|
|
quantum_factorial(j1 - j2 + j) *
|
||
|
|
quantum_factorial(-j1 + j2 + j) /
|
||
|
|
quantum_factorial(j1 + j2 + j + 1))
|
||
|
|
|
||
|
|
prefactor *= sqrt(quantum_factorial(j + m) * quantum_factorial(j - m) *
|
||
|
|
quantum_factorial(j1 - m1) * quantum_factorial(j1 + m1) *
|
||
|
|
quantum_factorial(j2 - m2) * quantum_factorial(j2 + m2))
|
||
|
|
|
||
|
|
# Sum over k
|
||
|
|
sum_result = 0.0
|
||
|
|
k_min = max(0, max(j2 - j - m1, j1 + m2 - j))
|
||
|
|
k_max = min(j1 + j2 - j, min(j1 - m1, j2 + m2))
|
||
|
|
|
||
|
|
# Ensure k_min and k_max are integers for the range
|
||
|
|
k_min = int(np.ceil(k_min))
|
||
|
|
k_max = int(np.floor(k_max))
|
||
|
|
|
||
|
|
for k in range(k_min, k_max + 1):
|
||
|
|
denominator = (quantum_factorial(k) *
|
||
|
|
quantum_factorial(j1 + j2 - j - k) *
|
||
|
|
quantum_factorial(j1 - m1 - k) *
|
||
|
|
quantum_factorial(j2 + m2 - k) *
|
||
|
|
quantum_factorial(j - j2 + m1 + k) *
|
||
|
|
quantum_factorial(j - j1 - m2 + k))
|
||
|
|
if np.isclose(denominator, 0, atol=1e-10):
|
||
|
|
continue
|
||
|
|
term = (-1)**k / denominator
|
||
|
|
sum_result += term
|
||
|
|
|
||
|
|
return prefactor * sum_result
|
||
|
|
|