ChiEFTint

Calc_Deuteron(chiEFTobj::ChiralEFTobject,to;io=stdout)

Function to calculate deuteron binding energy.

Gauss_Legendre(xmin,xmax,n;eps=3.e-16)

Calculating mesh points x and weights w for Gauss-Legendre quadrature. This returns x, w.

Legendre(n,x)

function to calculate Legendre polynomial $P_n(x)$

QL(z,J::Int64,ts,ws)

To calculate Legendre functions of second kind, which are needed for pion-exchange contributions, by Gauss-Legendre quadrature.

Rnl_all_ab(lmax,br,n_mesh,xr_fm)

Returns array for radial functions (prop to generalized Laguerre polynomials) HO w.f. in momentum space. Rnlk(l,n,k)=sqrt(br) * R(n,L,Z) *Z with Z=br*k (k is momentum in fm${}^{-1}$)

Vrel(chiEFTobj::ChiralEFTobject,to)

To define V in two-particle HO basis.

calc_Vmom!(pnrank,V12mom,tdict,xr,LEC,LEC2,l,lp,S,J,pfunc,n_reg,to;is_3nf=false)

calc. NN-potential for momentum mesh points

calculating coulomb contribution in HO base. Results are overwritten to Vcoulomb

construct_chiEFTobj(do2n3ncalib,itnum,optimizer,MPIcomm,io,to;fn_params="optional_parameters.jl")

It returns

• chiEFTobj::ChiralEFTobject parameters and arrays to generate NN (+2n3n) potentials. See also struct ChiralEFoObject.
• OPTobj (mutable) struct for LECs calibrations. It can be LHSobject/BOobject/MCMCobject/MPIMCMCobject struct.

def_sps_snt(emax,target_nlj)

Defining dicts for single particle states. One can truncate sps by specifying target_nlj, but it is usually empty so this function is usually not used.

freg(p,pp,n;Lambchi=500.0)

the regulator function used for NN or 2n3n contribution, Eq.(4.63) in EM's review: $f(p,p') = \exp [ -(p/\Lambda)^{2n} -(p'/\Lambda)^{2n} ]$

genLaguerre(n::Int,alpha,x)

returns generalized Laguaerre polynomials, $L^{\alpha}_n(x)$

get_twq_2b(emax,kh,kn,kl,kj,maxsps,jab_max)

returns infos, izs_ab, nTBME

• infos::Vector{Vector{Int64}} information of two-body channeles: { [$T_z$,parity,J,dim] }
• izs_ab::Vector{Vecrtor{Int64}} two-body kets: { [iza,ia,izb,ib] } where iz* stuff is isospin (-1 or 1) and i* stuff is label of sps.

For example, [-1,1,1,1] corresponds to |p0s1/2 n0s1/2>.

• nTBME::Int # of TBMEs

hw_formula(A,fnum)

empirical formula for harmonis oscillator parameter hw by mass number A fnum=2: for sd-shell, Ref. J. Blomqvist and A. Molinari, Nucl. Phys. A106, 545 (1968).

make_chiEFTint(;is_show=false,itnum=1,writesnt=true,nucs=[],optimizer="",MPIcomm=false,corenuc="",ref="nucl",Operators=[],fn_params="optional_parameters.jl",write_vmom=false,do_svd=false)

The interface function in chiEFTint. This generates NN-potential in momentum space and then transforms it in HO basis to give inputs for many-body calculations. The function is exported and can be simply called make_chiEFTint() in your script.

Optional arguments: Note that these are mainly for too specific purposes, so you do not specify these.

• is_show::Bool to show TimerOutputs
• itnum::Int number of iteration for LECs calibration with HFMBPT
• writesnt::Bool, to write out interaction file in snt (KSHELL) format. julia writesnt = false case can be usefull when you iteratively calculate with different LECs.
• nucs target nuclei used for LECs calibration with HFMBPT
• optimizer::String method for LECs calibration. "MCMC","LHS","BayesOpt" are available
• MPIcomm::Bool, to carry out LECs sampling with HF-MBPT and affine inveriant MCMC
• Operators::Vector{String} specifies operators you need to use in LECs calibrations
• fn_params::String path to file specifying the optional parameters
• write_vmom::Bool to write out in vmom partial wave channels

make_sp_state(chiEFTparams;io=stdout)

Defining the two-body channels or kets.

prepare_2b_pw_states(;io=stdout)

preparing two-body partical-wave channels, <Lp,S,J| |L,S,J>. For example, [pnrank,L,Lp,S,J] = [0, 0, 2, 1, 1] corresponds to proton-neutron 3S1-3D1 channel.

read_LECs(pottype)

read LECs for a specified potential type, em500n3lo,em500n3lo,emn500n4lo.

simply calculate $R_{nl}(p,b) = \sqrt{ \frac{2 n! b^3}{\Gamma(n+l+3/2)} (pb)^l e^{-p^2b^2/2} L^{l+1/2}_{n}(p^2b^2) }$

LO(chiEFTobj,to)

Calculate Leading Order (LO), Q^0` contact term.

N3LO(chiEFTobj,to)

Calculate Next-to-next-to-next-to Leading Order (N3LO), $Q^4$ contact term.

N4LO(chiEFTobj,to;n_reg=2)

Calculate Next-to-next-to-next-to-next-to Leading Order (N4LO) contact term.

NLO(chiEFTobj,to;n_regulator=2)

Calculate Next to Leading Order (NLO), $Q^2$ contact term. LECs=>C23S1,C23P0,C21P1,C23P1,C21S0,C23SD1,C23DS1,C23P2

The power for exponential regulator n_reg is 2 (default), For 3P1 channel, n_reg=3 is used as in the R.Machleidt's fortran code. I don't know the reason, but the LECs in EMN potential) may have been determined with n_reg=3 potential.

9xnthreads matrix to store sum of each tpe channal, C/T/S/LS/SigmaL

OPEP(chiEFTobj,to;pigamma=false)

calc. One-pion exchange potential in the momentum-space

Reference: R. Machleidt, Phys. Rev. C 63 024001 (2001).

Vc_term(chi_order,w,tw2,q2,Lq,Aq,nd_mpi,c1,c2,c3,Fpi2,LoopObjects;EMN=false,usingMachleidt=true,ImVerbose=false,calc_TPE_separately=false)

• NNLO: EM NNLO eq.(4.13), EKMN Eq.(C1) + relativistic correction Eq.(D7)
• N3LO:
  • f1term $c^2_i$: EM Eq.(D.1), EKMN Eq.(D1)
  • f6term $c_i/M_/N$: EM Eq.(D.4)

Vls_term(chi_order,w,tw2,q2,Lq,Aq,nd_mpi,c2,Fpi2;EMN=false)

Vt_term(chi_order,w,q2,k2,Lq,Aq,nd_mpi,r_d145,Fpi4)

• NLO: EM Eq.(4.10)
• NNLO: EM Eq.(4.15) & Eq.(4.22) => it_pi term
• N3LO: EM Eq.(D.11) (1/M^2_N` term), Eq.(D.23) (2-loop term)

Wc_term(chi_order,LoopObjects,w,tw2,q2,k2,Lq,Aq,nd_mpi,c4,r_d12,r_d3,r_d5,Fpi2;EMN=false,useMachleidt=true,calc_TPE_separately=false)

Ws_term(chi_order,LoopObjects,w,q2,Lq,Aq,nd_mpi,c4,Fpi2;EMN=false,useMachleidt=true,calc_TPE_separately=false)

• NNLO: EM Eq.(4.16) & Eq.(4.24) => it_pi term
• N3LO: EM Eq.(D.2) ($c^2_i$ term), Eq.(D.6) (c_i/M_Nterm), Eq.(D.12) (1/M^2_N`` term), Eq.(D.27) (2-loop term)

calc_LqAq(w,q,nd_mpi,usingSFR,LamSFR)

To calculate L(q) and A(q) appear in TPE contributions; L(q)=EM Eq.(4.11), A(q)=EM Eq.(D.28). If usingSFR, Lq&Aq are replaced by Eqs.(2.24)&(C3) in EKMN paper.

Ref: R. Machleidt, Phys. Rev. C 63, 024001 (2001).

n3lo_ImVc!(muvec,V,mpi,Fpi)

The expression can be found in eq.(D3) of EKMN paper.

n3lo_ImVsVt!(Vt,mpi,Fpi,r_d145,mudomain,ts,ws)

The expression can be found in eq.(D5) of EKMN paper and the second term in eq.(D5) is ommited. Note that LECs $\bar{d}_{14}-\bar{d}_{15}$ is given as d145 in NuclearToolkit.jl. Vt can be used for Vs=-q2*Vt

n3lo_ImWc!(V,mpi,Fpi,r_d12,r_d3,r_d5,r_d145,mudomain,ts,ws)

The expression can be found in eq.(D4) of EKMN paper.

n3lo_ImWsWt!(Wt,mpi,Fpi,mudomain,ts,ws)

The expression can be found in eq.(D6) of EKMN paper. Wt can be used for Ws=-q2*Wt

\[V_{C,S}(q) = -\frac{2q^6}{\pi} \int^{\tilde{\Lambda}}_{nm_\pi} d\mu \frac{\mathrm{Im} V_{C,S}(i\mu)}{\mu^5 (\mu^2+q^2)} \\ V_T(q) = \frac{2q^4}{\pi} \int^{\tilde{\Lambda}}_{nm_\pi} d\mu \]

tpe(chiEFTobj,to::TimerOutput)

calc. 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.

References

• EM: R. Machleidt and D.R. Entem Physics Reports 503 (2011) 1–7
• EMKN: D. R. Entem, N. Kaiser, R. Machleidt, and Y. Nosyk, Phys. Rev. C 91, 014002 (2015).

tpe_for_givenJT(chiEFTobj,LoopObjects,Fpi2,tmpLECs,J,pnrank,ts,ws,fff,dwn,nd_mpi,xr,pjs_para,gis_para,opfs_para,f_idx,tVs_para,lsj,tllsj_para,tdict,V12mom,tmpsum,to;calc_TPE_sep=true)

Calculating TPE contribution in a given momentum mesh point.

HObracket(nl, ll, nr, lr, n1, l1, n2, l2, Lam, d::Float64,dWs,tkey9j,dict9j_HOB,to)

To calc. generalized harmonic oscillator brackets (HOBs), $<<n_l\ell_l,n_r\ell_r:\Lambda|n_1\ell_1,n_2\ell_2:\Lambda>>_d$ from the preallocated 9j dictionary.

Reference:

• [1] B.Buck& A.C.merchant, Nucl.Phys. A600 (1996) 387-402
• [2] G.P.Kamuntavicius et al., Nucl.Phys. A695 (2001) 191-201

TMtrans(chiEFTobj::ChiralEFTobject,dWS,to;writesnt=true)

Function to carry out Talmi-Mochinsky transformation for NN interaction in HO space and to write out an sinput file.

get_nkey_from_abcdarr(tkey;ofst=1000)

To get integer key from an Int array (with length greater than equal 4)

kinetic_ob(nlj1, nlj2)

calc. kinetic one-body contribution \langle j_1 |T/\habr\omega| j_2 \rangle

kinetic_tb(nljtz1, nljtz2, nljtz3, nljtz4, J, dWS)

calc. kinetic two-body contribution $\langle j_1j_2|| -p_1 p_2/\hbar\omega ||j_3j_4\rangle_J$ using preallocated 6j-symbols.

dWS2n struct

it doesn't work for now

prep_wsyms()

preparing Clebsch-Gordan coefficients for some special cases: cg1s = (1,0,l,0|l',0), cg2s = (2,0,l,0|l',0)

red_nabla_j(nlj1,nlj2,d6j,key6j)

returns $b \langle j || \nabla || j_2\rangle$ Note that $l_1,l_2$ in nlj1&nlj2 are not doubled.

red_nabla_l(n1,l1,n2,l2)

returns $b \langle n_1,l_1|| \nabla ||n_2,l_2 \rangle$ in $l$-reduced matrix element.

associated Legendre Polynomials are calculated with AssociatedLegendrePolynomials.jl

vtrans(chiEFTobj,pnrank,izz,ip,Jtot,iza,ia,izb,ib,izc,ic,izd,id,nljsnt,V12ab,t5v,to)

Function to calculate V in pn-formalism:

\[\langle ab;JTz|V| cd;JTz \rangle = N_{ab} N_{cd} \sum_{\Lambda S \Lambda' S'} \sum_{n\ell N L}\sum_{n'\ell' N' L'} \sum_{J_\mathrm{rel}J'_\mathrm{rel}} [\Lambda][\Lambda'] \hat{S}\hat{S'} \hat{J}_\mathrm{rel}\hat{J}'_\mathrm{rel} \hat{j}_a \hat{j}_b \hat{j}_c \hat{j}_d (-1)^{\ell + S + J_\mathrm{rel} + L} (-1)^{\ell' + S' + J'_\mathrm{rel} + L'} \langle n N [ (\ell L)\Lambda S] J| n_a n_b [ (\ell_a \ell_b)\Lambda (\tfrac{1}{2}\tfrac{1}{2})S ]J \rangle_{d=1} \langle n' N' [ (\ell' L')\Lambda' S'] J| n_c n_d [ (\ell_c \ell_d)\Lambda' (\tfrac{1}{2}\tfrac{1}{2})S' ]J \rangle_{d=1} \left\{ \begin{matrix} \ell_a & \ell_b & \Lambda \\ 1/2 & 1/2 & S \\ j_a & j_b & J \end{matrix} \right\} \left\{ \begin{matrix} \ell_c & \ell_d & \Lambda' \\ 1/2 & 1/2 & S' \\ j_c & j_d & J \end{matrix} \right\} \left\{ \begin{matrix} L & \ell & \Lambda \\ S & J & J_\mathrm{rel} \end{matrix} \right\} \left\{ \begin{matrix} L' & \ell' & \Lambda' \\ S' & J & J'_\mathrm{rel} \end{matrix} \right\} \langle n\ell S J_\mathrm{rel} T|V_\mathrm{NN}|n'\ell' S' J'_\mathrm{rel} T\rangle\]

wigner9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)

calculate Wigner's 9j symbols, all j should be given as integer (0,1,...) or halfinteger (3/2, 3//2,...)

<τ2・τ3> = 6 (-1)^(t+t'+T+1/2) *wigner6j(t,t',1,1/2,1/2,1/2) * wigner6j(t,t',1,1/2,1/2,T)

anti_op_isospin(params,n12,l12,s12,j12,t12,n3,l3,j3,n45,l45,s45,j45,t45,n6,l6,j6,jtot,ttot,to)

Function to calc. matrix element of antisymmetrizer. Detailed derivation can be found in e.g., Eq.(3.119) of Master Thesis by Joachim Langhammer (2010), TUDarmstadt.

EGM => Eq.(19) in E.Epelbaum et al., PRC 66, 064001(2002). Navratil => Eq.(11) in P.Navratil Few Body Syst. (2007) 41:117-140 Lambda in given in MeV

prepare R_nl(r) with 2n12+l12 + 2N + L = e12 + e3 = N3 <= N3max Note that twice of reduced mass is used in this function

function used for proposals in Affine invariant MCMC

cand:: candidate point given by LatinHypercubeSampling observed:: list of index of cand which has been observed unobserved:: rest candidate indices history:: array of [logprior,logllh,logpost] for i=1:num_cand

for GP

Ktt,Kttinv,Ktp,L:: matrix needed for GP calculation yt:: mean value of training point, to be mean of initial random point yscale:: mean&std of yt that is used to scale/rescale data acquis:: vector of acquisition function values pKernel:: hypara for GP kernel, first one is tau and the other ones are correlation lengths adhoc=> tau =1.0, l=1/domain size

sample_AffineInvMCMC(numwalkers::Int, x0::Matrix{Float64},nstep::Integer, thinning::Integer,a::Float64=2.)

Function to carry out a parameter sampling with Affince invariant MCMC method.

Reference: Goodman & Weare, "Ensemble samplers with affine invariance", Communications in Applied Mathematics and Computational Science, DOI: 10.2140/camcos.2010.5.65, 2010.

RKstep(T,Ho,eta,R,faceta,fRK,Ht)

wrapper function to calc. a Runge-Kutta (RK) step

SRG(xr,wr,V12mom,dict_pwch,to)

Similarity Renormalization Group (SRG) transformation of NN interaction in CM-rel momentum space.

commutator(A,B,R,fac)

wrapper function to overwrite R by fac*(AB-BA)

srg_RK4(Ho,T,Ht,Hs,eta,R,sSRG,face,ds,numit,to; r_err=1.e-8,a_err=1.e-8,tol=1.e-6)

to carry out SRG transformation with RK4

prep_QWs(chiEFTobj,xr,ts,ws)

returns struct QLs, second kind Legendre functions, Eqs.(B1)-(B5) in [Kohno2013]. Note that QLs.QLdict is also used in OPEP to avoid redundant calculations. Reference: [Kohno2013] M.Kohno Phys. Rev. C 88, 064005(2013)

prep_integrals_for2n3n(chiEFTobj,xr,ts,ws,to)

preparing integrals and Wigner symbols for density-dependent 3NFs:

• Fis::Fis_2n3n vectors {F0s,F1s,F2s,F3s}, Eqs.(A13)-(A16) in [Kohno2013].
• QWs:: second kind Legendre functions, Eqs.(B1)-(B5) in [Kohno2013].
• wsyms:: Wingner symbols with specific j used for both 2n and 2n3n.

Reference: [Kohno2013] M.Kohno Phys. Rev. C 88, 064005(2013)

dict_em500n3lo

returns vector&dict. for EM500N3LO: Entem-Machleidt interaction upto N3LO with 500 MeV cutoff

dict_emn500n3lo

returns vector&dict. for EMN500N3LO: Entem-Machleidt-Nosyk interaction upto N3LO with 500 MeV cutoff

dict_emn500n4lo()

returns vector&dict. for EMN500N4LO: Entem-Machleidt-Nosyk interaction upto N4LO with 500 MeV cutoff

prep_valsidxs_dLECs!(vals,idxs,dLECs)

overwrites vals and idxs, which are to be used in e.g., MCMC algorithms, by dict for LECs dLECs.