写一个脚本，循环运行 Projected-Hartree-Fock程序，然后提取结果

知识点：

• 脚本可以传参数进去，就和 c/c++ 的代码可以传给 main 函数参数一样。

Calvin Johnson 和他以前的学生慢慢写出来的 Projected Hartree Fock 程序，Hartree Fock 的部分叫做 Sherpa，意思是夏尔巴人，是攀登珠穆朗玛峰的引路人（我就知道，Calvin Johnson起名字会费劲脑汁搞得有意思），投影的部分就做 LAMP，大概是灯光照耀道路。
下面的脚本叫做sd_nucleus_phf.sh，做 sd 壳中某个核在 usdb 下的 projected hartree fock：

#!/bin/bash

# These parameters should be passed to this script: sps, Np, Nn, int, Acore

sps=$1
Np=$2
Nn=$3
int=usdb  # isospin format, by default
Acore=16

nucleus="p"$Np"n"$Nn
A=$(( $Acore + 2 ))
B=$(( $Acore + $Np + $Nn ))
C=0.3
Jmax=20

echo "Start SHERPA"

rm $nucleus.sd   # 如果这些文件存在，下面的代码运行时会多问一个问题，导致准备的 input 对应不上。
rm $nucleus.res

./sherpa19.x << input
$sps    # name .sps file
$Np     # valence Z
$Nn     # valence N
n     # no external field
i     # isospin-format interaction
1     # 1 interaction file
$int     # name of interaction file
1 $A $B $C   # scaling of interaction matrix elements
n
test
19911
1000     # max number of iterations
15       # number of solutions, each starting from random
-1     # choose lowest energy
w   # write out to a file
$nucleus  # writes out to file xx.sd
$int
x   # leave HF menu
r   # go to RPA menu
e   # compute correlation energy
x
x
input

mv hf.basis $nucleus"hf.basis"
echo Hartree fock basis found in $nucleus"hf.basis"

echo "Start LAMP to project out states"

./lamp.x << INPUT
$sps 
$Np $Nn
1       # number of slater determinants
$nucleus
$Jmax   ! Jmax
n       ! write norm probabilities to a file or not
y       ! go on to energies or not
$int
1 $A $B $C      # scaling
end
0.0001  # tolerance for accept Js
20      # number of states to print onto screen
0.001   # tolerance for norm
y       # write out to file or not
$nucleus  # xx.res will be written into
n       # run with other parameters or not
INPUT

mv $nucleus.res sd/

echo "finished LAMP"

然后再写一个脚本sd_loop_phf.sh，循环遍历 (sd) 壳中所有 (N geq Z) 的偶偶核：

#!/bin/bash

# calculate spectrum of sd-shell even-even nuclei with PHF code

echo -n "Start running, "
echo $(date)
for(( pN=0; pN<=6; pN++ ))
do
        for(( nN=$pN; nN<=6; nN++ ))
        do
                Np=$(($pN * 2 ))
                Nn=$(($nN * 2 ))
                echo "Np = $Np, Nn = $Nn"
                bash sd_nucleus_phf.sh sd $Np $Nn usdb 16

        done
done

跑完以后会得到 sd/文件夹中的一系列pxnx.res文件。上一个帖子中有一个处理res文件的c++代码，稍作修改，即可对 PHF 的 res文件进行处理，可以得到 PHF 的 (sd) 壳所有偶偶核基态能量：

void GetPHFEgs( double ** PHFEgs ){
        double Egs;
        string reshead = "sd_phf/p", restail=".res", res;
        for(int Np = 0; Np <=12; Np +=2 ){
                for(int Nn = Np; Nn <= 12; Nn +=2 ){
                        res = reshead;
                        if( Np < 10 ) res += (char)(48+Np);
                        else res = res + (char)(48 + Np/10) + (char)(48 + Np%10);
                        res += "n";
                        if( Nn < 10 ) res += (char)(48+Nn);
                        else res = res + (char)(48 + Nn/10) + (char)(48 + Nn%10);
                        res += restail;
                        cout<<"Np = "<<Np<<", Nn = "<<Nn<<", res = "<<res<<endl;

                        double Egs;
                        ifstream fin(res);
                        if( !fin ) continue;
                        stringstream strm; string line, headword;
                        while( !fin.eof() && headword != "State" ){
                                getline(fin, line); strm.clear(); strm.str(line); strm>>headword;
                        }
                        getline(fin, line); // PHF的结果里有条分割线，读一下
                        getline(fin, line); strm.clear(); strm.str(line); strm>>headword >> Egs;
                        fin.close();

                        cout<<"Egs = "<<Egs<<endl;
                        PHFEgs[ Np/2 ][ Nn/2 ] = Egs;
                        PHFEgs[ Nn/2 ][ Np/2 ] = Egs;
                }
        }
}

有了上一个帖子的壳模型基态能量，和这个帖子的PHF基态能量，就可以两者相减，得知 PHF 基态比 SM 高多少。