ALLOSBRAC doesn't work

ho_basis.jl

@@ -2,6 +2,7 @@ using SparseArrays
 using NuclearToolkit
 using SpecialFunctions
 include("helper.jl")
+include("osbrackets.jl")
 
 # Gaussian potentials in HO space
 inv_factorial(n) = Iterators.prod(inv.(1:n))
@@ -36,7 +37,7 @@ function get_2p_basis(E_max)
     l2s = Int[]
 
     # E = 2*n1 + l1 + 2*n2 + l2
-    for E in 0 : 2*E_max
+    for E in 2*E_max : -2 : 0
         for n1 in 0 : E ÷ 2
             for n2 in 0 : (E - 2*n1) ÷ 2
                 for l1 in 0 : (E - 2*n1 - 2*n2)
@@ -84,14 +85,8 @@ function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float
 end
 
 function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
-    E_max = maximum(Es)
-    j_max = 2 * E_max + 1
-    l_max = j_max
-    to = 0 # unused
-    
-    dtri = NuclearToolkit.prep_dtri(l_max + 1);
-    dcgm0 = NuclearToolkit.prep_dcgm0(l_max);
-    d6j = nothing # NuclearToolkit.prep_d6j_int(E_max, j_max, to);
+    Emax = maximum(Es) ÷ 2      # TODO: Too many steps. Simplify.
+    BRAC = cal_BRAC(Emax, Λ)
 
     mat = spzeros(length(Es), length(Es))
     s = hcat(Es, n1s, l1s, n2s, l2s)
@@ -100,7 +95,8 @@ function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
         (Elhs, N, L, n, l) = s[i, :]
         (Erhs, n1, l1, n2, l2) = s[j, :]
         if Elhs == Erhs
-            mat[i, j] = NuclearToolkit.HObracket_d6j(N, L, n, l, n1, l1, n2, l2, Λ, 1.0, dtri, dcgm0, d6j)
+            phase = (-1)^(N + n + n1 + n2)
+            mat[i, j] = phase * get_bracket(BRAC, Emax, Λ, N, L, n, l, n1, l1, n2, l2)
         end
     end
 

osbrackets.jl

@@ -0,0 +1,34 @@
+# compile using
+# gfortran -shared -fPIC osbrac.f90 -o osbrac.so
+# gfortran -shared -fPIC allosbrac.f90 -o allosbrac.so
+
+function cal_BRAC(Emax, Λ)
+    LMIN = Λ
+    LMAX = Λ
+    NQMAX = 2 * Emax
+    CO = 1/sqrt(2)
+    SI = 1/sqrt(2)
+    # dimensions BRAC(0:LMAX,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:(NQMAX-LMIN)/2,0:LMAX,0:(NQMAX-LMIN)/2,LMIN:LMAX)
+    BRAC = zeros(Float64, 1 + LMAX, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + (NQMAX - LMIN) ÷ 2, 1 + LMAX, 1 + (NQMAX-LMIN) ÷ 2, 1)
+    @ccall "OSBRACKETS/allosbrac.so".allosbrac_(NQMAX::Ref{Int32},LMIN::Ref{Int32},LMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},BRAC::Ptr{Array{Float64}})::Cvoid
+
+    return BRAC
+end
+
+function get_bracket(BRAC, Emax, Λ, n1′, l1′, n2′, l2′, n1, l1, n2, l2)
+    # BRAC(NP,N1P,MP,N1,N2,N,M,L)
+
+    Nq = l1 + l2 + 2 * (n1 + n2)
+    Nq′ = l1′ + l2′ + 2 * (n1′ + n2′)
+    if Nq ≠ Nq′; return 0; end
+
+    ϵ = (2 * Emax - Λ) % 2
+    NP = (l1′ - l2′ + Λ - ϵ) ÷ 2
+    N1P = n1′
+    MP = (l1′ + l2′ - Λ - ϵ) ÷ 2
+    N1 = n1
+    N2 = n2
+    N = (l1 - l2 + Λ - ϵ) ÷ 2
+    M = (l1 + l2 - Λ - ϵ) ÷ 2
+    return BRAC[1 + NP, 1 + N1P, 1 + MP, 1 + N1, 1 + N2, 1 + N, 1 + M, 1]
+end