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,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, nnlosat.

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

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.

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.

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)

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.

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-Moshinsky transformation for NN interaction in HO space and to write out an sinput file.

Function to get canonical order of 6j-symbol arguments:

\[\left\{ \begin{matrix} j_1 & j_3 & j_5 \\ j_2 & j_4 & j_6\end{matrix} \right\} \]

This may not bring speed up though... The canonical order is defined as follows:

• Since the 6j-symbol is invariant under permutation of any two columns, we can always re-order the columns such that $j_1+j_2 \leq j_3+j_4 \leq j_5+j_6$.
• If sum of two columns are equal, we re-order columns based on j_col_score, using the minimum j value of the two columns.
• Since the 6j-symbol is invariant under permutation rows for any two columns, we can always re-order the rows such that $j_1 \leq j_2, j_3 \leq j_4, j_5 \leq j_6$ if any column have the same j.
• If the elements of any column are different from each other, the relation between upper and lower j will always be either $\vee \vee \vee \& \vee \land\land$ or $\land\land\land \& \land \vee \vee$. We can always re-order the rows such that the $\land\land\land$ or $\vee \vee \vee$ is satisfied.

Unlike other wigner symbols, constructing CG coefficients enconter the case with negative indices. One way to avoid this is to use the following function to get hash key for the CG coefficients. This function asssume all j is >= -3

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.

Function to construct dWS2n struct, CG-coefficients and Wigner symbols for the given parameters.

it doesn't work for now

Function to prepare CG coefficient for recoupling of spin (isospin). The resultant dictionary is used for e.g. to use Three-body matrix elements. Note that the all indices are doubled.

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)

Function to get the 3BME from the vector v3bme. To this end, we need to calculate the index of the 3BME in the vector v3bme from given Jab,Jde,J2,tab,tde,T,a,b,c,d,e,f.

Since the dict for 3BME, the order if set so that a>=b>=c, d>=e>=f, one needs to be careful to get the correct index for the 3BME.

Allocating 3BME vector and dictionary for the indices. Note that the dimension of v3bme is the one that is used for the 3BME (for modelspace) instead of the number of elements in the input ThBME file.

Counting the offset for the 3BME indices.

Since we gonna store in isospin formalism, we use only proton part.

While the snt format is given ascending order for $\ell$ and $j$ the 3BME is stored in descending order for $\ell$ and $j$. Hence, we need to reorder the indices for the 3BME. Note that this is ad hoc and not guaranteed to work for all cases.

read_me3jgz(fn, count_ME_file, to; verbose)

Function to read me3j.gz using GZip. The values are stored in a vector ThBME. The ordering of ThBME is not considered here.

The (re)order is done with only odd(tz=-1) indices. abc=0, bca=1, cab=2, acb=3, bac=4, cba=5

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.