diff --git a/docs/src/contributing.md b/docs/src/contributing.md index 60bd78e7..7f837ad5 100644 --- a/docs/src/contributing.md +++ b/docs/src/contributing.md @@ -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 diff --git a/example/log_sample_script.txt b/example/log_sample_script.txt index 9e2f0544..8f37535c 100644 --- a/example/log_sample_script.txt +++ b/example/log_sample_script.txt @@ -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 diff --git a/src/NuclearToolkit.jl b/src/NuclearToolkit.jl index 8beea5cd..0f33cbd8 100644 --- a/src/NuclearToolkit.jl +++ b/src/NuclearToolkit.jl @@ -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 diff --git a/src/ShellModel/betadecay.jl b/src/ShellModel/betadecay.jl index 99a14660..b26ccc28 100644 --- a/src/ShellModel/betadecay.jl +++ b/src/ShellModel/betadecay.jl @@ -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) @@ -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) diff --git a/src/ShellModel/lanczos_methods.jl b/src/ShellModel/lanczos_methods.jl index 68529c9f..a3bee17e 100644 --- a/src/ShellModel/lanczos_methods.jl +++ b/src/ShellModel/lanczos_methods.jl @@ -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] @@ -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 diff --git a/src/ShellModel/shellmodel_main.jl b/src/ShellModel/shellmodel_main.jl index bb113221..d9ed260e 100644 --- a/src/ShellModel/shellmodel_main.jl +++ b/src/ShellModel/shellmodel_main.jl @@ -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 @@ -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 @@ -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) @@ -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 @@ -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] @@ -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[] @@ -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] @@ -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 """ @@ -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}} @@ -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] @@ -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] @@ -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) @@ -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) ".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 \ No newline at end of file diff --git a/src/hartreefock.jl/io_input.jl b/src/hartreefock.jl/io_input.jl index e7d410d8..fa348fa9 100644 --- a/src/hartreefock.jl/io_input.jl +++ b/src/hartreefock.jl/io_input.jl @@ -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