Skip to content

Commit

Permalink
Merge pull request #65 from SotaYoshida/develop
Browse files Browse the repository at this point in the history
Develop
  • Loading branch information
SotaYoshida authored Dec 3, 2022
2 parents 6ca0ef4 + 29a6d9e commit eed67f8
Show file tree
Hide file tree
Showing 8 changed files with 193 additions and 29 deletions.
2 changes: 1 addition & 1 deletion docs/src/contributing.md
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ Function to calculate the lowest ``n`` wavefunctions and overwrite results to `
- `ψ::WaveFunction` the wavefunction, see struct `WaveFunction` for more details.
- `n::Int64` the number of states of interest.
"""
function Show_EigenStates!(H::Hamiltonian::WaveFunction,n::Int64)
function Eval_EigenStates!(H::Hamiltonian::WaveFunction,n::Int64)
# write your nice codes here

end
Expand Down
2 changes: 1 addition & 1 deletion example/log_sample_script.txt
Original file line number Diff line number Diff line change
Expand Up @@ -282,6 +282,6 @@ Starting VS-IMSRG flow
Operator:Rp2
HF point proton radius 2.045093 charge radius 2.206924 => HF+PT 2.219037
IMSRG point proton radius 2.091495 charge radius 2.249990
Mg24 Z,N=(12,12) c(8,8) v(4,4) mdim: 28503 ( 4.45 ) J [0.0, 2.0, 2.0, 3.0, 4.0, 4.0, 5.0, 1.0, 6.0, 2.0]
Mg24 Z,N=(12,12) c(8,8) v(4,4) mdim: 28503 ( 4.45 ) J [0.0, 2.0, 2.0, 3.0, 4.0, 4.0, 5.0, 1.0, -1.0, -1.0]
En. -109.065 -107.717 -106.056 -104.800 -104.748 -103.238 -101.177 -100.704 -100.138 -100.119
Ex. 0.000 1.348 3.010 4.265 4.317 5.828 7.889 8.361 8.928 8.946
4 changes: 3 additions & 1 deletion src/NuclearToolkit.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,12 +57,14 @@ include("ShellModel/input_int_snt.jl")
include("ShellModel/eigenvector_continuation.jl")
include("ShellModel/KSHELL.jl")
include("ShellModel/betadecay.jl")
include("ShellModel/trans_snt_msnt.jl")
export read_smsnt
export main_sm,samplerun_sm # from shellmodel.
export prepEC,solveEC,solveEC_UQ # from eigenvector_continuation.jl
export transit_main # from transit.jl
export read_kshell_summary
export eval_betadecay_from_kshell_log

export main_trans_msnt
end


Expand Down
5 changes: 2 additions & 3 deletions src/ShellModel/betadecay.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@ function get_quenching_factors(qtype="")
if qtype == "" || qtype == "SY"
# S. Yoshida et al., PRC 97, 054321(2018).
return quenching_factors(0.74,1.0,1.7,1.0,1.0,0.55)
#return quenching_factors(0.77,1.0,1.7,1.0,1.0,0.55)
elseif qtype == "W" || qtype == "Warburton"
# E. Warburton, J. Becker, B. Brown, and D. Millener, Ann. Phys. 187, 471 (1988).
return quenching_factors(1.0,1.1,1.5,1.0,1.0,0.51)
Expand Down Expand Up @@ -286,8 +285,8 @@ function read_bgtstrength_file!(fn::String,qfactors::quenching_factors,parent::k
if occursin("parity =",line)
prty = strip(split(line,"parity =")[end])
end
if occursin("N_EIGEN ",line)
n_eigen = parse(Int,split(split(line)[end],",")[1])
if occursin("N_EIGEN",line)
n_eigen = parse(Int,split(split(line,"=")[end],",")[1])
end
if !occursin("strength function ",line); continue; end
tl = split(line)
Expand Down
11 changes: 6 additions & 5 deletions src/ShellModel/lanczos_methods.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ function TRL(vks,uks,Tmat,k,
SPEs,ppinfo,nninfo,bis,bfs,block_tasks,
p_NiNfs,n_NiNfs,Mps,delMs,Vpn,
eval_jj,oPP,oNN,oPNu,oPNd,Jidxs,
tdims,num_ev,num_history,lm,ls,en,tol,to,doubleLanczos=false)
tdims,num_ev,num_history,lm,ls,en,tol,to,doubleLanczos=false,checkorth=false)
mdim = tdims[end]
TF=[false]

Expand Down Expand Up @@ -34,13 +34,14 @@ function TRL(vks,uks,Tmat,k,
end
talpha = dot(vk,vkp1)
Tmat[it,it] = talpha
diagonalize_T!(it,num_ev,Tmat,en,num_history,TF,tol)
diagonalize_T!(it,num_ev,Tmat,en,num_history,TF,tol)
if it == mdim; TF[1]=true;end
if TF[1];elit=it;break;end
svks = @views vks[1:it]
svks = @views vks[1:it]
@timeit to "ReORTH" ReORTH(it,vkp1,svks)
tbeta = sqrt(dot(vkp1,vkp1))
tmp = 1.0/tbeta
vkp1 .*= tmp
vkp1 .*= 1.0/tbeta
if checkorth; Check_Orthogonality(it,vks,en); end
Tmat[it+1,it] = tbeta; Tmat[it,it+1] = tbeta
if doubleLanczos
if tbeta < Jtol;TF[1]=true;elit=it;break;end
Expand Down
44 changes: 27 additions & 17 deletions src/ShellModel/shellmodel_main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ end

"""
main_sm(sntf,target_nuc,num_ev,target_J;
save_wav=false,q=1,is_block=false,is_show=false,num_history=3,lm=100,ls=20,tol=1.e-6,
save_wav=false,q=1,is_block=false,is_show=false,num_history=3,lm=100,ls=20,tol=1.e-8,
in_wf="",mdimmode=false,calc_moment = false,gfactors = [1.0,0.0,5.586,-3.826],effcharge=[1.5,0.5])
Digonalize the model-space Hamiltonian
Expand All @@ -64,7 +64,7 @@ Digonalize the model-space Hamiltonian
- `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 = 15` number of vectors to be used for Thick-Restart
- `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
Expand All @@ -83,7 +83,7 @@ function main_sm(sntf,target_nuc,num_ev,target_J;save_wav=false,
Anum = parse(Int64, match(reg,target_nuc).match)
lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsmsnt(sntf,Anum)
massformula = 1
if 16<= Anum <= 40; massformula = 2;end
if 16 <= Anum <= 40; massformula = 2;end
## massformula=2: J. Blomqvist and A. Molinari, Nucl. Phys. A106, 545 (1968).
## we use this for the sd-shell nuclei
hw, bpar = init_ho_by_mass(Anum,massformula)
Expand Down Expand Up @@ -142,7 +142,7 @@ function main_sm(sntf,target_nuc,num_ev,target_J;save_wav=false,
en =[ [1.e+4 for j=1:num_ev] for i = 1:num_history]
Rvecs = [ zeros(Float64,mdim) for i=1:num_ev]
Tmat = zeros(Float64,lm,lm)
vks =[]; uks=[];itmin = 1; elit=1
vks = Vector{Float64}[]; uks=Vector{Float64}[];itmin = 1; elit=1
doubleLanczos = false
if tJ !=-1; doubleLanczos = true;end

Expand Down Expand Up @@ -261,7 +261,7 @@ To read interaction file in ".snt" format.
function readsmsnt(sntf,Anum)
f = open(sntf,"r");tlines = readlines(f);close(f)
lines = rm_comment(tlines)
line = lines[1]
line = lines[1]
lp,ln,cp,cn = map(x->parse(Int,x),rm_nan(split(line," ")))
p_sps = [zeros(Int64,4) for i=1:lp]
n_sps = [zeros(Int64,4) for i=1:ln]
Expand Down Expand Up @@ -329,7 +329,8 @@ end
"""
def_mstates(p_sps,n_sps)
to define the single particle states specified by m_z
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 j states, [n,l,j,tz].
"""
function def_mstates(p_sps,n_sps)
mstates_p = Vector{Int64}[] ; mstates_n = Vector{Int64}[]; mz_p = Int64[]; mz_n = Int64[]
Expand Down Expand Up @@ -426,7 +427,7 @@ function HbitT1(p_sps::Array{Array{Int64,1}},
loff = loffs[vrank]
vecs= [ [ [ false for i = 1:lp] for j=1:2],
[ [ false for i = 1:ln] for j=1:2]]
blist = bit2b[]
blist=bit2b[]
Vs=Float64[]
@inbounds for (i,ME) in enumerate(TBMEs[vrank])
a,b,c,d,totJ,dummy = labels[vrank][i]
Expand Down Expand Up @@ -474,7 +475,7 @@ function HbitT1(p_sps::Array{Array{Int64,1}},
bV1[vrank] = blist
V1[vrank] = Vs
end
return bV1,V1 #bVpp,Vpp,bVnn,Vnn
return bV1,V1
end

"""
Expand Down Expand Up @@ -865,6 +866,19 @@ function ReORTH(it,vtarget,vks)
return nothing
end

function Check_Orthogonality(it::Int,vks,en)
print_vec("it = $it",en[1])
svks = @views vks[1:it+1]
for i = 1:it+1
for j = i:it+1
tdot = dot(svks[i],svks[j])
if (abs(tdot) > 1.e-10 && i != j) || (i==j && abs(1-tdot) > 1.e-10)
println("dot(",@sprintf("%3i",i),",",@sprintf("%3i",j),") = ",@sprintf("%15.4e",tdot))
end
end
end
return nothing
end

function myQR!(Q,R::FA2,
d1::Int64,d2::Int64) where{FA2<:Array{Float64,2}}
Expand All @@ -889,8 +903,7 @@ function myQR!(Q,R::FA2,
return nothing
end

function bl_QR!(Q,R::FA2,
d1::Int64,d2::Int64) where{FA2<:Array{Float64,2}}
function bl_QR!(Q,R::FA2,d1::Int64,d2::Int64) where{FA2<:Array{Float64,2}}
R .= 0.0
@inbounds for j = 1:d2
q = @views Q[:,j]
Expand All @@ -909,8 +922,7 @@ function bl_QR!(Q,R::FA2,
return nothing
end

function bisearch!(v::Array{Int64,1}, target::Int64,
ret::Array{Int64,1})
function bisearch!(v::Array{Int64,1}, target::Int64,ret::Array{Int64,1})
hi = length(v); lo = 1
ret[2] = 1; ret[3]=length(v)
@inbounds while ret[2] <= ret[3]
Expand All @@ -927,8 +939,7 @@ function bisearch!(v::Array{Int64,1}, target::Int64,
return nothing
end

function bisearch_ord!(v::Array{Int64,1}, target::Int64,
ret::Array{Int64,1})
function bisearch_ord!(v::Array{Int64,1}, target::Int64, ret::Array{Int64,1})
ret[3] = length(v); ret[2] = 1
@inbounds while ret[2] <= ret[3]
ret[1] = div(ret[2]+ret[3],2)
Expand All @@ -948,7 +959,7 @@ function func_j(j1::I,j2::I,m1::I,m2::I) where{I<:Int64}
sqrt( (0.5*j1*(0.5*j1+1.0)-0.5*m1*(0.5*m1-1.0))*(0.5*j2*(0.5*j2+1.0)-0.5*m2*(0.5*m2+1.0)) )
end

function J_from_JJ1(JJ,tol=1.e-6)
function J_from_JJ1(JJ,tol=1.e-8)
for J = 0:100 ## ad hoc J<=50
hJ = 0.5*J
if abs(hJ*(hJ+1.0)-JJ) <tol
Expand Down Expand Up @@ -1088,8 +1099,7 @@ function calc_1b_jumps!(bi::Int64,bf::Int64,j::Int64,
return hits
end

function add_bl_T!(q::Int64,k::Int64,
Tmat::FA2,R::FA2) where{FA2<:Array{Float64,2}}
function add_bl_T!(q::Int64,k::Int64,Tmat::FA2,R::FA2) where{FA2<:Array{Float64,2}}
Tmat[q*k+1:q*k+q,q*k-q+1:q*k] .= R
Tmat[q*k-q+1:q*k,q*k+1:q*k+q] .= R'
return nothing
Expand Down
152 changes: 152 additions & 0 deletions src/ShellModel/trans_snt_msnt.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,152 @@
function write_msnt_sps_spe(io,p_sps,n_sps,mstates_p,mstates_n,SPEs)
println(io, "!single particle space and SPEs")
println(io, "!num n l j tz mz SPE(MeV)")
ii = 0
for pn = 1:2
sps = ifelse(pn==1,p_sps,n_sps)
msps = ifelse(pn==1,mstates_p,mstates_n)
SPE = SPEs[pn]
for i = 1:length(msps)
n,l,j,tz,mz,idx = msps[i]
ii += 1
println(io,@sprintf("%4i",ii),@sprintf("%4i",n), @sprintf("%4i",l), @sprintf("%4i",j), @sprintf("%4i",tz), @sprintf("%5i",mz),@sprintf("%15.6f",SPE[idx]))
end
end
return nothing
end

function trans_snt_msnt(fn,Anum,Mtot,p_sps,n_sps,mstates_p,mstates_n,SPEs,olabels,oTBMEs)
io = open(replace(fn,".snt"=>".msnt"),"w")
println(io,"!snt:$fn")
println(io,"!totalM = $Mtot")
write_msnt_sps_spe(io,p_sps,n_sps,mstates_p, mstates_n,SPEs)
write_msnt_tbmes(io,Mtot,p_sps,n_sps,mstates_p,mstates_n,SPEs,olabels,oTBMEs)
close(io)
return nothing
end

struct ket_abJ
a::Int64
b::Int64
J::Int64
end

function write_msnt_tbmes(io,Mtot,p_sps,n_sps,mstates_p,mstates_n,SPEs,olabels,oTBMEs)
lp = length(p_sps)
lpm = length(mstates_p)
m_sps = vcat(mstates_p,mstates_n)
allSPEs = vcat(SPEs[1],SPEs[2])
dict_msps2sps = Dict{Int64,Int64}()
for (idx,tmp) in enumerate(mstates_p)
dict_msps2sps[idx] = tmp[end]
end
for (idx,tmp) in enumerate(mstates_n)
dict_msps2sps[idx+lpm] = tmp[end] + lp
end

dictTBMEs = Dict{Vector{Int64},Float64}()
for (i,tkey) in enumerate(olabels)
a,b,c,d,J,oidx = tkey
dictTBMEs[[a,b,c,d,J]] = oTBMEs[i]
dictTBMEs[[c,d,a,b,J]] = oTBMEs[i]
end

kets = [ ket_abJ[ ] for pnrank=1:3]
for pnrank = 1:3
ctxt = "pp:"
if pnrank==2; ctxt = "pn:";end
if pnrank==3; ctxt = "nn:";end
msps_a = ifelse(pnrank<=2,mstates_p,mstates_n)
msps_b = ifelse(pnrank>=2,mstates_n,mstates_p)
ofst_a = ifelse(pnrank<=2,0,lpm)
ofst_b = ifelse(pnrank>=2,lpm,0)
for (ia,a) in enumerate(msps_a)
na,la,ja,tza,mza,idx_a = a
ia += ofst_a
for (ib,b) in enumerate(msps_b)
ib += ofst_b
if ia > ib; continue;end
nb,lb,jb,tzb,mzb,idx_b = b
if Mtot != mza + mzb; continue;end
Tz = tza + tzb
for J = abs(div(ja-jb,2)):div(ja+jb,2)
if pnrank % 2 == 1 && J % 2 == 1; continue; end
push!(kets[pnrank],ket_abJ(ia,ib,J))
end
end
end
dim = length(kets[pnrank])
mat = zeros(Float64,dim,dim)
idxs = [ Int64[ ] for J =0:10]
for (idx_bra,bra) in enumerate(kets[pnrank])
a = bra.a
b = bra.b
sps_a = dict_msps2sps[a]
sps_b = dict_msps2sps[b]
ja,ma = m_sps[a][3:2:5]
jb,mb = m_sps[b][3:2:5]
Jbra = bra.J
push!(idxs[Jbra+1],idx_bra)
for (idx_ket,ket) in enumerate(kets[pnrank])
if idx_bra > idx_ket; continue;end
c = ket.a
d = ket.b
sps_c = dict_msps2sps[c]
sps_d = dict_msps2sps[d]
jc,mc = m_sps[c][3:2:5]
jd,md = m_sps[d][3:2:5]
Jket = ket.J
if Jbra != Jket; continue;end
tbme = dictTBMEs[ [sps_a,sps_b,sps_c,sps_d,Jket] ]
CG1 = clebschgordan(Float64,ja//2,ma//2,jb//2,mb//2, Jket, Mtot)
CG2 = clebschgordan(Float64,jc//2,mc//2,jd//2,md//2, Jket, Mtot)
tbme_M = tbme * sqrt( (1.0+deltaf(sps_a,sps_b)) *(1.0+deltaf(sps_c,sps_d)) ) * CG1 * CG2
println(io, ctxt, @sprintf("%5i",a),@sprintf("%5i",b),@sprintf("%5i",c),@sprintf("%5i",d), @sprintf("%5i",Jket), @sprintf("%15.6f",tbme_M))
mat[idx_bra,idx_ket] = mat[idx_ket,idx_bra] = tbme_M
end
end
if pnrank == 2
for J = 0:10
idx = idxs[J+1]
if length(idx) == 0; continue; end
submat = mat[idx,idx]
evals,evecs = eigen(submat)
println("J $J ")
for nn = 1:length(evals)
if abs(evals[nn]) < 1.e-12; continue;end
println(@sprintf("%15.6f", evals[nn]) )
# vec = evecs[:,nn]
# E1b = 0.0
# for tmp in idx
# tket = kets[pnrank][tmp]
# a = dict_msps2sps[tket.a]
# b = dict_msps2sps[tket.b]
# E1b += (allSPEs[a] + allSPEs[b])
# end
# E = E1b + evals[nn]
# println("E1b ", @sprintf("%15.6f",E1b), " E ",@sprintf("%15.6f",E))
end
end
end
end
return nothing
end

function main_trans_msnt(fn,target_nuc,target_Js=[])
Anum = parse(Int64, match(reg,target_nuc).match)
Mtot = 0;if Anum % 2 != 0; Mtot = 1;end
if Anum % 2 != Mtot % 2
@error "invalid target_Js = $target_Js"
end
if length(target_Js) > 0
Mtot = minimum(target_Js)
tJ=target_Js[1]
eval_jj = 0.5*tJ*(tJ/2+1)
end
lp,ln,cp,cn,massop,Aref,p,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsmsnt(fn,Anum)

target_el = replace.(target_nuc, string(Anum)=>"")
Z,N,vp,vn = getZNA(target_el,Anum,cp,cn)
mstates_p,mstates_n,mz_p,mz_n = def_mstates(p_sps,n_sps)
trans_snt_msnt(fn,Anum,Mtot,p_sps,n_sps,mstates_p,mstates_n,SPEs,olabels,oTBMEs)
end
2 changes: 1 addition & 1 deletion src/hartreefock.jl/io_input.jl
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ function rm_comment(lines)
nlines = []
for line in lines
line = strip(line)
if length(line)>1
if length(line) > 0
if startswith(line,"!")
continue
end
Expand Down

0 comments on commit eed67f8

Please sign in to comment.