Skip to content

Commit

Permalink
add CG3spin
Browse files Browse the repository at this point in the history
  • Loading branch information
@0382
0382 committed Apr 1, 2024
1 parent d167b35 commit 93d9f8a
Show file tree
Hide file tree
Showing 4 changed files with 83 additions and 13 deletions.
4 changes: 3 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,9 @@ double fast_binomial(int n, int k);
// CG coefficient
double CG(int dj1, int dj2, int dj3, int dm1, int dm2, int dm3);
// CG coefficient for two spin-1/2, equivalent to `CG(1, 1, 2*S, ds1, ds2, ds1+ds2)`, and faster
double CGspin(int ds1, int ds2, int S);
double CGspin(int dm1, int dm2, int S);
// <S12,M12|1/2,m1;1/2,m2><S,M|S12,M12;1/2,m3>
double CG3spin(int dm1, int dm2, int dm3, int S12, int dS)
// CG coefficient with m1 == m2 == m3 == 0
double CG0(int j1, int j2, int j3);
// Wigner 3j symbol
Expand Down
4 changes: 3 additions & 1 deletion README_zh.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,8 @@ double fast_binomial(int n, int k);
double CG(int dj1, int dj2, int dj3, int dm1, int dm2, int dm3);
// 两个 1/2 自旋的CG系数
double CGspin(int ds1, int ds2, int S);
// 三个 1/2 自旋两次耦合 <S12,M12|1/2,m1;1/2,m2><S,M|S12,M12;1/2,m3>
double CG3spin(int dm1, int dm2, int dm3, int S12, int dS)
// CG 系数特殊情况 m1 == m2 == m3 == 0
double CG0(int j1, int j2, int j3);
// 3j系数
Expand Down Expand Up @@ -98,7 +100,7 @@ wigner_init(21, "2bjmax", 6);
wigner_init(21, "Jmax", 9);
```
表示我们体系中最大可能的角动量为`Jmax = 21`，只计算到9j系数。
表示我们体系中最大可能的角动量为`Jmax = 21`，计算到9j系数。
`"Jmax"`模式没有任何假定，它总是安全的。即使考虑到三体耦合，只要使用三体的总`Jmax`来`wigner_init`，就能够保证不会溢出，尽管这可能会浪费一些空间。
Expand Down
52 changes: 48 additions & 4 deletions WignerSymbol.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -566,15 +566,59 @@ inline void wigner_init(int num, std::string type, int rank) { wigner.reserve(nu
inline double fast_binomial(int n, int k) { return wigner.binomial(n, k); }

// CG coefficient for two spin-1/2
inline double CGspin(int ds1, int ds2, int S)
inline double CGspin(int dm1, int dm2, int S)
{
static constexpr double inv_sqrt_2 = 0.70710678118654752;
static constexpr double inv_sqrt_2 = 0.7071067811865476;
static constexpr double values[2][2][2] = {0.0, 1.0, -inv_sqrt_2, inv_sqrt_2, inv_sqrt_2, inv_sqrt_2, 0.0, 1.0};
if (unsigned(S) > 1)
return 0;
if (std::abs(ds1) != 1 || std::abs(ds2) != 1)
if (std::abs(dm1) != 1 || std::abs(dm2) != 1)
return 0;
return values[ds1 > 0][ds2 > 0][S > 0];
return values[dm1 > 0][dm2 > 0][S];
}

// <S12,M12|1/2,m1;1/2,m2><S,M|S12,M12;1/2,m3>
inline double CG3spin(int dm1, int dm2, int dm3, int S12, int dS)
{
static constexpr double values[2][2][2][2][2] = {
0.0, // (-1/2, -1/2, -1/2, 0, 1/2) -> 0
0.0, // (-1/2, -1/2, -1/2, 0, 3/2) -> 0
0.0, // (-1/2, -1/2, -1/2, 1, 1/2) -> 0
1.0, // (-1/2, -1/2, -1/2, 1, 3/2) -> 1
0.0, // (-1/2, -1/2, 1/2, 0, 1/2) -> 0
0.0, // (-1/2, -1/2, 1/2, 0, 3/2) -> 0
-0.816496580927726, // (-1/2, -1/2, 1/2, 1, 1/2) -> -sqrt(2/3)
0.5773502691896257, // (-1/2, -1/2, 1/2, 1, 3/2) -> sqrt(1/3)
-0.7071067811865476, // (-1/2, 1/2, -1/2, 0, 1/2) -> -sqrt(1/2)
0.0, // (-1/2, 1/2, -1/2, 0, 3/2) -> 0
0.408248290463863, // (-1/2, 1/2, -1/2, 1, 1/2) -> sqrt(1/6)
0.5773502691896257, // (-1/2, 1/2, -1/2, 1, 3/2) -> sqrt(1/3)
-0.7071067811865476, // (-1/2, 1/2, 1/2, 0, 1/2) -> -sqrt(1/2)
0.0, // (-1/2, 1/2, 1/2, 0, 3/2) -> 0
-0.408248290463863, // (-1/2, 1/2, 1/2, 1, 1/2) -> -sqrt(1/6)
0.5773502691896257, // (-1/2, 1/2, 1/2, 1, 3/2) -> sqrt(1/3)
0.7071067811865476, // ( 1/2, -1/2, -1/2, 0, 1/2) -> sqrt(1/2)
0.0, // ( 1/2, -1/2, -1/2, 0, 3/2) -> 0
0.408248290463863, // ( 1/2, -1/2, -1/2, 1, 1/2) -> sqrt(1/6)
0.5773502691896257, // ( 1/2, -1/2, -1/2, 1, 3/2) -> sqrt(1/3)
0.7071067811865476, // ( 1/2, -1/2, 1/2, 0, 1/2) -> sqrt(1/2)
0.0, // ( 1/2, -1/2, 1/2, 0, 3/2) -> 0
-0.408248290463863, // ( 1/2, -1/2, 1/2, 1, 1/2) -> -sqrt(1/6)
0.5773502691896257, // ( 1/2, -1/2, 1/2, 1, 3/2) -> sqrt(1/3)
0.0, // ( 1/2, 1/2, -1/2, 0, 1/2) -> 0
0.0, // ( 1/2, 1/2, -1/2, 0, 3/2) -> 0
0.816496580927726, // ( 1/2, 1/2, -1/2, 1, 1/2) -> sqrt(2/3)
0.5773502691896257, // ( 1/2, 1/2, -1/2, 1, 3/2) -> sqrt(1/3)
0.0, // ( 1/2, 1/2, 1/2, 0, 1/2) -> 0
0.0, // ( 1/2, 1/2, 1/2, 0, 3/2) -> 0
0.0, // ( 1/2, 1/2, 1/2, 1, 1/2) -> 0
1.0 // ( 1/2, 1/2, 1/2, 1, 3/2) -> 1
};
if (unsigned(S12) > 1 || unsigned(dS) > 3 || (dS % 2 == 0))
return 0;
if (std::abs(dm1) != 1 || std::abs(dm2) != 1 || std::abs(dm3) != 1)
return 0;
return values[dm1 > 0][dm2 > 0][dm3 > 0][S12][dS / 2];
}

inline double CG(int dj1, int dj2, int dj3, int dm1, int dm2, int dm3)
Expand Down
36 changes: 29 additions & 7 deletions test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,19 +190,19 @@ void test_Moshinsky()
void test_CGspin()
{
std::mt19937 gen(0);
std::uniform_int_distribution<int> dist(-1, 1);
constexpr int N = 10'000;
std::uniform_int_distribution<int> dist(-1, 1); // has invalid input
constexpr int N = 100'000;
bool has_error = false;
for (int i = 0; i < N; ++i)
{
int ds1 = dist(gen);
int ds2 = dist(gen);
int dm1 = dist(gen);
int dm2 = dist(gen);
int S = dist(gen);
double x = CGspin(ds1, ds2, S);
double y = CG(1, 1, 2 * S, ds1, ds2, ds1 + ds2);
double x = CGspin(dm1, dm2, S);
double y = CG(1, 1, 2 * S, dm1, dm2, dm1 + dm2);
if (std::abs(x - y) > 1e-12)
{
std::cerr << "ds1 = " << ds1 << ", ds2 = " << ds2 << ", S = " << S << std::endl;
std::cerr << "dm1 = " << dm1 << ", dm2 = " << dm2 << ", S = " << S << std::endl;
std::cerr << "CGspin = " << x << ", CG = " << y << std::endl;
has_error = true;
std::exit(1);
Expand All @@ -212,6 +212,28 @@ void test_CGspin()
{
std::cout << "test CGspin passed" << std::endl;
}
for (int i = 0; i < N; ++i)
{
int dm1 = dist(gen);
int dm2 = dist(gen);
int dm3 = dist(gen);
int S12 = dist(gen);
int dS = dist(gen) + 2;
double x = CG3spin(dm1, dm2, dm3, S12, dS);
double y = CG(1, 1, 2 * S12, dm1, dm2, dm1 + dm2) * CG(2 * S12, 1, dS, dm1 + dm2, dm3, dm1 + dm2 + dm3);
if (std::abs(x - y) > 1e-12)
{
std::cerr << "dm1 = " << dm1 << ", dm2 = " << dm2 << ", dm3 = " << dm3 << ", S12 = " << S12
<< ", dS = " << dS << std::endl;
std::cerr << "CG3spin = " << x << ", CG = " << y << std::endl;
has_error = true;
std::exit(1);
}
}
if (!has_error)
{
std::cout << "test CG3spin passed" << std::endl;
}
}

void test_lsjj()
Expand Down

0 comments on commit 93d9f8a

Please sign in to comment.