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

脚本做两件事:

  1. 读取 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 中。
  2. 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的手册中有,这里只记录以下代码的思路。

  1. spo -> onebody
    spo 是 single particle orbit 类
    onebody 是 单体算符类
    spo的一个成员函数 gnr_1body 会生成一个 onebody 的 vector,存储可能存在的各种非集体单体算符
  2. 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 的代码:

  1. 输入文件 input/GME.int, input/pn.sp 可以在 STDIN 中指定文件名,而不是默认的固定名字。
  2. 输出文件可以取默认名字,但是代码得在 STDOUT 中提示用户。
  3. 输入、输出格式可以弄得更 compatible,输入一般没有 monopole,那就根 upn 的格式弄得 compatible吧,好记。

PandaWarrior 计算代码:

  1. 输入、输出文件同上1、2。
  2. 可以由用户指定在屏幕上显示多少个态。
  3. 兼容 upn/xpn/isospin 格式,即格式由用户指定。
  4. 还需要check Pandya transformation,在我印象中 particle-particle, hole-hole转的是对的,但是 particle-hole, hole-particle 好像有点问题。
  5. 以 pn_int 类为核心,把相互作用有关的功能都写成它的成员函数,然后把这个类移植进 mNPA, PVPC。
  6. STDIN 每行允许注释
上一篇:OS L5-3: Contiguous Scheme and Fragmentation


下一篇:CFD杂谈-效率