J-scheme 配对近似代码 PandaWarrior 对壳模型有效相互作用的处理
1. isospin 转 pn 格式
参考上一篇博客，里面有 isospin 格式与 pn 格式的定义。因为历史原因，PandaWarrior代码使用的 pn 格式不同于 xpn 或者 upn，需要自己从 isospin 格式转换过来。上一篇博客也记录了 pn/xpn/upn 的不同，所以也容易写一个代码，让代码兼容其中任何一种格式。以后有时间了可以做这件事，这件事现在优先级不是很高。
做这个转换的代码在 PandaWarrior/src/interaction-handler/interactions/ 中，运行为
sh turn-int-isospin-pn.sh
脚本做两件事：
- 读取 input/GME.int 中的相互作用，每个两体格式为：
a b c d J T V
代码会做一点 sorting，比如 \(V(abcd;JT) = (-1)^{j_a + j_b + J + T}V(bacd;JT)\)，这种对称性在上一篇博客中有。
sort 之后的结果存在 input/GME_sort.int 中。 - isospin -> pn 格式转换。
这是靠 PandaWarrior/src/interaction-handler/interactions/code/turn-int-isospin-pn.cpp 中的代码实现的，这个代码把 T=1 的相互作用分别 copy 成 pp 和 nn 的相互作用，大部分代码是转 pn 相互作用，关键部分摘录如下。旧代码写的不好看，但是能用就行。
for(int i=0;i<num_GME;i++){
int a=GME_a[i];
int b=GME_b[i];
int c=GME_c[i];
int d=GME_d[i];
int I=2*GME_I[i];
int T=2*GME_T[i];
double V=GME_V[i];
for(int k=0;k<num_pn;k++){
if(GME_a_pn[k]==a //V(abcd;I)
&&GME_b_pn[k]==b+nj_p
&&GME_c_pn[k]==c
&&GME_d_pn[k]==d+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V;
}
if(a!=b){
if(GME_a_pn[k]==b //V(bacd;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==a+nj_p
&&GME_c_pn[k]==c
&&GME_d_pn[k]==d+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[a-1]+sh[b-1]+I+T)/2);
}
if(c!=d){
if(GME_a_pn[k]==b //V(badc;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==a+nj_p
&&GME_c_pn[k]==d
&&GME_d_pn[k]==c+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[a-1]+sh[b-1]+sh[c-1]+sh[d-1])/2);
}
if(GME_a_pn[k]==a //V(abdc;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==b+nj_p
&&GME_c_pn[k]==d
&&GME_d_pn[k]==c+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[c-1]+sh[d-1]+I+T)/2);
}
}
}
else if(c!=d){
if(GME_a_pn[k]==a //V(aadc;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==b+nj_p
&&GME_c_pn[k]==d
&&GME_d_pn[k]==c+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[c-1]+sh[d-1]+I+T)/2);
}
}
if( min(a,b) != min(c,d) || max(a,b) != max(c,d) ){
if(GME_a_pn[k]==c //V(cdab;I)
&&GME_b_pn[k]==d+nj_p
&&GME_c_pn[k]==a
&&GME_d_pn[k]==b+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V;
}
if(a!=b){
if(GME_a_pn[k]==c //V(cdba;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==d+nj_p
&&GME_c_pn[k]==b
&&GME_d_pn[k]==a+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[a-1]+sh[b-1]+I+T)/2);
}
if(c!=d){
if(GME_a_pn[k]==d //V(dcba;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==c+nj_p
&&GME_c_pn[k]==b
&&GME_d_pn[k]==a+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[a-1]+sh[b-1]+sh[c-1]+sh[d-1])/2);
}
if(GME_a_pn[k]==d //V(dcab;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==c+nj_p
&&GME_c_pn[k]==a
&&GME_d_pn[k]==b+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[c-1]+sh[d-1]+I+T)/2);
}
}
}
else if(c!=d){
if(GME_a_pn[k]==d //V(dcaa;I) omitted in usdb.int is collected in GME_V_pn
&&GME_b_pn[k]==c+nj_p
&&GME_c_pn[k]==a
&&GME_d_pn[k]==a+nj_p
&&GME_I_pn[k]==I/2){
GME_V_pn[k]+=V*theta((sh[c-1]+sh[d-1]+I+T)/2);
}
}
}
}
}
在这一块代码之前，已经定好了有哪些 \(V_{pn}(abcd;J V=?)\)，这段代码从isospin相互作用中解析出来\(V_{pn}\)中的那些\(V\)值。
过程有些trivial，就是挨个找 \(V(abcd;JT), V(bacd;JT), V(abdc;JT), ..., V(dcba;JT)\)，因为isospin相互作用中有些省略，所以需要挨个确认。
2. pn -> P+Q形式
\(V_{pp}, V_{nn}\)天然就是 Pairing 形式，主要是 \(V_{pn}\)需要转换，这个转换公式在 PandaWarrior的手册中有，这里只记录以下代码的思路。
- spo -> onebody
spo 是 single particle orbit 类
onebody 是 单体算符类
spo的一个成员函数 gnr_1body 会生成一个 onebody 的 vector，存储可能存在的各种非集体单体算符 - xpn_int::gnr_Vcoef
上一篇博客说了，xpn_int 这个名字完全是误解，应该是 pn_int，以后会更正。gnr_Vcoef是一个成员函数，解析 pn 相互作用，生成 QQ 形式相互作用的强度。
3. 效果演示：
isospin 格式的 jun45.int 转成 pn 格式是
!SP 1 0 3 5/2 1
!SP 2 1 1 3/2 1
!SP 3 1 1 1/2 1
!SP 4 0 4 9/2 1
!SP 5 0 3 5/2 -1
!SP 6 1 1 3/2 -1
!SP 7 1 1 1/2 -1
!SP 8 0 4 9/2 -1
constant = 0.0
-8.708700, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, // monopole terms, including s.p.e.
0.0, -9.828000, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, -7.838800, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, -6.261700, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, -8.708700, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, -9.828000, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -7.838800, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.261700,
334 -8.708700 -9.828000 -7.838800 -6.261700 -8.708700 -9.828000 -7.838800 -6.261700
2 2 2 2 0 1 -0.649200
2 2 2 2 2 1 0.245900
2 1 2 2 2 1 -0.353600
2 3 2 2 2 1 -0.593200
1 1 2 2 0 1 -0.840400
1 1 2 2 2 1 -0.394900
1 3 2 2 2 1 -0.231200
3 3 2 2 0 1 -1.215300
...
后面还有好多行，一共334个两体相互作用。计算 cu58 (p1n1) 的能谱得到
States E Ex J pi ?th
1 -22.1772 0 1+ 1th
2 -21.7893 0.387923 3+ 1th
3 -21.4466 0.730632 0+ 1th
4 -20.7007 1.47658 1+ 2th
5 -20.4653 1.71193 2+ 1th
6 -20.1488 2.02846 2+ 2th
7 -20.1132 2.06404 4+ 1th
8 -19.7789 2.39839 3+ 2th
9 -19.2654 2.91183 2+ 3th
10 -19.1301 3.04713 5+ 1th
11 -18.9386 3.23868 4+ 2th
12 -18.7368 3.44042 1+ 3th
13 -18.4782 3.69906 2+ 4th
14 -18.4426 3.73467 0+ 2th
15 -18.3726 3.80464 2- 1th
16 -18.2126 3.96463 3+ 3th
17 -18.154 4.02323 1+ 4th
18 -18.1407 4.03651 3+ 4th
19 -18.1381 4.03917 1+ 5th
20 -17.9441 4.2331 6- 1th
...
后面还有些能态，一共62个态。这些结果都与 BigStick(用 upn 格式的 jun45pn.int) 对照了，完全是一致的。另外我也测试了p1n2, p2n3 两个 cases (with scaling)，结果都与 BigStick 完全对上。
这说明咱的代码弄对了（过程中更正了PandaWarrior中 src/class.cpp Ln31 关于求余的一个bug: (-1)%2 =-1 而不是1）。但是凭良心说，咱这代码可真难用。可以考虑以下调整：
isospin -> pn 的代码：
- 输入文件 input/GME.int, input/pn.sp 可以在 STDIN 中指定文件名，而不是默认的固定名字。
- 输出文件可以取默认名字，但是代码得在 STDOUT 中提示用户。
- 输入、输出格式可以弄得更 compatible，输入一般没有 monopole，那就根 upn 的格式弄得 compatible吧，好记。
PandaWarrior 计算代码：
- 输入、输出文件同上1、2。
- 可以由用户指定在屏幕上显示多少个态。
- 兼容 upn/xpn/isospin 格式，即格式由用户指定。
- 还需要check Pandya transformation，在我印象中 particle-particle, hole-hole转的是对的，但是 particle-hole, hole-particle 好像有点问题。
- 以 pn_int 类为核心，把相互作用有关的功能都写成它的成员函数，然后把这个类移植进 mNPA, PVPC。
- STDIN 每行允许注释
