test cg and RI functions

Project.toml

@@ -4,20 +4,24 @@ authors = ["Christoph Ortner <christophortner@gmail.com> and contributors"]
 version = "0.0.3"
 
 [deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+PartialWaveFunctions = "793d2195-304b-438e-bbb1-bc33c872ac39"
+Polynomials4ML = "03c4bcba-a943-47e9-bfa1-b1661fc2974f"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SpheriCart = "5caf2b29-02d9-47a3-9434-5931c85ba645"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-julia = "1"
 StaticArrays = "1.5"
+julia = "1"
 
 [extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
 Polynomials4ML = "03c4bcba-a943-47e9-bfa1-b1661fc2974f"
-WignerD = "87c4ff3e-34df-11e9-37a7-516cea4e0402"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
-BlockDiagonals = "0a1fb500-61f7-11e9-3c65-f5ef3456f9f0"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+WignerD = "87c4ff3e-34df-11e9-37a7-516cea4e0402"
 
 [targets]
-test = ["Test", "Polynomials4ML", "WignerD", "Rotations","BlockDiagonals"]
+test = ["Test", "Polynomials4ML", "WignerD", "Rotations", "BlockDiagonals"]

src/O3_alternative.jl

@@ -0,0 +1,92 @@
+# Alternative to the computation of rotation equivariant coupling coefficients
+
+using PartialWaveFunctions
+using Combinatorics
+using LinearAlgebra
+
+export re_basis_new
+
+function CG(l,m,L,N) 
+    M=m[1]+m[2]
+    if L[2]<abs(M)
+        return 0.
+    else
+        C=PartialWaveFunctions.clebschgordan(l[1],m[1],l[2],m[2],L[2],M) 
+    end
+    for k in 3:N
+        if L[k]<abs(M+m[k])
+            return 0.
+        elseif L[k-1]<abs(M)
+            return 0.
+        else
+            C*=PartialWaveFunctions.clebschgordan(L[k-1],M,l[k],m[k],L[k],M+m[k])
+            M+=m[k]
+        end
+    end
+    return C
+end
+
+function SetLl0(l,N)
+    set = Vector{Int64}[]
+    for k in abs(l[1]-l[2]):l[1]+l[2]
+        push!(set, [0; k])
+    end
+    for k in 3:N-1
+        setL=set
+        set=Vector{Int64}[]
+        for a in setL
+            for b in abs(a[k-1]-l[k]):a[k-1]+l[k]
+                push!(set, [a; b])
+            end
+        end
+    end
+    setL=set
+    set=Vector{Int64}[]
+    for a in setL
+        if a[N-1]==l[N]
+            push!(set, [a; 0])
+        end
+    end
+    return set
+end
+
+# Function that computes the set ML0
+function ML0(l,N)
+    setML = [[i] for i in -abs(l[1]):abs(l[1])]
+    for k in 2:N-1
+        set = setML
+        setML = Vector{Int64}[]
+        for m in set
+            append!(setML, [m; lk] for lk in -abs(l[k]):abs(l[k]) )
+        end
+    end
+    setML0=Vector{Int64}[]
+    for m in setML
+        s=sum(m)
+        if  abs(s) < abs(l[N])+1
+            push!(setML0, [m; -s])
+        end
+    end
+    return setML0
+end
+
+function re_basis_new(l)
+    N=size(l,1)
+    L=SetLl0(l,N)
+    r=size(L,1)
+    if r==0 
+        return zeros(Float64, 0, 0)
+    else 
+        setML0=ML0(l,N)
+        sizeML0=length(setML0)
+        U=zeros(Float64, r, sizeML0)
+        M = Vector{Int64}[]
+        for (j,m) in enumerate(setML0)
+            push!(M,m)
+            for i in 1:r
+                U[i,j]=CG(l,m,L[i],N)
+            end
+        end
+    end
+    return U,M
+end

src/RepLieGroups.jl

@@ -4,4 +4,6 @@ include("yyvector.jl")
 
 include("obsolete/obsolete.jl")
 
+include("O3_alternative.jl")
+
 end

src/obsolete/O3.jl

@@ -250,6 +250,23 @@ function Ctran(l::Int64,m::Int64,μ::Int64)
    end
 end
 
+
+# function Ctran(l::Int64,m::Int64,μ::Int64)
+# 	if abs(m) ≠ abs(μ)
+# 	   return 0
+# 	elseif abs(m) == 0
+# 	   return 1
+# 	elseif m < 0 && μ < 0
+# 	   return - im * (-1)^m # 1/sqrt(2)
+# 	elseif m < 0 && μ > 0
+# 	   return im # (-1)^m/sqrt(2)
+# 	elseif m > 0 && μ < 0
+# 	   return (-1)^m # - im * (-1)^m/sqrt(2)
+# 	else
+# 	   return 1. # im/sqrt(2)
+# 	end
+#  end
+
 Ctran(l::Int64) = sparse(Matrix{ComplexF64}([ Ctran(l,m,μ) for m = -l:l, μ = -l:l ])) |> dropzeros
 
 ## NOTE: Ctran(L) is the transformation matrix from rSH to cSH. More specifically, 

test/runtests.jl

@@ -4,6 +4,7 @@ using Test
 @testset "RepLieGroups.jl" begin
     # Write your tests here.
     @testset "SYYVector" begin include("test_yyvector.jl"); end
-    @testset "O3-ClebschGordan" begin include("test_obsolete_cg.jl"); end
-    @testset "O3" begin include("test_obsolete_o3.jl"); end
+    # @testset "O3-ClebschGordan" begin include("test_obsolete_cg.jl"); end
+    # @testset "O3" begin include("test_obsolete_o3.jl"); end
+    @testset "CGcoef vs PartialWaveFunctions" begin include("test_cg_vs_partialwavefunctions.jl"); end
 end

test/test_RI_basis.jl

@@ -0,0 +1,125 @@
+using SpheriCart, StaticArrays, LinearAlgebra, RepLieGroups, WignerD, 
+      Combinatorics
+using RepLieGroups.O3: Rot3DCoeffs, Rot3DCoeffs_real
+O3 = RepLieGroups.O3
+using Test
+
+# Evaluation of spherical harmonics
+function eval_cY(rbasis::SphericalHarmonics{LMAX}, 𝐫) where {LMAX}  
+    Yr = rbasis(𝐫)
+    Yc = zeros(Complex{eltype(Yr)}, length(Yr))
+    for l = 0:LMAX
+       # m = 0 
+       i_l0 = SpheriCart.lm2idx(l, 0)
+       Yc[i_l0] = Yr[i_l0]
+       # m ≠ 0 
+       for m = 1:l 
+          i_lm⁺ = SpheriCart.lm2idx(l,  m)
+          i_lm⁻ = SpheriCart.lm2idx(l, -m)
+          Ylm⁺ = Yr[i_lm⁺]
+          Ylm⁻ = Yr[i_lm⁻]
+          Yc[i_lm⁺] = (-1)^m * (Ylm⁺ + im * Ylm⁻) / sqrt(2)
+          Yc[i_lm⁻] = (Ylm⁺ - im * Ylm⁻) / sqrt(2)
+       end
+    end 
+    return Yc
+ end
+
+ function rand_sphere() 
+    u = @SVector randn(3)
+    return u / norm(u) 
+ end
+
+ function rand_rot() 
+    K = @SMatrix randn(3,3)
+    return exp(K - K') 
+ end
+
+ function f(Rs, q; coeffs=coeffs, MM=MM, ll=ll) 
+    Lmax = maximum(ll)
+    real_basis = SphericalHarmonics(Lmax)
+    YY = []
+    for i in 1:length(ll)
+        push!(YY, eval_cY(real_basis, Rs[i]))
+    end
+    out = zero(eltype(YY[1]))
+    for (c, mm) in zip(coeffs[q, :], MM) 
+        ind = Int[]
+        for i in 1:length(ll)
+            push!(ind, SpheriCart.lm2idx(ll[i], mm[i]))
+        end
+        out += c * prod(YY[i][ind[i]] for i in 1:length(ll))
+    end
+    return out 
+end
+
+
+# for the moment the code with generalized CG only works with L=0
+L = 0
+cc = Rot3DCoeffs(L)
+ll = SA[2,2,2,3,3]
+
+# version with svd
+@time coeffs1, MM1 = O3.re_basis(cc, ll)
+nbas = size(coeffs1, 1)
+
+#version with gen CG coefficients
+@time coeffs2, MM2 = re_basis_new(ll)
+
+# simple test on size
+@test size(coeffs1) == size(coeffs2)
+@test size(MM1) == size(MM2)
+
+P1 = sortperm(MM1)
+P2 = sortperm(MM2)
+MMsorted1 = MM1[P1]
+MMsorted2 = MM2[P2]
+# check that same mm values
+@test MMsorted1 == MMsorted2
+
+coeffsp1 = coeffs1[:,P1]
+coeffsp2 = coeffs2[:,P2]
+
+# test that full rank
+@test rank(coeffsp1) == size(coeffsp1,1)
+@test rank(coeffsp2) == size(coeffsp2,1)
+
+# check that the coef span the same space - test fails
+@test nbas == rank([coeffsp; coeffsp2], rtol = 1e-12)
+
+
+Rs = [rand_sphere() for _ in 1:length(ll)]
+Q = rand_rot() 
+QRs = [Q*Rs[i] for i in 1:length(Rs)]
+fRs1 = [ f(Rs, q; coeffs=coeffs1, MM=MM1, ll=ll) for q = 1:nbas ]
+fRs1Q = [ f(QRs, q; coeffs=coeffs1, MM=MM1, ll=ll) for q = 1:nbas ]
+
+# check invariance (for now)
+@test norm(fRs1 .- fRs1Q) < 1e-12
+
+fRs2 = [ f(Rs, q; coeffs=coeffs2, MM=MM2, ll=ll) for q = 1:nbas ]
+fRs2Q = [ f(QRs, q; coeffs=coeffs2, MM=MM2, ll=ll) for q = 1:nbas ]
+
+# check invariance (for now)
+@test norm(fRs2 .- fRs2Q) < 1e-12
+
+# Test on batch
+ntest = 1000
+A1 = zeros(nbas, ntest)
+A2 = zeros(nbas, ntest)
+for i = 1:ntest 
+   Rs = [rand_sphere() for _ in 1:length(ll)]
+   for q = 1:nbas
+       fRs = f(Rs, q; coeffs = coeffs1, MM=MM1, ll=ll)
+       @assert abs.(imag(fRs)) < 1e-16
+       A1[q, i] = real(fRs)
+
+       fRs2 = f(Rs, q; coeffs = coeffs2, MM=MM2, ll=ll)
+       @assert abs.(imag(fRs2)) < 1e-16
+       A2[q, i] = real(fRs2)
+   end
+end
+
+# check that functions span same space
+rk = rank([A1;A2]; rtol = 1e-12)  
+@test rk == nbas