Merge pull request #111 from SotaYoshida/dev

Dev

src/chiEFTint/dict_LECs.jl

@@ -191,7 +191,7 @@ function dict_emn500n4lo()
     dLECs["D_3P2"] =   5.342585336
     dLECs["hD_3S1"] = -2.949089421
     dLECs["D_3S1"] = -20.793199632
-    dLECs["hD_3SD1"] = 1.3545478412
+    dLECs["hD_3SD1"] = 1.345478412 
     dLECs["D_3SD1"] = 2.176852098
     dLECs["D_3D1"] = -6.01826561
     dLECs["D_1D2"] = -1.545851484

src/chiEFTint/main_chiEFTint.jl

@@ -25,17 +25,21 @@ function make_chiEFTint(;is_show=false,itnum=1,writesnt=true,nucs=[],optimizer="
     BE_d_bare = Calc_Deuteron(chiEFTobj,to;io=io)
     @timeit to "renorm." SRG(chiEFTobj,to)
     BE_d_srg = Calc_Deuteron(chiEFTobj,to;io=io)
-    println("E(2H): bare = ",@sprintf("%8.5f", BE_d_bare),
-            " srg = ", @sprintf("%8.5f", BE_d_srg), " Diff.", @sprintf("%8.3e", BE_d_bare - BE_d_srg))
+    if chiEFTobj.params.srg
+        println("E(2H): bare = ",@sprintf("%9.6f", BE_d_bare),
+                " srg = ", @sprintf("%9.6f", BE_d_srg), " Diff.", @sprintf("%8.3e", BE_d_bare - BE_d_srg))
+    else
+        println("E(2H): bare = ",@sprintf("%9.6f", BE_d_bare))        
+    end
     HFdata = prepHFdata(nucs,ref,["E"],corenuc)
 
     if do_svd; target_LSJ = [[0,0,0,0],[0,2,1,1],[1,1,1,0],[2,2,0,2]]; svd_vmom(chiEFTobj,target_LSJ); end
 
     if write_vmom
         target_LSJ = [[0,0,1,1],[1,1,1,0],[1,1,0,1],[1,1,1,1],[0,0,0,0],[0,2,1,1],[3,3,1,3]]
         write_onshell_vmom(chiEFTobj,2,target_LSJ;label="pn"); write_onshell_vmom(chiEFTobj,3,target_LSJ;label="nn")
-        momplot(chiEFTobj,2,target_LSJ; fnlabel=ifelse(chiEFTobj.params.srg,"srg","bare"))
-        momplot(chiEFTobj,3,target_LSJ; fnlabel=ifelse(chiEFTobj.params.srg,"srg","bare"))
+        #momplot(chiEFTobj,2,target_LSJ; fnlabel=ifelse(chiEFTobj.params.srg,"srg","bare"))
+        #momplot(chiEFTobj,3,target_LSJ; fnlabel=ifelse(chiEFTobj.params.srg,"srg","bare"))
     end
 
     if do2n3ncalib #calibrate 2n3n LECs by HFMBPT
@@ -83,6 +87,7 @@ function construct_chiEFTobj(do2n3ncalib,itnum,optimizer,MPIcomm,io,to;fn_params
     f_ss!(opfs[2]);f_c!(opfs[3])
     ## prep. Gauss point for integrals
     ts, ws = Gauss_Legendre(-1,1,40)
+    ts, ws = Gauss_Legendre(-1,1,96)
     ## prep. for TBMEs        
     infos,izs_ab,nTBME = make_sp_state(params;io=io)
     println(io,"# of channels 2bstate ",length(infos)," #TBME = $nTBME")

src/chiEFTint/pionexchange.jl

@@ -5,7 +5,7 @@ calc. One-pion exchange potential in the momentum-space
 
 Reference: R. Machleidt, Phys. Rev. C 63 024001 (2001).
 """
-function OPEP(chiEFTobj,to;pigamma=false)
+function OPEP(chiEFTobj,to;pigamma=true)
     ts = chiEFTobj.ts; ws = chiEFTobj.ws; xr = chiEFTobj.xr; V12mom = chiEFTobj.V12mom
     dict_pwch = chiEFTobj.dict_pwch; 
     lsjs = chiEFTobj.lsjs; tllsj = chiEFTobj.tllsj
@@ -34,13 +34,13 @@ function OPEP(chiEFTobj,to;pigamma=false)
                     nfac = 1.0 / (xdwn * ydwn)
                     ree =  1.0 /sqrt(ex*ey) *freg(x,y,4)
                     f_sq!(opfs,xdwn,ydwn)
-                    t_fc = fff * ree *hc3 * nfac
+                    t_fc = fff  *hc3 * nfac * ree
                     tVs .= 0.0
                     cib_lsj_opep(t_fc,opfs,x,y,mpi02,1,J,pnrank,ts,ws,tVs,QLdict)
-                    if pnrank == 2 ## pp/nn
+                    if pnrank == 2 # charged pion exchange
                         cib_lsj_opep(t_fc,opfs,x,y,mpipm2,2,J,pnrank,ts,ws,tVs,QLdict;additive=true)
                         if pigamma
-                            t_fc_pig = t_fc 
+                            t_fc_pig = t_fc
                             cib_lsj_opep(t_fc,opfs,x,y,mpipm2,2,J,pnrank,ts,ws,tVs,QLdict,pigamma;additive=true,factor_pig=t_fc_pig)
                         end
                     end
@@ -49,7 +49,7 @@ function OPEP(chiEFTobj,to;pigamma=false)
                         tl,tlp,S,tJ = lsj[idx] 
                         if pnrank%2 == 1 && (tl+S+1)%2 != 1;continue;end
                         V12idx = get(tdict,tllsj,-1)
-                        if V12idx == -1;continue;end                        
+                        if V12idx == -1;continue;end
                         V12mom[V12idx][i,j] += tVs[idx]
                     end
                 end
@@ -59,36 +59,78 @@ function OPEP(chiEFTobj,to;pigamma=false)
     return nothing
 end
 
+"""
+pi-gamma correction Eq.(4) of U. van Kolck et al., Phys. Rev. Lett. 80, 4386 (1998).
+
+Note that 1/(1+beta^2) is ommited, since it is also included in the OPEP contribution,
+i.e. one needs to consider only the difference other than this factor.
+"""
 function fac_pig(beta,c5=0.0)
     return - (1.0-beta)^2 / (2*beta^2) * log(1+beta) +(1.0+beta)/(2*beta) -2.0*c5
 end
 
 """
 Ref: R. Machleidt, Phys. Rev. C 63, 024001 (2001).
+
+For pi-gamma correction term, which is introduced in e.g. EM500 interaction,
+one needs to evaluate integrals in (B11)~ terms in an explicit manner.
 """
-function cib_lsj_opep(fac_in,opfs,x,y,mpi2,nterm,J,pnrank,ts,ws,tVs,QLdict,pigamma=false;additive=false,factor_pig=1.0)
+function cib_lsj_opep(fac_in,opfs,x,y,mpi2,nterm,J,pnrank,ts,ws,tVs,QLdict,pigamma=false;additive=false,factor_pig=1.0)    
     x2 = x^2; y2 = y^2; z = (mpi2+x2+y2) / (2*x*y)
     QJ = QJm1 = 0.0
-    nfac = fac_in
+    nfac = fac_in 
     q2s = Float64[ ]
+    qfac = 1.0
     if pigamma
         q2s = zeros(Float64,length(ts))        
-        nfac = ifelse(J==0,0.0,factor_pig)
         for (i,t) in enumerate(ts)
             q2 = x2 + y2 -2.0*x*y*t
             beta = q2 / mpi2
-            q2s[i] = (fsalpha/pi) * fac_pig(beta)
+            nfac = fac_in  * fsalpha/pi      
+            q2s[i] = fac_pig(beta) 
         end
         QJ = QL_numeric_fac(z,J,ts,ws,q2s)
         if J>0;QJm1=QL_numeric_fac(z,J-1,ts,ws,q2s);end    
     else
         QJ = QL(z,J,ts,ws,QLdict)
         if J>0;QJm1=QL(z,J-1,ts,ws,QLdict);end
     end
+
+    # Using Eq.(B17)-(B20), which is valid for OPEP (not pi-gamma correction term)
     IJ0 = nfac * QJ
-    IJ1 = nfac * (z * QJ -delta(J,0)) #Eq. (B18)
-    IJ2 = nfac * (J*z* QJ + QJm1) /(J+1) #Eq. (B19)
-    IJ3 = nfac * sqrt(J/(J+1)) * (z* QJ - QJm1) #Eq. (B20)
+    IJ1 = nfac * (z * QJ  - delta(J,0) )  
+    IJ2 = nfac * (J*z* QJ + QJm1) /(J+1) 
+    IJ3 = nfac * sqrt(J/(J+1)) * (z* QJ - QJm1) 
+
+    if pigamma
+        QJ = QJm1 = 0.0 
+        for i=1:length(ws)
+            t = ts[i]
+            q2 = x2 + y2 -2.0*x*y*t
+            denom = (mpi2 +q2)
+            pj = Legendre(J,t) 
+            QJ += pj * q2s[i] * ws[i] / denom 
+            if J > 0 
+                QJm1 += Legendre(J-1,t) * q2s[i] * ws[i] / denom
+            end
+        end
+        QJ *= x * y
+        QJm1 *= x * y
+        IJ0 = nfac * QJ
+        QJ_ = 0.0 
+        for i=1:length(ws)
+            t = ts[i]
+            q2 = x2 + y2 -2.0*x*y*t
+            denom = mpi2 +q2
+            pj = Legendre(J,t) 
+            QJ_ += ts[i]* pj * q2s[i] * ws[i] / denom            
+        end
+        QJ_ *= x * y
+        IJ1 = nfac * QJ_
+        IJ2 = nfac * ( J * QJ_ + QJm1)/ (J+1) 
+        IJ3 = nfac * sqrt(J/(J+1)) * (QJ_ - QJm1)
+    end
+
     #Eq. (B28)
     v1 = opfs[1] * IJ0 + opfs[2] *IJ1
     v2 = opfs[3] * IJ0 + opfs[4] *IJ2
@@ -103,6 +145,7 @@ function cib_lsj_opep(fac_in,opfs,x,y,mpi2,nterm,J,pnrank,ts,ws,tVs,QLdict,pigam
     v34 = -sqrt(J*(J+1)) * (v3-v4)
     v56 =  sqrt(J*(J+1)) * (v5+v6)
     d2j1 = 1.0/(2*J+1)
+
     if nterm == 1
         phase = ifelse(pnrank==2,-1.0,1.0)
         tVs[1] = additive_sum(additive,tVs[1],v1 *phase)
@@ -119,11 +162,11 @@ function cib_lsj_opep(fac_in,opfs,x,y,mpi2,nterm,J,pnrank,ts,ws,tVs,QLdict,pigam
         tVs[1] = additive_sum(additive,tVs[1],ttis * v1)
         tVs[2] = additive_sum(additive,tVs[2],ttit * v2)
         tVs[3] = additive_sum(additive,tVs[3],d2j1 * ((J+1)* (ttis*v3) + J*(ttis*v4)-(ttis*v56)))
-        tVs[4] = additive_sum(additive,tVs[4],d2j1 * ( J*(v3*ttis) + (J+1)*(ttis*v4) +(ttis*v56)))
+        tVs[4] = additive_sum(additive,tVs[4],d2j1 * ttis * ( J*v3 + (J+1)*v4 +v56) )
         tVs[5] = additive_sum(additive,tVs[5],-d2j1 * ((ttis*v34)-(J+1)*(ttis*v5)+J*(ttis*v6)))
         tVs[6] = additive_sum(additive,tVs[6],-d2j1 * ((ttis*v34)+J*(ttis*v5)-(J+1)*(ttis*v6)))
     end
-    return nothing 
+     return nothing 
 end
 
 function additive_sum(TF::Bool,retv,inv)
@@ -206,7 +249,7 @@ end
 """
     tpe(chiEFTobj,to::TimerOutput)
 
-calc. two-pion exchange terms up to N3LO(EM) or N4LO(EMN)
+Function to calculate two-pion exchange terms up to N3LO(EM) or N4LO(EMN)
 
 The power conting schemes for EM/EMN are different;
 The ``1/M_N`` correction terms appear at NNLO in EM and at N4LO in EMN.
@@ -485,7 +528,7 @@ function Vs_term(chi_order,LoopObjects,w,tw2,q2,k2,Lq,Aq,nd_mpi,r_d145,Fpi2;EMN=
         else
             f_N3LO_Vs = gA4 / (32.0 * pi^2 * Fpi4)
             f_N3LO_2l_Vs = -gA2 * r_d145 /(32.0*pi^2 *Fpi4)            
-            # There is type in EM review (D.11) 5/8 -> 3/8.
+            # There is typo in EM review (D.11) 5/8 -> 3/8.
             # See also N. Kaiser, Phys. Rev. C 65 (2002) 017001.
             tmp_s += Lq * (k2 + 3.0*q2/8.0 + nd_mpi4 /w2) *f_N3LO_Vs
             tmp_s += w2 * Lq * f_N3LO_2l_Vs
@@ -533,7 +576,7 @@ function Ws_term(chi_order,LoopObjects,w,q2,Lq,Aq,nd_mpi,c4,Fpi2;EMN=false,useMa
         f_N3LO_a_Ws = c4^2 / (96.0*pi^2 *Fpi4)
         tmp_s += Lq * w2 * f_N3LO_a_Ws
         f_N3LO_b_Ws = - c4 / (192.0*pi^2 *Fpi4)
-        tmp_s += Lq * ( gA2 *(16.0*nd_mpi2 + 7.0 * q2) -w2) * f_N3LO_b_Ws
+        tmp_s += Lq * ( gA2 *(16.0*nd_mpi2 + 7.0 * q2) -w2) * f_N3LO_b_Ws #ci/MN
         f_N3LO_c_Ws = 1.0/ (1536.0* pi^2 * Fpi4)
         f_N3LO_2l_Ws= gA4 / (2048.0 * pi^2 * Fpi6)
         if EMN
@@ -603,7 +646,7 @@ function Vc_term(chi_order,w,tw2,q2,Lq,Aq,nd_mpi,c1,c2,c3,Fpi2,LoopObjects;EMN=f
         tmp_s += Lq *( brak6 * f_N3LO_f6_Vc)
         ### 2-loop corrections
         if EMN 
-           obj = LoopObjects.n3lo; ImV = obj.ImVc; tmp_s += Integral_ImV(obj.mudomain,q2,obj,ImV)
+            obj = LoopObjects.n3lo; ImV = obj.ImVc; tmp_s += Integral_ImV(obj.mudomain,q2,obj,ImV)
         else
             ## EM eq.(D.9)
             Mm2cor = Lq *(2.0*nd_mpi2^4 / w^4 + 8.0*nd_mpi2^3 / w2 -q2^2 -2.0*nd_mpi2^2 ) + nd_mpi2^3 /(2.0 *w2)
@@ -712,7 +755,7 @@ function Vls_term(chi_order,w,tw2,q2,Lq,Aq,nd_mpi,c2,Fpi2;EMN=false)
     end
     if  chi_order >= 3
         f_N3LO_a_Vls = c2 *gA2 /(8.0* pi^2 * Fpi4)
-        tmp_s += f_N3LO_a_Vls * w2 * Lq
+        tmp_s += f_N3LO_a_Vls * w2 * Lq 
         f_N3LO_b_Vls = gA4 /(4.0* pi^2 * Fpi4)
         if !EMN
             tmp_s += f_N3LO_b_Vls * Lq * (11.0/32.0 * q2 + nd_mpi4 / w2)
@@ -734,7 +777,7 @@ function Wls_term(chi_order,w,q2,Lq,Aq,nd_mpi,c4,Fpi2;EMN=false)
     if  chi_order >= 3
         f_N3LO_a_Wls = -c4 /(48.0*pi^2 * Fpi4)
         f_N3LO_b_Wls = 1.0 /(256.0*pi^2 * Fpi4) 
-        tmp_s += f_N3LO_a_Wls*Lq * (gA2 *(8.0*nd_mpi2 + 5.0*q2) + w2)
+        tmp_s += f_N3LO_a_Wls*Lq * (gA2 *(8.0*nd_mpi2 + 5.0*q2) + w2) #ci/MN 
         if !EMN
             tmp_s += f_N3LO_b_Wls*Lq * (16.0*gA2 *(nd_mpi2 + 3.0/8.0 *q2)
                      + 4.0/3.0 * gA4 *(4.0*nd_mpi4 /w2 -11.0/4.0 * q2 -9.0*nd_mpi2) -w2)

src/chiEFTint/struct_const_io.jl

@@ -713,6 +713,14 @@ end
 plot nn potential in partial wave over relative momentum space
 """
 function momplot(chiEFTobj,pnrank,llpSJ_s;fnlabel="",ctext="",fpath="",write_hdf5=true)
+    cpnrank = "pp"
+    if pnrank == 1
+        cpnrank = "pp"
+    elseif pnrank == 2
+        cpnrank = "pn"
+    elseif pnrank == 3
+        cpnrank = "nn"
+    end
     xr = chiEFTobj.xr_fm
     V12mom = chiEFTobj.V12mom
     tdict = chiEFTobj.dict_pwch[pnrank]
@@ -735,7 +743,6 @@ function momplot(chiEFTobj,pnrank,llpSJ_s;fnlabel="",ctext="",fpath="",write_hdf
         end
         for i =1:n_mesh;  tv[i] = V[i,i]; end
         if fpath != ""; tfdat = [xf,yfs[vidx]];end
-        #pw_plt(tx,xr,V,tv,pnrank,fnlabel;fdat=tfdat)
         if write_hdf5
             fn = "vmom_"*fnlabel*"_chiEFT_"*tx*"_"*string(pnrank)*".h5"
             h5open(fn,"w") do file
@@ -753,8 +760,8 @@ end
 function momplot_from_file(pnrank,tx)
     fn_bare = "vmom_bare_chiEFT_"*tx*"_"*string(pnrank)*".h5"
     fn_srg = "vmom_srg_chiEFT_"*tx*"_"*string(pnrank)*".h5"
-    if !isfile(fn_bare);println("file $fn_bare is not found");return nothing;end
-    if !isfile(fn_srg);println("file $fn_srg is not found");return nothing;end
+    if !isfile(fn_bare);println("file $fn_bare is not found, you need to run srg=false calculations to create 'vmom' files");return nothing;end
+    if !isfile(fn_srg);println("file $fn_srg is not found , you need to run srg=true calculations to create 'vmom' files");return nothing;end
     file = h5open(fn_bare,"r"); x_b = read(file,"x"); V_b = read(file,"V"); V_b_diag = diag(V_b);  close(file)    
     file = h5open(fn_srg,"r");  x_s = read(file,"x"); V_s = read(file,"V"); V_s_diag = diag(V_s); close(file)
     tls = ["pp", "pn", "nn"]

test/runtests.jl

@@ -1,7 +1,6 @@
 using NuclearToolkit
 using Test
 
-
 @testset  "NuclearToolkit.jl" begin
     @testset "generate NN potential" begin
         @test make_chiEFTint(;fn_params="optional_parameters.jl")
@@ -14,7 +13,7 @@ using Test
         sntf = "tbme_em500n3lo_barehw20emax2.snt.bin"
         ## HF-MBPT from snt/snt.bin
         @testset "HFMBPT results under bare EM500,hw20,e2,nmesh50" begin
-            Eref = [1.47561,-5.80114,0.39324] 
+            Eref = [1.493, -5.805, 0.395]
             HFobj1 = hf_main(nucs,sntf,hw,emax;return_obj=true,verbose=true)
             Es1 = [HFobj1.E0, HFobj1.EMP2, HFobj1.EMP3]
             println("Eref $Eref")
@@ -35,7 +34,7 @@ using Test
                 @test ((HFobj1.E0-HFobj2.E0)^2 + (HFobj1.EMP2-HFobj2.EMP2)^2 + (HFobj1.EMP3-HFobj2.EMP3)^2) < 1.e-6
             end
         end
-        Eref = -4.06623902 
+        Eref = -4.05224909
         @testset "IMSRG results under bare EM500,hw20,e2,nmesh50" begin
             IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,return_obj=true)
             Es = IMSRGobj.H.zerobody[1]
@@ -47,7 +46,7 @@ using Test
             @test abs(Eref - Es[1]) < 1.e-6
         end
         @testset "shell model calculation" begin
-            Eref = [ -10.737, -8.426]
+            Eref = [ -10.720, -8.410]
             vs_sntf = "vsimsrg_p-shell_coreHe4refHe4_He4_hw20e2_Delta0.0.snt";  n_eigen=2;targetJ=[]
             Ens = main_sm(vs_sntf,"Be8",n_eigen,targetJ)
             @test ((Eref[1]-Ens[1])^2 + (Eref[2] - Ens[2])^2) < 1.e-6