add wigner9j symbols
#18
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,3 +1,4 @@ | ||
| Manifest.toml | ||
| *.jl.cov | ||
| *.jl.*.cov | ||
| *.jl.mem |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,25 +1,25 @@ | ||
| name = "WignerSymbols" | ||
| uuid = "9f57e263-0b3d-5e2e-b1be-24f2bb48858b" | ||
| authors = ["Jutho Haegeman"] | ||
| version = "2.1.0" | ||
|
|
||
| [deps] | ||
| HalfIntegers = "f0d1745a-41c9-11e9-1dd9-e5d34d218721" | ||
| LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" | ||
| Primes = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae" | ||
| RationalRoots = "308eb6b3-cc68-5ff3-9e97-c3c4da4fa681" | ||
|
|
||
| [compat] | ||
| HalfIntegers = "1" | ||
| LRUCache = "1.3" | ||
| Primes = "0.4 - 0.9" | ||
| RationalRoots = "0.1 - 0.9" | ||
| julia = "1.5" | ||
|
|
||
| [extras] | ||
| LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" | ||
| Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" | ||
| Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" | ||
|
|
||
| [targets] | ||
| test = ["Test", "LinearAlgebra", "Random"] |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,5 +1,5 @@ | ||
| module WignerSymbols | ||
| export δ, Δ, clebschgordan, wigner3j, wigner6j, wigner9j, racahV, racahW, HalfInteger | ||
|
|
||
| using HalfIntegers | ||
| using RationalRoots | ||
|
|
@@ -14,16 +14,21 @@ convert(BigInt, primefactorial(401)) # trigger compilation and generate some fix | |
|
|
||
| const Key3j = Tuple{UInt,UInt,UInt,Int,Int} | ||
| const Key6j = NTuple{6,UInt} | ||
| const Key9j = NTuple{9,HalfInteger} | ||
|
|
||
| const Wigner3j = LRU{Key3j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) | ||
| const Wigner6j = LRU{Key6j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) | ||
| const Wigner9j = LRU{Key9j,Tuple{Rational{BigInt},Rational{BigInt}}}(; maxsize = 10^6) | ||
|
|
||
| function set_buffer3j_size!(; maxsize) | ||
| resize!(Wigner3j; maxsize = maxsize) | ||
| end | ||
| function set_buffer6j_size!(; maxsize) | ||
| resize!(Wigner6j; maxsize = maxsize) | ||
| end | ||
| function set_buffer9j_size!(; maxsize) | ||
| resize!(Wigner9j; maxsize = maxsize) | ||
| end | ||
|
|
||
| # check integerness and correctness of (j,m) angular momentum | ||
| ϵ(j, m) = (abs(m) <= j && ishalfinteger(j) && isinteger(j-m) && isinteger(j+m)) | ||
|
|
@@ -242,6 +247,99 @@ function racahW(T::Type{<:Real}, j₁, j₂, J, j₃, J₁₂, J₂₃) | |
| end | ||
| end | ||
|
|
||
| """ | ||
| wigner9j(T::Type{<:Real} = RationalRoot{BigInt}, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) -> ::T | ||
|
|
||
| Compute the value of the Wigner-9j symbol | ||
| ⎧ j₁ j₂ j₃ ⎫ | ||
| ⎨ j₄ j₅ j₆ ⎬ | ||
| ⎩ j₇ j₈ j₉ ⎭ | ||
| as a type `T` real number. By default `T = RationalRoot{BigInt}`. | ||
|
|
||
| Returns `zero(T)` if any of the triangle conditions TODO | ||
| are not satisfied, but throws a `DomainError` if | ||
| the `jᵢ`s are not integer or halfinteger. | ||
| """ | ||
| wigner9j(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) = wigner9j(RRBig, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) | ||
| function wigner9j(T::Type{<:Real}, j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) | ||
| # check validity of `jᵢ`s | ||
| for jᵢ in (j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) | ||
| (ishalfinteger(jᵢ) && jᵢ >= zero(jᵢ)) || throw(DomainError("invalid jᵢ", jᵢ)) | ||
| end | ||
| return _wigner9j(T, HalfInteger.((j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉))...) | ||
| end | ||
|
|
||
| const _perms9j = [([1, 2, 3], [1, 2, 3]) ([1, 2, 3], [1, 3, 2]) ([1, 2, 3], [2, 1, 3]) ([1, 2, 3], [2, 3, 1]) ([1, 2, 3], [3, 1, 2]) ([1, 2, 3], [3, 2, 1]); | ||
| ([1, 3, 2], [1, 2, 3]) ([1, 3, 2], [1, 3, 2]) ([1, 3, 2], [2, 1, 3]) ([1, 3, 2], [2, 3, 1]) ([1, 3, 2], [3, 1, 2]) ([1, 3, 2], [3, 2, 1]); | ||
| ([2, 1, 3], [1, 2, 3]) ([2, 1, 3], [1, 3, 2]) ([2, 1, 3], [2, 1, 3]) ([2, 1, 3], [2, 3, 1]) ([2, 1, 3], [3, 1, 2]) ([2, 1, 3], [3, 2, 1]); | ||
| ([2, 3, 1], [1, 2, 3]) ([2, 3, 1], [1, 3, 2]) ([2, 3, 1], [2, 1, 3]) ([2, 3, 1], [2, 3, 1]) ([2, 3, 1], [3, 1, 2]) ([2, 3, 1], [3, 2, 1]); | ||
| ([3, 1, 2], [1, 2, 3]) ([3, 1, 2], [1, 3, 2]) ([3, 1, 2], [2, 1, 3]) ([3, 1, 2], [2, 3, 1]) ([3, 1, 2], [3, 1, 2]) ([3, 1, 2], [3, 2, 1]); | ||
| ([3, 2, 1], [1, 2, 3]) ([3, 2, 1], [1, 3, 2]) ([3, 2, 1], [2, 1, 3]) ([3, 2, 1], [2, 3, 1]) ([3, 2, 1], [3, 1, 2]) ([3, 2, 1], [3, 2, 1])] | ||
|
|
||
| const _signs9j = [1 -1 -1 1 1 -1; | ||
| -1 1 1 -1 -1 1; | ||
| -1 1 1 -1 -1 1; | ||
| 1 -1 -1 1 1 -1; | ||
| 1 -1 -1 1 1 -1; | ||
| -1 1 1 -1 -1 1] | ||
|
|
||
| function _wigner9j(T::Type{<:Real}, j₁::HalfInteger, j₂::HalfInteger, j₃::HalfInteger, | ||
| j₄::HalfInteger, j₅::HalfInteger, j₆::HalfInteger, | ||
| j₇::HalfInteger, j₈::HalfInteger, j₉::HalfInteger) | ||
|
|
||
| # check triangle conditions | ||
| if !(δ(j₁, j₂, j₃) && δ(j₄, j₅, j₆) && δ(j₇, j₈, j₉) && δ(j₁, j₄, j₇) && δ(j₂, j₅, j₈) && δ(j₃, j₆, j₉)) | ||
| return zero(T) | ||
| end | ||
|
|
||
| # dictionary lookup, check all 72 permutations | ||
| k = [j₁ j₂ j₃; j₄ j₅ j₆; j₇ j₈ j₉] | ||
| for (p, m) in zip(_perms9j, _signs9j) | ||
| kk = Tuple(reshape(k[p...], 9)) | ||
| kkT = Tuple(reshape(transpose(k[p...]), 9)) | ||
| if haskey(Wigner9j, kk) | ||
| r, s = Wigner9j[kk] | ||
| return m * _convert(T, s) * convert(T, signedroot(r)) | ||
| elseif haskey(Wigner9j, kkT) | ||
| r, s = Wigner9j[kkT] | ||
| return m * _convert(T, s) * convert(T, signedroot(r)) | ||
| end | ||
| end | ||
|
|
||
| # order irrelevant: product remains the same | ||
| n₁, d₁ = Δ²(j₁, j₂, j₃) | ||
| n₂, d₂ = Δ²(j₄, j₅, j₆) | ||
| n₃, d₃ = Δ²(j₇, j₈, j₉) | ||
| n₄, d₄ = Δ²(j₁, j₄, j₇) | ||
| n₅, d₅ = Δ²(j₂, j₅, j₈) | ||
| n₆, d₆ = Δ²(j₃, j₆, j₉) | ||
|
|
||
| snum, rnum = splitsquare(n₁ * n₂ * n₃ * n₄ * n₅ * n₆) | ||
| sden, rden = splitsquare(d₁ * d₂ * d₃ * d₄ * d₅ * d₆) | ||
|
|
||
| rnum, rden = divgcd!(rnum, rden) | ||
| rr = Base.unsafe_rational(convert(BigInt, rnum), convert(BigInt, rden)) | ||
|
|
||
| snum, sden = divgcd!(snum, sden) | ||
| snum2 = compute9jseries(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉) | ||
| for n = 1:length(sden.powers) | ||
| p = bigprime(n) | ||
| while sden.powers[n] > 0 | ||
| q, r = divrem(snum2, p) | ||
| if iszero(r) | ||
| snum2 = q | ||
| sden.powers[n] -= 1 | ||
| else | ||
| break | ||
| end | ||
| end | ||
| end | ||
| s = Base.unsafe_rational(convert(BigInt, snum)*snum2, convert(BigInt, sden)) | ||
|
|
||
| Wigner9j[(j₁, j₂, j₃, j₄, j₅, j₆, j₇, j₈, j₉)] = (rr, s) | ||
| return _convert(T, s) * convert(T, signedroot(rr)) | ||
| end | ||
|
|
||
| # COMPUTATIONAL ROUTINES | ||
| #------------------------ | ||
| # squared triangle coefficient | ||
|
|
@@ -312,7 +410,6 @@ function compute3jseries(β₁, β₂, β₃, α₁, α₂) | |
| end | ||
| den = commondenominator!(nums, dens) | ||
| totalnum = sumlist!(nums) | ||
tgorordo marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| for n = 1:length(den.powers) | ||
| p = bigprime(n) | ||
| while den.powers[n] > 0 | ||
|
|
@@ -365,13 +462,84 @@ function compute6jseries(β₁, β₂, β₃, α₁, α₂, α₃, α₄) | |
| return Base.unsafe_rational(totalnum, totalden) | ||
| end | ||
|
|
||
| # compute the sum appearing in the 9j symbol | ||
| function compute9jseries(a, b, c, d, e, f, g, h, j) | ||
| I1 = max(abs(h - d), abs(b - f), abs(a - j)) | ||
| I2 = min(h + d, b + f, a + j) | ||
| krange = I1:I2 | ||
|
|
||
| s = sum(krange) do k | ||
| p = iseven(2k) ? big(2k + 1) : -big(2k + 1) | ||
|
|
||
| b₁ = let (m₁, m₂, m₃, m₄, m₅, m₆) = (a, b, c, f, j, k) | ||
| α₁ = convert(BigInt, m₁ + m₅ - m₆) | ||
| α₂ = convert(BigInt, m₁ - m₅ + m₆) | ||
| α₃ = convert(BigInt, -m₁ + m₅ + m₆) | ||
| β₁ = convert(BigInt, m₁ + m₅ + m₆) | ||
| β₂ = convert(BigInt, m₂ + m₄ + m₆) | ||
| β₃ = convert(BigInt, m₃ + m₄ + m₅) | ||
| β₀ = convert(BigInt, m₁ + m₂ + m₃) | ||
| wei9jbracket(α₁, α₂, α₃, β₁, β₂, β₃, β₀) | ||
| end | ||
|
|
||
| b₂ = let (m₁, m₂, m₃, m₄, m₅, m₆) = (f, d, e, h, b, k) | ||
| α₁ = convert(BigInt, m₁ + m₅ - m₆) | ||
| α₂ = convert(BigInt, m₁ - m₅ + m₆) | ||
| α₃ = convert(BigInt, -m₁ + m₅ + m₆) | ||
| β₁ = convert(BigInt, m₁ + m₅ + m₆) | ||
| β₂ = convert(BigInt, m₂ + m₄ + m₆) | ||
| β₃ = convert(BigInt, m₃ + m₄ + m₅) | ||
| β₀ = convert(BigInt, m₁ + m₂ + m₃) | ||
| wei9jbracket(α₁, α₂, α₃, β₁, β₂, β₃, β₀) | ||
| end | ||
|
|
||
| b₃ = let (m₁, m₂, m₃, m₄, m₅, m₆) = (h, j, g, a, d, k) | ||
| α₁ = convert(BigInt, m₁ + m₅ - m₆) | ||
| α₂ = convert(BigInt, m₁ - m₅ + m₆) | ||
| α₃ = convert(BigInt, -m₁ + m₅ + m₆) | ||
| β₁ = convert(BigInt, m₁ + m₅ + m₆) | ||
| β₂ = convert(BigInt, m₂ + m₄ + m₆) | ||
| β₃ = convert(BigInt, m₃ + m₄ + m₅) | ||
| β₀ = convert(BigInt, m₁ + m₂ + m₃) | ||
|
There was a problem hiding this comment. Was there a specific reason for calling the last one There was a problem hiding this comment. Also, the construction of the b₁ = wei9jbracket(a, b, c, f, j, k)
b₂ = wei9jbracket(f, d, e, h, b, k)
b₃ = wei9jbracket(h, j, g, a, d, k)There was a problem hiding this comment.
No real reason - loosely the last one stood out to me as a break in a pattern since it doesn't have a matching αᵢ (wrt. the There was a problem hiding this comment.
wrt. the |
||
| wei9jbracket(α₁, α₂, α₃, β₁, β₂, β₃, β₀) | ||
| end | ||
|
|
||
| p * b₁ * b₂ * b₃ | ||
| end | ||
|
|
||
| return convert(BigInt, s) | ||
| end | ||
|
|
||
| # Wei square bracket terms appearing the 9j series | ||
| function wei9jbracket(α₁::BigInt, α₂::BigInt, α₃::BigInt, β₁::BigInt, β₂::BigInt, β₃::BigInt, β₀::BigInt) | ||
| p = max(β₁, β₂, β₃, β₀) | ||
| q = min(α₁ + β₂, α₂ + β₃, α₃ + β₀) | ||
| lrange = p:q | ||
| T = PrimeFactorization{eltype(eltype(factorialtable))} | ||
|
|
||
| terms = Vector{T}(undef, length(lrange)) | ||
| for (j, l) in enumerate(lrange) | ||
| b₁ = iseven(l) ? copy(primebinomial(big(l + 1), l - β₁)) : neg!(copy(primebinomial(big(l + 1), l - β₁))) | ||
| b₂ = primebinomial(α₁, l - β₂) | ||
| b₃ = primebinomial(α₂, l - β₃) | ||
| b₄ = primebinomial(α₃, l - β₀) | ||
|
|
||
| terms[j] = mul!(mul!(mul!(b₁, b₂), b₃), b₄) | ||
| end | ||
| total = sumlist!(terms) | ||
| return total | ||
| end | ||
|
|
||
| function _precompile_() | ||
| @assert precompile(wigner3j, (Type{Float64}, Int, Int, Int, Int, Int, Int)) | ||
| @assert precompile(wigner6j, (Type{Float64}, Int, Int, Int, Int, Int, Int)) | ||
| @assert precompile(wigner9j, (Type{Float64}, Int, Int, Int, Int, Int, Int, Int, Int, Int)) | ||
| @assert precompile(wigner3j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int})) | ||
| @assert precompile(wigner6j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int})) | ||
| @assert precompile(wigner9j, (Type{BigFloat}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int}, Rational{Int})) | ||
| @assert precompile(wigner3j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt)) | ||
| @assert precompile(wigner6j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt)) | ||
| @assert precompile(wigner9j, (Type{RRBig}, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt, HalfInt)) | ||
| end | ||
| _precompile_() | ||
|
|
||
|
|
||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have to admit that I am confused by how the permutations are encoded here and so how this code functions exactly. Would it be possible to follow the 3j and 6j philosophy and try to impose a canonical order. For example, given all the symmetries of the 9j symbol, it seems like one could always impose:
j₁ <= j₂ <= j₄ <= min(j₃, j₅, j₆, j₇, j₈, j₉)When the inequalities are strict, I believe this fixes the order completely, i.e. any permuation or reflection would destroy this order. However, when there are equalities, further conditions could be imposed.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was following the approach mentioned by the wigxjpf authors in this slide-deck (3rd slide from the end), where the symmetry checks for the 9j are done explicitly for the 72 equivalent permutations+transposition (up to sign) against the cache, rather than using a canonical ordering like the 3j or 6j.
The rationale there being (I believe) that unlike what happens in the 3j or 6j cases, the semimagic square symmetries give an under-constraining 6 conditions on 9 linearly independent entries rather than over-constraining the system to give a single canonical ordering.
My implementation could probably be cleaner or made more readable, but I'm using that Julia arrays can be indexed with vectors to get a permutation, e.g.
and then hard-coded all the column and row permutations accordingly (in lexicographic order) as well as the parity of each permutation. The transpositions get checked too but aren't encoded in the
_perms9jtable of indices, hence thekkandkkTcases.I think you're right though, and like your suggested condition as a starting point; it might be plausible to make a choice of ordering that amounts to a choice of the last few constraints, and handle a few cases to do better than checking all permutations. I'm out of time to spare today, so I'll sit down with this one later in the week to work out the cases (for when there are equalities) and how to cast an arbitrary symbol into the resulting form to see exactly how this can be improved!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Great. Looking at the slides, I indeed remember that for 3j and 6j you can go further and really compute a unique linear index for each of the symbols (such that those that are related up to a sign by symmetries have indeed the same linear index). That's probably too much to ask for the 9j symbols, but I don't use that anyway. I simply use the symmetries to reduce them to canonical order, but then still rely on a hash table.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've been busier than planned, so I've only been able to take a brief of a stab at a canonical form. The condition you suggest as a starting point unfortunately doesn't work: if a symbol can be canonicalized to meet it, then it's a good condition in that any row or col swap or transposition would break it - however, not all symbols can be canonicalized to satisfy that condition in the first place.
The issue is that anything that sets a particular entry over-constrains the permutations that can follow. e.g. in the suggested canonicalization
j₁must be made the smallest number in the symbol but doing so disallows the 1st row and 1st column from being involved in future row or column swaps to set the other entries, but in order to move the 2nd smallest number intoj₂it must have landed in the 1st row or the 1st column during the permutation that setj₁- which there's no guarantee it will have. Because of patterns like this I'm finding most canonicalization schemes I can come up with in the same vein require specifying many cases/conditions and get a bit ugly (and aren't looking much better than just checking the permutations).I'll see if I can find some better way to approach canonicalization over the week here - apologies for the delay.