ShellModel

Struct to store 2b-jumps (T=1).

This includes the following information: f::Int64 the final index f for a given index i/ coef::Float64 the coefficient of the 2b-jump.

all_perm!(ln::Int64,num_valence::Int64,occs::Array{Array{Bool,1}})

make all possible permutation of 'bits'

Example: If 2 protons and 1 neutron are in the 0p-shell space, valence orbits(0p1/2,0p3/2) => -1/2, 1/2, -3/2, -1/2, 1/2, 3/2

configurations are represented like:

proton: 000011, 000101, ..., 110000

neutron: 000001, 000010, ..., 100000

Function to check if the specified configuration tocc_j is valid within a given truncation scheme.

Function to check occupations in proton and neutron are only in time-reversal pairs.

constructmsps(psps,n_sps)

Function to define the single particle states specified by [n,l,j,tz,mz,p(n)idx]. The last elements pidx and nidx to represent original index of p_sps/n_sps, which is specified by [n,l,j,tz]`.

mainsm(sntf,targetnuc,numev,targetJ; savewav=false,q=1,isblock=false,isshow=false,numhistory=3,lm=300,ls=20,tol=1.e-8, inwf="",mdimmode=false,calcmoment=false, visualize_occ=false, gfactors=[1.0,0.0,5.586,-3.826],effcharge=[1.5,0.5])

Digonalize the model-space Hamiltonian

Arguments:

• sntf: path to input interaction file (.snt fmt)
• target_nuc: target nucleus
• num_ev: number of eigenstates to be evaluated
• target_J: target total J (specified by e.g. [0]). Set to [] if you want lowest states with any J. Note that J should be doubled (J=0=>[0], J=1/2=>[1], J=1=>[2],...)

Optional arguments:

• q=1 block size for Block-Lanczos methods
• is_block=false whether or not to use Block algorithm
• save_wav=false whether or not to save wavefunction file
• is_show = true to show elapsed time & allocations
• lm = 100 number of Lanczos vectors to store
• ls = 20 number of vectors to be used for Thick-Restart
• tol= 1.e-8 tolerance for convergence check in the Lanczos method
• in_wf="" path to initial w.f. (for preprocessing)
• mdimmode=false true => calculate only the M-scheme dimension
• calc_moment=false true => calculate mu&Q moments
• calc_entropy=false true => calculate entropy, which is not optimized yet.
• visualize_occ=false true => visualize all configurations to be considered
• gfactors=[1.0,0.0,5.586,-3.826] angular momentum and spin g-factors
• effcgarge=[1.5,0.5] effective charges
• truncation_scheme=="" option to specify a truncation scheme, "jocc" and "pn-pair" is supported and you can use multiple schemes by separating them with camma like "jocc,pn-pair".
• truncated_jocc=Dict{String,Vector{Int64}}() option to specify the truncation scheme for each orbit, e.g. Dict("p0p1"=>[1],"n0p3"=>[2,3]) means that the occupation number for proton 0p1 is truncated upto 1, and for neutron 0p3 min=2 and max=3"

occ(p_sps::Array{Array{Int64,1}},msps_p::Array{Array{Int64,1}},mzp::Array{Int64,1},num_vp::Int64,
    n_sps::Array{Array{Int64,1}},msps_n::Array{Array{Int64,1}},mzn::Array{Int64,1},num_vn::Int64,Mtot::Int64;pnpair_truncated=false)

prepare bit representations of proton/neutron Slater determinants => pbits/nbits

joccp, joccn: corresponding occupation numbers for a "j" state, which is used for one-body operation and OBTDs.

Mps/Mns: total M for proton/neutron "blocks"

For 6Li in the p shell and M=0, Mps = [-3,-1,1,3] & Mns = [3,1,-1,-3] blocks => [ (Mp,Mn)=(-3,3),(Mp,Mn)=(-1,1),...]

tdims: array of cumulative number of M-scheme dimensions for "blocks"

tdims =[ # of possible configurations of (-3,3), # of possible configurations of (-1,1),...]

function prep_Hamil_T1(p_sps::Array{Array{Int64,1}},n_sps::Array{Array{Int64,1}},
            msps_p::Array{Array{Int64,1},1},msps_n::Array{Array{Int64,1},1},
            labels::Array{Array{Array{Int64,1},1},1},TBMEs::Array{Array{Float64,1}})

First, this makes bit representations of T=1 (proton-proton&neutron-neutron) interactions for each {m_z}.

prep_Hamil_pn(p_sps::Array{Array{Int64,1}},n_sps::Array{Array{Int64,1}},
       msps_p::Array{Array{Int64,1},1},msps_n::Array{Array{Int64,1},1},
       labels::Array{Array{Int64,1}},TBMEs::Array{Float64,1},zeroME=false)

make bit representation of T=0 (proton-neutron) interactions for each {m_z}

readsmsnt(sntf,Anum)

To read interaction file in ".snt" format.

• sntf: path to the interaction file
• Anum: mass number (used for "scaling" of TBMEs)

visualize all configurations in a given model-space

Thick-Restart Block Lanczos method

Set HO parameter by an empirical formula.

The default choice is 41A^(-1/3) MeV. For sd-shell nuclei (16<=A<=40), we use 45A^(-1/3)-25A^(-2/3) MeV by J. Blomqvist and A. Molinari, Nucl. Phys. A106, 545 (1968).

transit_main(sntf,target_nuc,jl2,jr2,in_wfs;num_ev_l=100,num_ev_r=100,q=1,is_block=false,is_show=true,calc_EM=true,gfactors=[1.0,0.0,5.586,-3.826],eff_charge=[1.5,0.5])

Function tocalculate the M1&E2 transitions from two given wavefunctions

Arguments

• sntf:path to the interaction file
• target_nuc:target nucleus in string e.g., "Si28"
• jl2:J*2 for the left w.f.
• jr2:J*2 for the right w.f.
• in_wfs:["path to left wf","path to right wf"]

Optional arguments

• num_ev_l=100:upper limit of the number of eigenvectors for the left w.f.
• num_ev_r=100:upper limit of the number of eigenvectors for the right w.f.
• is_show=true:to display the TimerOutput
• gfactors=[1.0,0.0,5.586,-3.826]:g factors [glp,gln,gsp,gsn]
• eff_charge=[1.5,0.5]:effective charges [ep,en]

Optional arguments (not used, but for future developments)

• q=1:block size
• is_block=false:to use Block algorithm

solveEC(Hs,target_nuc,tJNs;
        write_appwav=false,verbose=false,calc_moment=true,wpath="./",is_show=false,
        gfactors = [1.0,0.0,5.586,-3.826],effcharge=[1.5,0.5],exact_logf="")

To solve EC (generalized eigenvalue problem) to approximate the eigenpairs for a given interaction.

\[H \vec{v} = \lambda N \vec{v}\]

Transition densities and overlap matrix for H and N are read from "tdmat/" directory (To be changed to more flexible)

Arguments:

• Hs:array of paths to interaction files (.snt fmt)
• target_nuc: target nucleus
• tJNs:array of target total J (doubled) and number of eigenstates to be evaluated e.g., [ [0,5], [2,5] ], when you want five lowest J=0 and J=1 states.

Optional arguments:

• write_appwav=false:write out the approximate wavefunction
• verbose=false:to print (stdout) approx. energies for each interaction
• calc_moment=true: to calculate mu&Q moments
• wpath="./": path to sample eigenvectors to construct approx. w.f.
• is_show=false: to show TimerOutput
• gfactors=[1.0,0.0,5.586,-3.826]: g-factors to evaluate moments
• effcharge=[1.5,0.5]:effective charges to evaluate moments

Optional arguments for author's own purpose

• exact_logf="":path to logfile for E(EC) vs E(Exact) plot

read_kshell_summary(fns::Vector{String};targetJpi="",nuc="")

Reading KSHELL summary file from specified paths, fns. In some beta-decay studies, multiple candidates for parent g.s. will be considered. In that case, please specify optional argument targetJpi e.g., targetJpi="Si36", otherwise the state havbing the lowest energy in summary file(s) will be regarded as the parent ground state. Example:

Energy levels

N    J prty N_Jp    T     E(MeV)  Ex(MeV)  log-file

1     0 +     1     6   -319.906    0.000  log_Ar48_SDPFSDG_j0p.txt 

Ecoulomb_SM(Z,N)

Additional Coulomb energy (MeV)

\[E_C(Z,N) = 0.7 \frac{Z(Z-1) -0.76 [Z(Z-1)]^{2/3}}{R_c} \\ R_c = e^{1.5/A} A^{1/3} \left( 0.946 -0.573 \left( \frac{2T}{A} \right)^2 \right) \\ T = |Z-N|/2 \]

used in references.

• S.Yoshida et al., Phys. Rev. C 97, 054321 (2018).
• J. Duflo and A. P. Zuker, Phys. Rev. C 52, R23 (1995).
• E. Caurier et al., Phys. Rev. C 59, 2033 (1999).

Fermi_function(Z,We,R)

\[F(Z,W) = 4 (2p_e R)^{-2(1-\gamma_1)} \exp \left( \pi y \frac{|\Gamma(\gamma_1 + iy)|^2}{ |\Gamma(2\gamma_1 + 1)|^2} \right) \\ \gamma_k = \sqrt{k^2 - \alpha^2 Z^2} \\ y = \alpha Z W / p_e\]

Fermi_integral(Qval,Ex,pZ,A,nmesh=40)

Fermi integrals are evaluated with Eqs. in Chapter 5 of "Weak Interactions and Nuclear Beta Decay" by H. F. Schopper.]

\[f_0 = \int^{W_0}_{1} F(Z,W) \sqrt{W^2-1} W(W_0-W)^2 dW, \\ W_0 = Q(\beta^-) + 1 - E_\mathrm{ex.}\]

Note that Ws are in natural unit, divided by $m_ec^2$.

ME_FF_func!(chan,qfactors,C_J,Mred,lambda_Ce,mass_n_nu,dict)

Function to calc FF matrix elements. Since M0rs,M1r,M1rs,M2rs (M0sp,M1p) are evaluated in fm (1/fm) in KSHELL, it is needed to multiply 1/lambdaCe (lambdaCe) to get transition matrix element in natural unit.

eval_betadecay_from_kshell_log(fns_sum_parent::Vector{String},fns_sum_daughter::Vector{String},fns_GT::Vector{String},fns_FF::Vector{String},parentJpi::String;pnuc="",dnuc="",verbose=false)

Arguments

• fns_sum_parent::Vector{String} vector of paths to KSHELL-summary files for the parent nucleus
• fns_sum_daughter::Vector{String} vector of paths to KSHELL-summary files for the daughter nucleus
• fns_GT::Vector{String} vector of paths to KSHELL-log files for Gamow-Teller transitions.
• fns_FF::Vector{String} vector of paths to KSHELL-log files for First-Forbidden transitions.

Optional Arguments

• pnuc::String to specify the parent nuclei explicitly
• dnuc::String to specify the daughter nuclei explicitly
• verbose::Bool option to show GT/FF transitions for each state