joeydumont/wignersymbols

A C++ library to compute the Wigner 3j- and 6j- symbols (with FORTRAN backend).

Wigner Symbols

A C++ ensemble of functions to compute the Wigner 3j- and 6j- symbols. It reimplements the algorithm designed by Schulten and Gordon in C++, but also contains the original Fortran implementation. It can act either as a complete C++ replacement to the original Fortran implementation, or a C++ interface to them. See the API docs for details. It can either compute an array of Wigner 3j or 6j symbols, or a single coefficient. It also computes the Clebsch-Gordan coefficients.

Note that there is a third party Python wrapper available on GitHub.

Compilation Instructions

This library uses CMake to help the build process. First, download the source code. It is recommended to create a separate directory for building, i.e.

mkdir build/
cd build/

Then, run

cmake .. && make && sudo make install

By default, the library is installed to /usr/lib/ and the include files are in /usr/include/. To install to another directory, say /usr/local/, use the command-line argument

cmake -DCMAKE_INSTALL_PREFIX:PATH=/usr/local && make && sudo make install

API documentation

We list the user-facing functions that compute the Wigner symbols. The functions are behind the namespace WignerSymbols. Both the C++ implementation (my own) and the original Fortran implementation (by Schulten and Gordon) are provided by this package. In my own tests, I had found that the original Fortran implementation provided only single precision. This might be compiler-dependent, so you might have better luck. The Fortran version is also somewhat faster (10% faster, approximately). In any case, this program provides either a complete C++ replacement, or a C++ interface to the evaluation of Wigner symbols.

C++ implementation

• std::vector<double> wigner3j(double l2, double l3, double m1, double m2, double m3)
  Computes Wigner 3j symbols with all possible values of l1. Returns an std::vector<double> with the coefficients sorted by increasing values of l1.
• double wigner3j(double l1, double l2, double l3, double m1, double m2, double m3)
  Computes a specific Wigner 3j symbol.
• double clebschGordan(double l1, double l2, double l3, double m1, double m2, double m3)
  Computes a specific Clebch-Gordan coeffcient.
• std::vector<double> wigner6j(double l2, double l3, double l4, double l5, double l6)
  Computes Wigner 6j symbols with all possible values of l1. Returns an std::vector<double> with the coefficients sorted by increasing values of l1.
• double wigner6j(double l1, double l2, double l3, double l4, double l5, double l6)
  Computes a specific Wigner 6j symbol.

Fortran implementation

• std::vector<double> wigner3j_f(double l2, double l3, double m1, double m2, double m3)
  Computes Wigner 3j symbols with all possible values of l1.
• double wigner3j_f(double l1, double l2, double l3, double m1, double m2, double m3)
  Computes a specific Wigner 3j symbol.
• double clebschGordan_f(double l1, double l2, double l3, double m1, double m2, double m3)
  Computes a specific Clebch-Gordan coeffcient.
• std::vector<double> wigner6j_f(double l2, double l3, double l4, double l5, double l6)
  Computes Wigner 6j symbols with all possible values of l1.
• double wigner6j_f(double l1, double l2, double l3, double l4, double l5, double l6)
  Computes a specific Wigner 6j symbol.

Bibliography

• K. Schulten and R. G. Gordon, Recursive evaluation of 3j and 6j coefficients, Comput. Phys. Commun. 11, 269–278 (1976). DOI: 10.1016/0010-4655(76)90058-8
• K. Schulten, Exact recursive evaluation of 3j- and 6j-coefficients for quantum-mechanical coupling of angular momenta, J. Math. Phys. 16, 1961 (1975). DOI: 10.1063/1.522426.
• J. Luscombe and M. Luban, Simplified recursive algorithm for Wigner 3j and 6j symbols, Phys. Rev. E 57, 7274–7277 (1998). DOI: 10.1103/PhysRevE.57.7274.

wignersymbols's Issues

Algorithm underflows when computing Wigner 3j symbols with large quantum numbers

When the quantum numbers become large, the relative sizes of the Wigner symbols can span tens of orders of magnitude. This frequently leads to overflow in the recurrence relations, which we deal with by dividing by a large number. However, this leads to underflows in the algorithm, resulting in the algorithm returning 0s instead of valid Wigner symbols. See

Wigner 3j symbol                 Mathematica       WignerSymbols
[ 751  856 1200  464 -828  364] -9.41731061215e-58 -0.0
[ 841  379 1011 -631  313  318] -2.44096504011e-41 -0.0
[ 570 1007 1392  327 -933  606] -1.74376347733e-98 0.0
[ 970  727 1202  533 -663  130] -6.93009562166e-12 0.0
[ 905  919 1670  869 -594 -275] -3.48516309858e-195 -0.0
[ 895  574 1392  793 -365 -428] -1.41868655509e-146 -0.0

for examples of this.

Wigner 3j returns 0s instead of actual input, but yields good inputs on permutation of the l's

Thanks for Jeff Zheng for pointing this out to me via email (@jeffzhen). This comes mostly verbatim for his e-mail to me.

I found some quite confusing bug where it returns 0 (where Mathematica returns non-zero). In all cases, the whole list of values for all l1 are 0 or -0. However, if I permute the 1,2,3 indices to say 2,3,1, I always get the correct answer. Here's one example:
[l1, l2, l3, m1, m2, m3] = [ 325 999 1221 280 899 -1179] Mathematica answer: 0.000163542255419
call (l1, l2, l3, m1, m2, m3): 0
call (l2, l3, m1, m2, m3): all 0's returned
call (l2, l3, l1, m2, m3, m1): 0.000163542255418

I appended more examples at the end. Would you mind take a quick look and see what is causing this? Since your code can clearly handle l numbers much larger than these, and it can return correct answer by just permuting the indices, my guess is that it is some benign bug?

More examples:
[ 693 896 1371 513 -838 325] 3.82641344761e-140 -0.0 0.0175918809377 3.82641344761e-140
[ 772 874 1231 442 756 -1198] 0.00171146251356 0.0 0.0 0.00171146251356
[ 842 996 1714 376 979 -1355] -9.12730017838e-14 -0.0 0.0 -9.12730017838e-14
[ 101 987 1084 -80 -935 1015] -0.00274486614091 -0.0 0.0 -0.00274486614091
[ 51 1003 978 -32 993 -961] -0.00561491486489 -0.0 0.0 -0.00561491486489
[ 217 1008 1107 -76 -987 1063] -0.00129163704248 -0.0 0.0 -0.00129163704248
[ 408 894 954 -114 -840 954] 0.000141186572123 0.0 0.0 0.000141186572123
[ 808 980 1734 341 954 -1295] -1.13581014918e-35 -0.0 0.0 -1.13581014918e-35
[ 827 1020 1570 590 902 -1492] -0.00091099434249 -0.0 0.0 -0.00091099434249

Add topic tags

I suggest adding the topics quantum-mechanics, wigner-symbols, physics in the About section, as explained at https://docs.github.com/en/repositories/managing-your-repositorys-settings-and-features/customizing-your-repository/classifying-your-repository-with-topics

Wigner 3j returns 0s instead of actual output

Thanks again to @jeffzhen for pointing this out.

In the computation of the Wigner symbol WignerSymbols::wigner3j(529, 992, 1243, 196, -901, 705), the function returns 0 (and the vector function for all allowed l1s returns a string of 0s) while the actual output is close to 1.97986e-18.

This is due to an overflow of the normalization condition. Will need to be fixed.

License version in copyright statements

The copyright statements in many of the source files simply say "This source code is subject to the terms of the GNU Lesser Public License" but do not say which version. Also, the official name of the license is "GNU Lesser General Public License".

The LICENSE file in the project root is for the LGPL version 3, so I presume that is what was intended, but you still don't say whether the license is "version 3 only" or "version 3 or any later version"?