diff --git a/ho_basis.jl b/ho_basis.jl
index 40ffd11..4c9bb8a 100644
--- a/ho_basis.jl
+++ b/ho_basis.jl
@@ -1,5 +1,4 @@
 using SparseArrays
-using NuclearToolkit
 using SpecialFunctions
 include("helper.jl")
 include("osbrackets.jl")
@@ -86,7 +85,15 @@ end
 
 function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
     Emax = maximum(Es) ÷ 2      # TODO: Too many steps. Simplify.
-    BRAC = cal_BRAC(Emax, Λ)
+
+    ul1s = unique(l1s)
+    ul2s = unique(l2s)
+    BRACs = Matrix{Array}(undef, length(ul1s), length(ul2s))
+    for l1 in ul1s
+        for l2 in ul2s
+            BRACs[l1, l2] = cal_BRAC(Emax, Λ, l1, l2)
+        end
+    end
 
     mat = spzeros(length(Es), length(Es))
     s = hcat(Es, n1s, l1s, n2s, l2s)
@@ -96,7 +103,7 @@ function Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ)
         (Erhs, n1, l1, n2, l2) = s[j, :]
         if Elhs == Erhs
             phase = (-1)^(N + n + n1 + n2)
-            mat[i, j] = phase * get_bracket(BRAC, Emax, Λ, N, L, n, l, n1, l1, n2, l2)
+            mat[i, j] = phase * get_bracket(BRACs[l1, l2], Emax, Λ, N, L, n, l, n1, n2)
         end
     end
 
diff --git a/osbrackets.jl b/osbrackets.jl
index f64ca23..1a02122 100644
--- a/osbrackets.jl
+++ b/osbrackets.jl
@@ -2,33 +2,41 @@
 # gfortran -shared -fPIC osbrac.f90 -o osbrac.so
 # gfortran -shared -fPIC allosbrac.f90 -o allosbrac.so
 
-function cal_BRAC(Emax, Λ)
-    LMIN = Λ
-    LMAX = Λ
+using Libdl
+
+function cal_BRAC(Emax, Λ, l1, l2)
+    ϵ = (2 * Emax - Λ) % 2
+    N = (l1 - l2 + Λ - ϵ) ÷ 2
+    M = (l1 + l2 - Λ - ϵ) ÷ 2
+    L = Λ
     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
+    FIRSTCALL = true
+    # dimensions BRAC(0:L,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2,0:(NQMAX-L)/2)
+    BRAC = zeros(Float64, 1 + L, 1 + (NQMAX - L) ÷ 2, 1 + (NQMAX - L) ÷ 2, 1 + (NQMAX - L) ÷ 2, 1 + (NQMAX - L) ÷ 2)
+
+    lib = Libdl.dlopen("OSBRACKETS/osbrac.so") # Open the library explicitly.
+    sym = Libdl.dlsym(lib, :osbrac_)   # Get a symbol for the function to call.
+    # call signature OSBRAC(N,M,L,NQMAX,CO,SI,FIRSTCALL,BRAC)
+    @ccall $sym(N::Ref{Int32},M::Ref{Int32},L::Ref{Int32},NQMAX::Ref{Int32},CO::Ref{Float64},SI::Ref{Float64},FIRSTCALL::Ref{UInt8},BRAC::Ptr{Array{Float64}})::Cvoid
+    Libd.dlclose(lib)
 
     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
-
+function get_bracket(BRAC, Emax, Λ, n1′, l1′, n2′, l2′, n1, n2)
     ϵ = (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]
+    N2P = n2′
+    if MP+N1P+N2P == M+N1+N2
+        # BRAC(NP,N1P,MP,N1,N2)
+        return BRAC[1 + NP, 1 + N1P, 1 + MP, 1 + N1, 1 + N2]
+    else
+        return 0
+    end
 end