Hartreefock

struct Dict1b

Fields

• snt2ms::Dict{Int64,Int64} map from snt idx to modelspace(ms) idx
• ms2snt::Dict{Int64,Int64} map from ms idx to snt idxdef_struct.jl

struct HamiltonianNormalOrdered

Fields

• H::Operator Hamiltonian operator
• E0::Float64 NO0B of H
• EMP2::Float64 PT2 correction to E0
• EMP3::Float64 PT3 correction to E0
• Cp::Matrix{Float64} eigenvectors of hp, used unitary trans. HO=>HF basis (proton part)
• Cn::Matrix{Float64} eigenvectors of hn, unitary trans. HO=>HF basis (neutron part)
• e1b_p::Vector{Float64} eigenvalues of hp
• e1b_n::Vector{Float64} eigenvalues of hn
• modelspace::ModelSpace

mutable struct IMSRGobject

Fields

• H0::Operator Hamiltonian for starting point of BCH product
• H::Operator Hamiltonian $H(s)$
• s::Vector{Float} current $s$ and $ds$
• smax::Float maximum $s$
• dsmax::Float maximum $ds$
• maxnormOmega::Float maximum ||Omega||
• eta::Operator generator of IMSRG flow (antihermite Operator)
• Omega::Operator generator of IMSRG flow (antihermite Operator)
• eta_criterion::Float ||eta|| to check convergence
• denominatorDelta::Float64 parameter for multi-major shell decoupling
• n_written_omega::Int # of written Omega by splitting to solve IMSRGflow
• Ncomm::Vector{Int} # of commutator evaluated during IMSRG flow

struct ModelSpace

Fields

• p_sps::Vector{SingleParticleState} proton single particle states (only odd integer)
• n_sps::Vector{SingleParticleState} neutron single particle states
• sps::Vector{SingleParticleState} single particle states (odd number ones=>proton even=>neutron)
• occ_p::Matrix{Float64} matrix representing the occupation number of proton (needed for density matrix)
• occ_n::Matrix{Float64} matrix representing the occupation number of neutron
• holes::Vector{Vector{Int64}} idx list of holes
• particles::Vector{Vector{Int64}} idx list of particles
• spaces::space_channel space_channel (mutable struct)

mutable struct Operator

Fields

• zerobody::Vector{Float64} zerobody part of the operator
• onebody::Vector{Matrix{Float64}} one-body matrix elements ([1]=>proton, [2]=>neutron)
• twobody::Vector{Matrix{Float64}} two-body matrix elements, having array structure [ch]
• hermite::Bool whether it is hermitian operator or not
• antihermite::Bool antihermitian or not

struct PandyaObject, used for Pandya transformation (especially in comm222ph_ss!)

mutable struct SingleParticleState

Fields

• n::Int64 principal quantum number of the single particle state(sps)
• l::Int64 azimuthal quantum number of the sps
• j::Int64 angular momentum
• tz::Int64 z-component of isospin (doubled) tz=-1 => proton & tz=1 => neutron
• occ::Float64 occupation number (can be fractional) of the sps
• c::Bool indicating whether the single-particle state belongs to "core" or not
• v::Bool whether belongs to "valence" or not
• q::Bool whether belongs to "q-space" or not

struct VdictCh

Fields

• Vch::Int64 two-body channel (specified by JPT)
• Vdict::Dict{Int64,Int64} dict to get idx from ket, which is used in only vPandya function for HFMBPT

struct basedat contains base infomation of the calculation

Fields

• nuc::nuclei information of target/core nucleus
• sntf::String filename/path to input interaction
• hw::Int64 hbar omega parameter used to define single particle states
• emax::Int64 emax truncation for the entire calculations
• ref::String to specify ref="core" or ref="nucl"

struct chan1b

Fields

• chs1b::Vector{Dict{Int64,Vector{Int64}}} dict of single particle states with non-zero contribution (having same l,j) [dict for proton sps, dict for neutron sps]
• chs1b_redundant::Vector{Dict{Int64,Vector{Int64}}} redundant version of chs1b (with i>j case)
• snt2ms::Dict{Int64,Int64} map from snt idx to modelspace(ms) idx
• ms2snt::Dict{Int64,Int64} map from ms idx to snt idx

struct chan2b referred to as "tbc" (two-body channel) in some functions

Fields

• Tz::Int64 total tz, -2(pp),0(pn),2(n)
• prty::Int64 parity
• J::Int64 total J
• kets::Vector{Vector{Int64}} vector of ket (e.g. [1,1], [1,3],...)

struct Chan2bD

Fields

• Chan2b::Vector{chan2b} array of chan2b (ch=1,...,nchan)
• dict_ch_JPT::Dict{Vector{Int64},VdictCh} dict to get VdictCh by given key [J,prty,T]
• dict_ch_idx_from_ket::Vector{Vector{Dict{Int64,Vector{Vector{Int64}}}}} dict to get [ch,idx], having array structure [pnrank(=1/2/3)][J+1], key=ket
• dict_idx_from_chket::Vector{Dict{Vector{Int64},Int64}} dict to get idx from ket, having array structure [ch]

struct dWS2n, Wigner symbols used in PreCalcHOB

Fields

• dtri::Dict{Vector{Int64},Float64} dict for trinomial
• dcgm0::Dict{Int64,Float64} dict for special CG coefficients (l0l'0|L0)

struct dictTBMEs contains dictionaries for TBME/monopole

Fields

• dictTBMEs::Vector{Dict{Vector{Int64},Float64}} one can get pp/pn/nn dict by dictTBMEs[pnrank] (pnrank=1,2,3)
• dictMonopole::Vector{Dict{Vector{Int64},valDictMonopole}} one can get monopole component of two-body interaction by dictMonopole[pnrank][key], key to be ket array like [1,1]

struct hfdata, used to calculate multiple nuclei in a single runscript

Fields

• nuc::nuclei information of target/core nucleus
• data::Vector{Vector{Float64}} will be experimental data from AME2020 (if available)
• datatype::Vector{String} supposed to be ["E"] for now

struct nuclei

Fields

• Z::Int64 proton number of the reference nucleus
• N::Int64 neutron number of the ref.
• A::Int64 mass number of the ref.
• el::String element (e.g., "He")
• cnuc::String string element nameA (e.g., "He8")
• cZ::Int64 proton number of core nucleus
• cN::Int64 neutron number of core
• corenuc::String core nucleus (e.g., "He4")

struct single_util122, used to make operation related commutator122 matrix manipulation

mutable struct space_channel; dictionaries to get the two-body channels that have kets (specified by pp,ph, etc.)

Fields

• pp::Dict{Int64,Vector{Int64}} particle-particle
• ph::Dict{Int64,Vector{Int64}} particle-hole
• hh::Dict{Int64,Vector{Int64}} hole-hole
• cc::Dict{Int64,Vector{Int64}} core-core
• vc::Dict{Int64,Vector{Int64}} valence-core
• qc::Dict{Int64,Vector{Int64}} qspace-core
• vv::Dict{Int64,Vector{Int64}} valence-valence
• qv::Dict{Int64,Vector{Int64}} qspace-valence
• qq::Dict{Int64,Vector{Int64}} qspace-qspace

HF_MBPT2(binfo,modelspace,fp,fn,e1b_p,e1b_n,Chan2b,Gamma)

Calculate 2nd order correction to HF energy

\[E^{(2)} = \frac{1}{4}\sum_{abij} \frac{\bar{H}^{[2]}_{abij} \bar{H}^{[2]}_{ijab}}{\epsilon^{ab}_{ij}} = \frac{1}{4} \sum_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}}\sum_{\{m\}}\sum_{JJ'MM'} \frac{{}^J\bar{H}^{[2]}_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}} {}^J\bar{H}^{[2]}_{\tilde{i}\tilde{j}\tilde{a}\tilde{b}}}{\epsilon^{ab}_{ij}} (j_a j_b m_a m_b|J M) (j_a j_b m_a m_b|J' M') (j_i j_j m_i m_j|J M) (j_i j_j m_i m_j|J' M')\]

\[= \frac{1}{4} \sum_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}}\sum_{JJ'MM'} \frac{{}^J\bar{H}^{[2]}_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}} {}^{J'}\bar{H}^{[2]}_{\tilde{i}\tilde{j}\tilde{a}\tilde{b}}}{\epsilon^{ab}_{ij}} \delta_{JJ'} \delta_{MM'} = \frac{1}{4} \sum_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}}\sum_{J}(2J+1) \frac{{}^J\bar{H}^{[2]}_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}} {}^J\bar{H}^{[2]}_{\tilde{i}\tilde{j}\tilde{a}\tilde{b}}}{\epsilon^{ab}_{ij}}\]

HF_MBPT3(binfo,modelspace,e1b_p,e1b_n,Chan2b,dict_2b_ch,dWS,Gamma,to;io=stdout)

Calculate 2nd order correction to HF energy

\[E^{(3)}= \frac{1}{8} \sum_{\tilde{a}\tilde{b}\tilde{c}\tilde{i}\tilde{j}\tilde{k}} \sum_{J}(2J+1) \frac{ {}^{J}\bar{H}^{[2]}_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}} {}^{J}\bar{H}^{[2]}_{\tilde{i}\tilde{j}\tilde{c}\tilde{d}} {}^{J}\bar{H}^{[2]}_{\tilde{c}\tilde{d}\tilde{a}\tilde{b}} } {\epsilon^{\tilde{a}\tilde{b}}_{\tilde{i}\tilde{j}}\epsilon^{\tilde{c}\tilde{d}}_{\tilde{i}\tilde{j}}} + \frac{1}{8} \sum_{\tilde{a}\tilde{b}\tilde{i}\tilde{j}\tilde{k}\tilde{l}} \sum_{J}(2J+1) \frac{ {}^{J}\bar{H}^{[2]}_{\tilde{i}\tilde{j}\tilde{a}\tilde{b}} {}^{J}\bar{H}^{[2]}_{\tilde{a}\tilde{b}\tilde{k}\tilde{l}} {}^{J}\bar{H}^{[2]}_{\tilde{k}\tilde{l}\tilde{i}\tilde{j}} } {\epsilon^{\tilde{a}\tilde{b}}_{\tilde{i}\tilde{j}}\epsilon^{\tilde{c}\tilde{d}}_{\tilde{i}\tilde{j}}} -\sum_{\tilde{a}\tilde{b}\tilde{c}\tilde{i}\tilde{j}\tilde{k}} \sum_{J}(2J+1) \frac{ {}^JH^{XC}_{\tilde{a}\tilde{i}\tilde{j}\tilde{b}} {}^JH^{XC}_{\tilde{j}\tilde{b}\tilde{k}\tilde{c}} {}^JH^{XC}_{\tilde{k}\tilde{c}\tilde{a}\tilde{i}} }{ \epsilon^{\tilde{a}\tilde{b}} \epsilon_{\tilde{i}\tilde{j}} \epsilon^{\tilde{a}\tilde{c}} \epsilon_{\tilde{k}\tilde{j}} }\]

Ref. Many-Body Methods in Chemistry and Physics by Isaiah Shavitt and Rodney J. Bartlett (2009, Cambridge Molecular Science). More details can be found in e.g. Dr. thesis by A.Tichai (2017, TU Darmstadt).

vPandya(a,b,c,d,ja,jb,jc,jd,totJ,dict_2b_ch,d6j_lj,Gamma,keych,key6j,keyab;verbose=false)

returns generalized Pandya transformed matrix element:

\[\tilde{V}^J_{ajib} = -\sum_{J'} [J'] \begin{Bmatrix} j_a & j_j & J \\ j_i & j_d & J' \end{Bmatrix} V^{J'}_{abij}\]

cZN_from_corenuc(rZ,rN,corenuc)

get $Z$ and $N$ of the core nucleus

def_chan1b(dim1b,sps,dicts1b)

define Chan1b: dict. to get 1b-channels to be coupled to a given channel Chan1b = [ dictforprotonsps, dictforneutronsps] dicts1b Basically ToBeCoupled will be used, but redundant one is needed in some cases

def_chan2b(binfo,dicts,sps)

define two-body utils and returns them as Chan2bD struct

def_nuc(Z,N,ref,corenuc)

constructor of nuclei strict from given cnuc,ref,corenuc

def_nuc(nuc::Vector{Int},ref,corenuc)

constructor of nuclei strict from given Z,N,ref,corenuc

get_ZNref(ref,Z,N,corenuc)

get $Z$ and $N$ of the target reference

get_abcdarr_from_intkey!(nkey,abcdarr;ofst=1000)

A destructive function to get original array from integer array.

make_sps_and_dict_isnt2ims(p_sps,n_sps,lp)

make dicts1b, snt-idx(sntidx) = 1-lp (proton) & lp+1~lp+ln (neutron), modelspace-idx(msidx) = odd(1,3,...)-> proton, even(2,4,...) -> neutron

returns:

• dict_snt2ms: from sntidx to msidx
• dict_ms2snt: from msidx to sntidx

make_sps_from_pnsps(p_sps,n_sps,Chan1b)

construct sps from p_sps and n_sps

readsnt(sntf,binfo,to)

Function to read snt file. Note that it is slightly different from readsnt() in ShellModel.jl.

readsnt(sntf,binfo,to)

Function to read snt.bin file.

rm_comment(lines)

remove fortran like comment from input (snt fmt) strings

store_1b2b(sps,dicts1b,dicts,binfo)

update_1b!(binfo,sps,Hamil)

Update one-body(1b) part of Hamiltonian for different target nuclei

update_2b!(binfo,sps,Hamil,dictTBMEs,Chan2bD,dicts)

Update two-body(2b) kinetic part for different target nuclei

ReorderHFSPS!(h_p,h_n,Cp,Cn,e1b_p,e1b_n,Chan1b)

"reorder" HF single particle space. Since we diagonalize the h_p,h_n (istead of subblock mat), we need to specify the correspondance between ordering of sps and that of HFSPEs obtained by solving HF eigenvalue problem

add_ch_ket!(ch,iket,tdict)

add ch & idx for kets in spaces::space_channel (pp/hh/etc.)

calc_Gamma!(Gamma,sps,Cp,Cn,V2,Chan2b,maxnpq)

calculate $\Gamma$ (two-body HF interaction)

def_holeparticle(Chan1b,occ_p,occ_n,p_sps,n_sps)

define hole/particle space by occ_p, occ_n

def_sps(emax)

Function to define sps::Vector{SingleParticleState} from emax.

getHNO(binfo,tHFdata,E0,p_sps,n_sps,occ_p,occ_n,h_p,h_n,e1b_p,e1b_n,Cp,Cn,V2,Chan1b,Chan2b::tChan2b,Gamma,maxnpq,dict_2b_ch,dWS,to) where{tChan2b <: Vector{chan2b}}

obtain spherical HF solution and calc. MBPT correction (upto 2nd&3rd order) to g.s. energy

get_space_chs(sps,Chan2b)

define hole/particle single particle states. In this function, only the hh/pp/ph (needed for IMSRG) are defined, and other channels will be updated later for target normal ordering or VS-IMSRG flow.

hf_iteration(binfo,tHFdata,sps,Hamil,dictTBMEs,Chan1b,Chan2bD,Gamma,maxnpq,dWS,to;itnum=100,verbose=false,HFtol=1.e-14,inttype="snt")

solve HF equation

This function returns object with HamiltonianNormalOrdered (HNO) struct type, which contains...

• E0,EMP2,EMP3 HF energy and its MBPT corrections
• fp/fn::Matrix{Float64} one-body int.
• Gamma:: Vector{Matrix{Float64}} two-body int.

hf_main(nucs,sntf,hw,emax;verbose=false,Operators=String[],is_show=false,doIMSRG=false,valencespace=[],corenuc="",ref="nucl")

main function to carry out HF/HFMBPT calculation from snt file

Arguments

• nucs::Vector{String} target nuclei
• sntf path to input interaction file
• hw hbar omega
• emax_calc emax for HF/IMSRG

Optional Arguments

• verbose=false to see detailed stdout for HF
• Operators=String[] target observables other than Hamiltonian
• is_show=false to show TimerOutput log (summary of run time and memory allocation)
• doIMSRG=false to carry out IMSRG/VSIMSRG calculation
• valencespace=[] to spacify the valence space (e.g., ["sd-shell"]), if this is not empty, it tries to do vsimsrg calculations
• corenuc="" core nucleus, example=> "He4"
• ref="nucl" to specify target reference state, "core" or "nucl" is supported

hf_main_mem(chiEFTobj,nucs,dict_TM,dWS,to;verbose=false,Operators=String[],valencespace=[],corenuc="",ref="core")

"without I/O" version of hf_main

ini_occ!(pconf,occ_p,nconf,occ_n)

initialize occupation number matrices (occ_p&occ_n) by naive filling configurations pconf&nconf

naive_filling(sps,n_target,emax,for_ref=false)

calculate naive filling configurations by given sps and proton/neutron number (n_target)

For some nuclei, especially for heavier ones, carrying out naive filling is ambiguous (e.g., neutron occupation of 22O can be both 0s1(2),0p1(2),0p3(4),0d5(6) and 0s1(2),0p1(2),0p3(4),1s1(2), 0d3(4)). In this function, "naive filling" means to try fill orbits with lower $2n+l$ and then "lower" $j$.

prepHFdata(nucs,ref,datatype,corenuc)

Constructor of an array of hfdata struct.

printEHF(Es)

print HF energy and its break down Es=[E1b,E2bpp,E2bnn,E2bpn]

recalc_v!(A,dicts)

Function to calculate two-body interaction from snt file. This is needed because in the readsnt/readsnt_bin function, the interaction part and the kinetic term are stored separately to avoid multiple reads of the input file for calculation of multiple nuclei.

Calc_Expec(binfo,Chan1b,Chan2b,HFobj,Op_Rp2,dict_2b_ch,dWS,MatOp,to;hfmbptlevel=true,verbose=false)

Calculate expectation value of Rp2 and its HFMBPT corrections.

Details about HFMBPT correction can be found in Many-Body Methods in Chemistry and Physics by Isaiah Shavitt and Rodney J. Bartlett (2009, Cambridge Molecular Science) or Appendix in T. Miyagi et al., Phys. Rev. C 105, 0143022 (2022).

TCM/A

One doesn't have to calc. two-body part here, since the necessary part is already calculated in Vpp.

Calculate_RCM(binfo,Chan1b,Chan2b,sps,Op_Rp2,dWS,to;non0_cm=true,non0_ij=true)

calculate $R_{CM}$ term

Note that rirj term is also included here to avoid redundancy.

Calculate_Rp(binfo,Chan1b,Chan2b,HFobj,Op_Rp2,dWS,dict_2b_ch,to;hfmbptlevel=true)

To calculate squared point proton radius and its MBPT correction. The squared point proton radius is related to the charge radius as follows

\[R^2_{ch} = R^2_p + \langle r^2_p \rangle + \frac{N}{Z} \langle r^2_n \rangle + \frac{3}{4m^2_p c^4} + \langle r^2 \rangle_{SO},\]

where $\langle r^2_p \rangle = 0.769 \mathrm{fm}^2$, $\langle r^2_n \rangle = -0.116 \mathrm{fm}^2$, $\frac{3}{4m^2_p c^4} =0.033\mathrm{fm}^2$ is the so-called Darwin-Foldy term, and the last term is Spin-Orbit correction term.

Calculate_SOterm(binfo,Chan1b,HFobj,Op_Rp2)

Calculate Spin-Orbit Correction term for Rp2. We are using "simple expression for the correction to the mean-square charge radius due to the spin-orbit term" in the reference below:

\[\langle r^2 \rangle_{SO} =-\frac{1}{Z}\sum_i \frac{\mu_i}{M^2} (\kappa_i+1),\]

where $\mu_p = 2.793 \mu_N$, $\mu_n = −1.913\mu_N$, and $\kappa = \ell (\mathrm{for } j=\ell-1/2), -(\ell+1) (\mathrm{for} j=\ell+1/2)$

Ref: A.Ong, J.C.Berengut, and V.V.Flambaum, Phys. Rev. C 82, 014320 (2010).

Calculate_intR2p(binfo,Chan1b,HFobj,Op_Rp2)

calculate a part of squared point proton radius R2p. Note that this function overwrites the onebody part of given Op_Rp2.

InitOp(Chan1b,Chan2b)

initialize scaler Operator, note that hermite is true as default

aOp!(Op::Operator,a::Float64)

function to multiply scaler to an operator.

aOp1_p_bOp2!(Op1::Operator,Op2::Operator,a::Float64,b::Float64)

function to overwrite Op2 by a*Op1 + b*Op2

aOp1_p_bOp2(Op1::Operator,Op2::Operator,a,b)

function to return new Operator = aOp1 + bOp2

aOp1_p_bOp2_Op3!(Op1::Operator,Op2::Operator,Op3::Operator,a,b,c)

function to overwrite Op3 by c*Op3 + a*Op1 + b*Op2

calc_single_r1r2(bra,ket,sps,J,dWS,b2,to)

Calc $<r_1 \cdot r_2>$ for a given 2b-channel.

• bra: <ab| a&b: s.p.s. (n,l,j,tz)
• ket: |cd> c&d: s.p.s. (n,l,j,tz)

function to redefine RCM. Note that non0_ij for Calculate_RCM assumed to be false.

eval_rch_hfmbpt(binfo,Chan1b,Chan2bD,HFobj,Op_Rp2,dWS,to)

function to evaluate charge radii with HFMBPT

eval_rch_imsrg(binfo,Chan1b,Chan2bD,HFobj,IMSRGobj,PandyaObj,dWS,dictMono,to)

evaluate charge radii with IMSRG.

getNormalOrderedO(HFobj,targetOp,Chan1b,Chan2bD,to;verbose=false,undo=false,OpeqH=false,firstNO=false)

NormalOrdering for a target Operator. For now, it only supports scaler operators.

show_Hamil_norm(Op::Operator;tol=1.e-6,normtype="fro")

Function to show 1b/2b norm of a given Operator. It may be usuful for debug.