From 422df1b7f3973f65478c379afb59f6a79eee4235 Mon Sep 17 00:00:00 2001 From: SotaYoshida Date: Tue, 5 Mar 2024 12:55:29 +0900 Subject: [PATCH 1/4] redundant stdout removed --- src/ShellModel/eigenvector_continuation.jl | 18 +- src/ShellModel/eigenvector_continuation.jl~ | 1357 +++++++++++++++++ src/ShellModel/lanczos_methods.jl | 1 - src/ShellModel/lanczos_methods.jl~ | 871 +++++++++++ test/HFMBPT_IMSRG_test.jl | 4 +- test/HFMBPT_IMSRG_test.jl~ | 76 + test/ShellModel_test.jl | 2 +- test/parameters/optional_parameters_forDMD.jl | 4 +- .../parameters/optional_parameters_forDMD.jl~ | 16 + test/runtests.jl | 3 - 10 files changed, 2337 insertions(+), 15 deletions(-) create mode 100644 src/ShellModel/eigenvector_continuation.jl~ create mode 100644 src/ShellModel/lanczos_methods.jl~ create mode 100644 test/HFMBPT_IMSRG_test.jl~ create mode 100644 test/parameters/optional_parameters_forDMD.jl~ diff --git a/src/ShellModel/eigenvector_continuation.jl b/src/ShellModel/eigenvector_continuation.jl index 61469ffc..df918821 100644 --- a/src/ShellModel/eigenvector_continuation.jl +++ b/src/ShellModel/eigenvector_continuation.jl @@ -100,11 +100,11 @@ function readRvecs!(mdim,inpf,tJ,vecs,wfinfos;num=100) end function prepEC(Hs,target_nuc,num_ev,num_ECsample,tJ,mode; - num_history=3,lm=100,ls=15,tol=1.e-9,q=1, + num_history=3,lm=400,ls=15,tol=1.e-9,q=1, path_to_samplewav="",calc_moment=true, save_wav = false,tdmatdir="./tdmat/", gfactors = [1.0,0.0,5.586,-3.826], - effcharge=[1.5,0.5],debugmode=false,doublelanczos=true) + effcharge=[1.5,0.5],debugmode=false,doublelanczos=true,is_show=false) to = TimerOutput() sntf = Hs[1] Anum = parse(Int64, match(reg,target_nuc).match) @@ -292,7 +292,9 @@ function prepEC(Hs,target_nuc,num_ev,num_ECsample,tJ,mode; close(io) end end - show_TimerOutput_results(to) + if is_show + show_TimerOutput_results(to) + end return nothing end @@ -448,7 +450,9 @@ function solveEC(Hs,target_nuc,tJNs; if length(Hs) > 1 && exact_logf != "" @timeit to "scatter plot" plot_EC_scatter(target_nuc,Hs,sumV,tJNs,Dims,exlines) end - show_TimerOutput_results(to) + if is_show + show_TimerOutput_results(to) + end return nothing end @@ -583,8 +587,10 @@ function solveEC_UQ(Hs,target_nuc,tJNs,Erefs,errors; # "std: ", @sprintf("%10.4f ", std([iThetas[i][n] for i =1:Ns]))) #end println(tx) - write_history(iThetas,Ts,llhs,tJNs,Ens) - show_TimerOutput_results(to) + write_history(iThetas,Ts,llhs,tJNs,Ens) + if is_show + show_TimerOutput_results(to) + end println("") return nothing end diff --git a/src/ShellModel/eigenvector_continuation.jl~ b/src/ShellModel/eigenvector_continuation.jl~ new file mode 100644 index 00000000..8b10ddf2 --- /dev/null +++ b/src/ShellModel/eigenvector_continuation.jl~ @@ -0,0 +1,1357 @@ +function J2_to_J(J2) + J="" + if Int(J2) % 2 == 0 + J=string(Int(div(J2,2))) + else + J=string(Int(J2))*"/2" + end + return J +end + +function myCholesky!(tmpA,ln::Int64,cLL::LowerTriangular{Float64,Array{Float64,2}}) + l11 = sqrt(tmpA[1,1]) + cLL[1,1] = l11 + cLL[2,1] = tmpA[2,1]/l11; cLL[2,2] = sqrt( tmpA[2,2]-cLL[2,1]^2) + for i=3:ln + for j=1:i-1 + cLL[i,j] = tmpA[i,j] + for k = 1:j-1 + cLL[i,j] += - cLL[i,k]*cLL[j,k] + end + cLL[i,j] = cLL[i,j] / cLL[j,j] + end + cLL[i,i] = tmpA[i,i] + for j=1:i-1 + cLL[i,i] += -cLL[i,j]^2 + end + cLL[i,i] = sqrt(cLL[i,i]) + end + nothing +end + +function read_wfinfos(inpf) + f = open(inpf,"r");tlines = readlines(f);close(f) + lines = rm_comment(tlines) + ret = [] + for line in lines + tl = rm_nan(split(line," ")) + push!(ret,[tl[1],parse(Int64,tl[3]),parse(Int64,tl[5])]) + end + return ret +end + +function read_tdmat(inpf) + io=open(inpf,"r") + Dim=read(io,Int) + ln = read(io,Int) + lTD = read(io,Int) + TDs = [ Float64[] for i=1:ln] + infos = [ Int64[] for i=1:ln] + Nmat = zeros(Float64,Dim,Dim) + for nth = 1:ln + i = read(io,Int); j = read(io,Int); Nij = read(io,Int) + infos[nth] = [i,j,Nij] + tN = read(io,Float64) + Nmat[i,j] = tN + Nmat[j,i] = tN + TDs[nth] = [ read(io,Float64) for k=1:lTD] + end + close(io) + return Dim,ln,lTD,Nmat,infos,TDs +end + +function read_wf_from_info(wfinfos,mdim,vs,wpath) + ## e.g., ./wavsamples_sigma3/Mg24_sample_1.wav nth= 1 TDmatidx 1 + for tmp in wfinfos + inpf,nth,TDidx = tmp + if wpath != "./" && wpath != "" + tn = split(inpf,"/") + inpf = wpath* ifelse(wpath[end]=="/","","/") * tn[end] + end + io = open(inpf,"r") + num_ev = read(io,Int32) + mtot = read(io,Int32) + vals = [read(io,Float64) for i = 1:num_ev] + ojs = [read(io,Int32) for i=1:num_ev] + Css = [ [read(io,Float64) for i=1:mdim] for j=1:num_ev] + vs[TDidx] .= Css[nth] + close(io) + end + return nothing +end + +function readRvecs!(mdim,inpf,tJ,vecs,wfinfos;num=100) + io = open(inpf,"r") + num_ev = read(io,Int32) + mtot = read(io,Int32) + vals = [read(io,Float64) for i = 1:num_ev] + ojs = [read(io,Int32) for i=1:num_ev] + Css = [ [read(io,Float64) for i=1:mdim] for j=1:num_ev] + n = minimum([num,num_ev]) + for i=1:n + j = Int64(ojs[i]) + if j == tJ + push!(vecs,Css[i]) + push!(wfinfos,[inpf,i,length(vecs)]) + end + end + close(io) + return nothing +end + +function prepEC(Hs,target_nuc,num_ev,num_ECsample,tJ,mode; + num_history=3,lm=100,ls=15,tol=1.e-9,q=1, + path_to_samplewav="",calc_moment=true, + save_wav = false,tdmatdir="./tdmat/", + gfactors = [1.0,0.0,5.586,-3.826], + effcharge=[1.5,0.5],debugmode=false,doublelanczos=true) + to = TimerOutput() + sntf = Hs[1] + 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) + hw, bpar = init_ho_by_mass(Anum,1) # mass formula + if 16<= Anum <= 40 + hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell + end + Mtot = 0;if Anum % 2 != 0; Mtot = 1;end + if typeof(tJ) != Int64;println("targetJ must be integer!!");exit();end + Mtot = tJ + eval_jj = 0.5*tJ*(tJ/2+1) + if Anum % 2 != Mtot % 2; println("invalid target J");exit();end + + target_el = replace.(target_nuc, string(Anum)=>"") + Z,N,vp,vn = getZNA(target_el,Anum,cp,cn) + msps_p,msps_n,mz_p,mz_n = construct_msps(p_sps,n_sps) + pbits,nbits,jocc_p,jocc_n,Mps,Mns,tdims = occ(p_sps,msps_p,mz_p,vp, + n_sps,msps_n,mz_n,vn,Mtot) + oPP,oNN,oPNu,oPNd = prep_J(tdims,p_sps,n_sps,msps_p,msps_n, + pbits,nbits) + mdim = tdims[end] + lblock=length(pbits) + mdim_print(target_nuc,Z,N,cp,cn,vp,vn,mdim,tJ) + if mode == "sample" + for (iter,sntf) in enumerate(Hs) + @timeit to "make sample" begin + println("sntf: $sntf") + lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsmsnt(sntf,Anum) + pp_2bjump,nn_2bjump = prep_Hamil_T1(p_sps,n_sps,msps_p,msps_n,pbits,nbits,labels,TBMEs) + bVpn,Vpn,delMs = prep_Hamil_pn(p_sps,n_sps,msps_p,msps_n,labels[3],TBMEs[3]) + l_pbit = length(msps_p);l_nbit = length(msps_n) + bis,bfs,p_NiNfs,n_NiNfs,num_task = prep_pn(lblock,pbits,nbits,Mps,delMs,bVpn,Vpn) + bVpn=nothing + block_tasks = make_distribute(num_task) + + Js = [ 0.5*Mtot*(0.5*Mtot+1) for i = 1:num_ev] + Jtasks = zeros(Int64,lblock) + for i = 1:lblock + Jtasks[i] = length(pbits[i])*length(nbits[i]) + end + Jidxs = make_distribute_J(Jtasks) + + 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) + itmin = 1; elit=1 + vks = [ zeros(Float64,mdim) for i=1:lm] + uks = [ zeros(Float64,mdim) for i=1:ls] + initialize_wf(vks[1],"rand",tJ,mdim) + + @timeit to "Lanczos" elit = TRL(vks,uks,Tmat,itmin, + pbits,nbits,jocc_p,jocc_n,SPEs, + pp_2bjump,nn_2bjump,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=doublelanczos,debugmode=debugmode) + Rvecs = @views uks[1:num_ev] + vals,vecs = eigen(@views Tmat[1:elit*q,1:elit*q]) + @inbounds for (nth,Rvec) in enumerate(Rvecs) + Rvec .= 0.0 + @inbounds for k=1:length(vals) + Rvec .+= vecs[k,nth] .* vks[k] + end + Rvec .*= 1.0/sqrt(dot(Rvec,Rvec)) + end + Jv = zeros(Float64,mdim) + for (nth,Rv) in enumerate(Rvecs) + Jv .= 0.0 + operate_J!(Rv,Jv,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd) + Js[nth] += dot(Rv,Jv) + end + totJs = J_from_JJ1.(Js) + println("J $totJs")#;println("Js $Js") + print("En. ");map(x -> @printf("%24.15f ",x), en[1]);print("\n") + csnt = split(split(sntf,"/")[end],".")[1] + oupf= path_to_samplewav*target_nuc*"_"*csnt*".wav" + if tJ != -1;oupf=path_to_samplewav*target_nuc*"_"*csnt*"_j"*string(tJ)*".wav";end + if save_wav + writeRitzvecs(mdim,Mtot,en[1],totJs,Rvecs,oupf) + end + if calc_moment + @timeit to "moment" tx_mom = eval_moment(Mtot,Rvecs,totJs,p_sps,n_sps, + msps_p,msps_n,tdims, + jocc_p,jocc_n,pbits,nbits,bpar, + gfactors,effcharge) + if tx_mom != "";print(tx_mom);end + end + println("") + end + end + elseif mode =="TD" + try + mkdir(tdmatdir) + catch + nothing + end + #@timeit to "read sample wavs." begin + svecs = Vector{Float64}[ ] + wfinfos = [ ] + for fidx = 0:num_ECsample-1 + inpf = path_to_samplewav*target_nuc*"_tmp_$fidx"*"_j"*string(tJ)*".wav" + readRvecs!(mdim,inpf,tJ,svecs,wfinfos;num=num_ev) + end + + io=open(tdmatdir*"wfinfos_"*target_nuc*"_J"*string(tJ)*".dat","w") + for tmp in wfinfos + println(io,split(tmp[1],"/")[end], " nth= ",tmp[2], " TDmatidx ",tmp[3]) + end + close(io) + @timeit to "misc" begin + Dim = length(svecs) + Hmat = zeros(Float64,Dim,Dim) + Nmat = zeros(Float64,Dim,Dim) + println("Dim $Dim") + end + @timeit to "Nmat calc." begin + @threads for j = 1:Dim + v = svecs[j] + for i = 1:j + tN = dot(svecs[i],v) + Nmat[i,j] = tN; Nmat[j,i]=tN[1] + end + end + end + @timeit to "OBTD" begin + numv = length(svecs) + idxs = Tridx[] + for i = 1:length(svecs) + for j=i:length(svecs) + push!(idxs, Tridx(i,j,idx_from_ij(i,j,numv))) + end + end + OBTDs = [ [ zeros(Float64,length(p_sps)), + zeros(Float64,length(n_sps)) + ] for Nij = 1:length(idxs) ] + calcOBTD(OBTDs,idxs,p_sps,n_sps, + msps_p,msps_n,tdims, + jocc_p,jocc_n,pbits,nbits,svecs) + end + + @timeit to "prepTBTD" begin + opTBTD1,pjumps,njumps,bif_idxs = prepTBTD(tJ,idxs,p_sps,n_sps, + msps_p,msps_n, + pbits,nbits, + labels,TBMEs,svecs, + tdims,Mps) + TBTDs= [zeros(Float64,length(oTBMEs)) for Nij = 1:length(idxs)] + end + @timeit to "calcTBTD" begin + calcTBTD(TBTDs,opTBTD1,pjumps,njumps, + pbits,nbits,tdims,svecs, + idxs,bif_idxs, + olabels,oTBMEs,labels,to) + end + @timeit to "write TDs" begin + oupf=tdmatdir*target_nuc*"_J"*string(tJ)*".dat" + lTD = length(p_sps)+length(n_sps)+length(oTBMEs) + io=open(oupf,"w") + write(io,Dim) + write(io,length(idxs)) + write(io,lTD) + for i = 1:length(svecs) + for j=i:length(svecs) + Nij = idx_from_ij(i,j,length(svecs)) + write(io,i);write(io,j) + write(io,Nij) + write(io,Nmat[i,j]) + TDs = Float64[] + for nth=1:length(OBTDs[Nij][1]) + push!(TDs,OBTDs[Nij][1][nth]) + end + for nth=1:length(OBTDs[Nij][2]) + push!(TDs,OBTDs[Nij][2][nth]) + end + for nth=1:length(TBTDs[Nij]) + push!(TDs,TBTDs[Nij][nth]) + end + for TD in TDs + write(io,TD) + end + end + end + close(io) + end + end + if is_show + show_TimerOutput_results(to) + end + return nothing +end + +""" + 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. +```math +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 + +!!! note + All the effective interactions must be in "the same order" and must be consistent with interaction file from which the transition density matrices were made. + +""" +function solveEC(Hs,target_nuc,tJNs; + write_appwav=false,verbose=false, + calc_moment=true,wpath="./",tdmatpath="./", + is_show=false, + gfactors = [1.0,0.0,5.586,-3.826], + effcharge=[1.5,0.5], + exact_logf="") + + if length(tJNs)==0;println("specify targetJ and # of states");exit();end + to = TimerOutput() + sntf = Hs[1] + 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) + hw, bpar = init_ho_by_mass(Anum,1) # mass formula + if 16<= Anum <= 40 + hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell + end + target_el = replace.(target_nuc, string(Anum)=>"") + Z,N,vp,vn = getZNA(target_el,Anum,cp,cn) + msps_p,msps_n,mz_p,mz_n = construct_msps(p_sps,n_sps) + #println("nuc: $target_el") + exlines = "" + try + f = open(exact_logf,"r");exlines = readlines(f);close(f) + catch + nothing + #println("logfile of exact calc. is missing ") + end + + Dims = Int64[] + sumV=[] + wfinfos=[] + lTD = 0 + + for (jidx,tmp) in enumerate(tJNs) + tJ,num_ev_target = tmp + vals = zeros(Float64,num_ev_target) + Mtot = tJ + pbits,nbits,jocc_p,jocc_n,Mps,Mns,tdims = occ(p_sps,msps_p,mz_p,vp,n_sps,msps_n,mz_n,vn,Mtot) + mdim = tdims[end] + @timeit to "read TDmat" begin + inpf = tdmatpath*target_nuc*"_J"*string(tJ)*".dat" + Dim,ln,lTD,Nmat,infos,TDs = read_tdmat(inpf) + push!(Dims,Dim) + if Dim+1 <= num_ev_target; + println("Error!: sample number must be larger than $num_ev_target") + exit() + end + end + #println("for J=$tJ/2 states Dim.=$Dim") + @timeit to "Cholesky" begin + tMat = zeros(Float64,Dim,Dim) + tildH = zeros(Float64,Dim,Dim) + L = LowerTriangular(zeros(Float64,Dim,Dim)) + try + myCholesky!(Nmat,Dim,L) + catch + u = 1.1 * 1.e-16 + shift = tr(Nmat) * u * ( (Dim+2)/(1.0-((Dim+1)*(Dim+2))*u)) + #println("#epsilon prescription is used!! shift: $shift" ) + Nmat = Nmat + shift .* Matrix(1.0I, Dim,Dim) + myCholesky!(Nmat,Dim,L) + end + Linv = inv(L) + end + Hmat = zeros(Float64,Dim,Dim) + + for sntidx = 1:length(Hs) + sntf = Hs[sntidx] + tMat .= 0.0 + @timeit to "Hmat" begin + lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsmsnt(sntf,Anum) + MEs = [SPEs[1],SPEs[2],oTBMEs] + MEs = [(MEs...)...] + for (idx,TD) in enumerate(TDs) + i,j,Nij = infos[idx] + tmp = dot(TD,MEs) + Hmat[i,j] = tmp; Hmat[j,i] = tmp + end + @timeit to "Gen. Eigen" begin + mul!(tMat,Linv,Hmat) + mul!(tildH,Linv',tMat) + vals,vecs = Arpack.eigs(tildH,nev=num_ev_target,which=:SR,tol=1.e-8, maxiter=300) + vals = real.(vals) + println("vals $vals ln $(length(vals)) num_ev_target $num_ev_target") + if verbose + print_vec("$target_nuc 2J=$tJ En(EC)",vals) + end + end + push!(sumV,[sntf,tJ,vals]) + end + ## write out approx. wav. for only the first Hamiltonian + if write_appwav && sntidx == 1 + wfinfos = read_wfinfos("./tdmat/wfinfos_"*target_nuc*"_J"*string(tJ)*".dat") + @timeit to "read sample wavs." begin + vs = [zeros(Float64,mdim) for i=1:Dim] + read_wf_from_info(wfinfos,mdim,vs,wpath) + end + @timeit to "make .appwav" begin + Rvecs = [ zeros(Float64,mdim) for i=1:length(vals)] + for nth = 1:length(vals) + for k = 1:length(vs) + @. Rvecs[nth] += vecs[k,nth] * vs[k] + end + Rvecs[nth] ./= sqrt(sum(Rvecs[nth].^2)) + end + end + @timeit to "write .appwav" begin + csnt = split(split(sntf, "/")[end],".")[1] + wavname = "./appwavs/"*target_nuc*"_"*csnt*"_j"*string(tJ)*".appwav" + Jlist = [tJ for i=1:length(vals)] + writeRitzvecs(mdim,Mtot,vals,Jlist,Rvecs,wavname) + end + end + end + end + if length(Hs) > 1 && exact_logf != "" + @timeit to "scatter plot" plot_EC_scatter(target_nuc,Hs,sumV,tJNs,Dims,exlines) + end + if is_show + show_TimerOutput_results(to) + end + return nothing +end + + +function solveEC_UQ(Hs,target_nuc,tJNs,Erefs,errors; + itnum_MCMC = 5000,burnin=500,var_M=2.e-2, + calc_moment=true,wpath="./", + is_show=false, + gfactors = [1.0,0.0,5.586,-3.826], + effcharge=[1.5,0.5], + exact_logf="", + intf="", + num_replica=10) + UQ=true; write_appwav=false + if length(tJNs)==0;println("specify targetJ and # of states");exit();end + to = TimerOutput() + + sntf = Hs[1] + Anum = parse(Int64, match(reg,target_nuc).match) + lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsnt(sntf,Anum) + hw, bpar = init_ho_by_mass(Anum,1) # mass formula + if 16<= Anum <= 40 + hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell + end + target_el = replace.(target_nuc, string(Anum)=>"") + Z,N,vp,vn = getZNA(target_el,Anum,cp,cn) + msps_p,msps_n,mz_p,mz_n = construct_msps(p_sps,n_sps) + println("nuc: $target_nuc") + + Dims = Int64[] + lTD = 0 + lnJ = length(tJNs) + HNmats = [ [ zeros(Float64,2,2) ] for i =1:lnJ] + for i=1:lnJ;deleteat!(HNmats[i],1);end + AllTDs = [ [ zeros(Float64,2) ] ]; deleteat!(AllTDs,1) + Allinfos = [ [ [0,0] ] ]; deleteat!(Allinfos,1) + + for (jidx,tmp) in enumerate(tJNs) + tJ,num_ev_target = tmp + Mtot = tJ + @timeit to "read TDmat" begin + inpf = "./tdmat/"*target_nuc*"_J"*string(tJ)*".dat" + Dim,ln,lTD,Nmat,infos,TDs = read_tdmat(inpf) + push!(Dims,Dim) + if Dim+1 <= num_ev_target; + println("Error!: sample number must be larger than $num_ev_target") + exit() + end + end + @timeit to "Cholesky" begin + tMat = zeros(Float64,Dim,Dim) + tildH = zeros(Float64,Dim,Dim) + L = LowerTriangular(zeros(Float64,Dim,Dim)) + try + myCholesky!(Nmat,Dim,L) + catch + u = 1.1 * 1.e-16 + shift = tr(Nmat) * u * ( (Dim+2)/(1.0-((Dim+1)*(Dim+2))*u)) + println("#epsilon prescription is used!! shift: $shift" ) + Nmat = Nmat + shift .* Matrix(1.0I, Dim,Dim) + myCholesky!(Nmat,Dim,L) + end + Linv = inv(L) + end + Hmat = zeros(Float64,Dim,Dim) + push!(AllTDs,copy(TDs)) + push!(Allinfos,copy(infos)) + push!(HNmats[jidx],Hmat) + push!(HNmats[jidx],Linv) + end + + ## MCMC + in_snt =false + Theta = zeros(Float64,length(SPEs[1])+length(SPEs[2])+length(oTBMEs)) + lME = length(Theta) + tDim = Dims[1] + tMat = zeros(Float64,tDim,tDim) + tildH = zeros(Float64,tDim,tDim) + label_T1,label_T0 = make_int(p_sps,n_sps) + idx_s_from_i,facs = int2snt(p_sps,n_sps,label_T1,label_T0,olabels) + + ### random walk in .snt space (not used now) + #MEs = [SPEs[1],SPEs[2],oTBMEs] + #Theta .= [(MEs...)...] + #Thetas = [ zeros(Float64,lME) ]; deleteat!(Thetas,1) + + ### random walk in .int space + ### normal Metropolis-Hastings + #c_Theta = zeros(Float64,lME); c_Theta .= Theta + # iSPEs,V1s,V0s = readVint(intf,label_T1,label_T0) + # ln_int = length(iSPEs)+length(V1s)+length(V0s) + # Vint = zeros(Float64,ln_int) + # iThetas = [ zeros(Float64,ln_int)];deleteat!(iThetas,1) + # nSPEs = copy(iSPEs); nV1s = copy(V1s); nV0s = copy(V0s) + # evalVsnt!(Theta,iSPEs,V1s,V0s,idx_s_from_i,facs) + # tx, llhs, Ens = intMCMC(itnum_MCMC,burnin,var_M,iThetas,Theta,c_Theta,Vint, + # nSPEs,nV1s,nV0s,iSPEs,V1s,V0s, + # idx_s_from_i,facs, + # lnJ,tJNs,HNmats,AllTDs,Allinfos, + # tMat,tildH,Erefs,errors,to) + # plot_MCMC(iThetas,llhs,tJNs,Ens,Erefs) + + ### parallel tempering + oSPEs,oV1s,oV0s = readVint(intf,label_T1,label_T0) + iSPEs = [ copy(oSPEs) for i = 1:num_replica] + V1s = [ copy(oV1s) for i = 1:num_replica] + V0s = [ copy(oV0s) for i = 1:num_replica] + Ts = [ 1.0 + 0.25*(i-1) for i=1:num_replica] + #Ts = [ 1.0 *i for i=1:num_replica] + + ln_int = length(iSPEs[1])+length(V1s[1])+length(V0s[1]) + Vint = zeros(Float64,ln_int) + iThetas = Vector{Float64}[ ] + #AllThetas = [ [ zeros(Float64,ln_int)] for i = 1:num_replica] + #for i=1:num_replica; deleteat!(AllThetas[i],1); end + nSPEs = deepcopy(iSPEs); nV1s = deepcopy(V1s); nV0s = deepcopy(V0s) + evalVsnt!(Theta,iSPEs[1],V1s[1],V0s[1],idx_s_from_i,facs) + Thetas = [ copy(Theta) for i=1:num_replica] + c_Thetas = deepcopy(Thetas) + tx, llhs, Ens = intMCMC_PT(itnum_MCMC,burnin,var_M,Ts, + iThetas,Thetas,c_Thetas,Vint, + nSPEs,nV1s,nV0s,iSPEs,V1s,V0s, + oSPEs,oV1s,oV0s, + idx_s_from_i,facs, + lnJ,tJNs,HNmats,AllTDs,Allinfos, + tMat,tildH,Erefs,errors,to) + + Ns = length(iThetas) + #for n =1:ln_int + # println("i= ", @sprintf("%4i", n), + # "mean: ",@sprintf("%10.4f ", mean([iThetas[i][n] for i =1:Ns])), + # "std: ", @sprintf("%10.4f ", std([iThetas[i][n] for i =1:Ns]))) + #end + println(tx) + write_history(iThetas,Ts,llhs,tJNs,Ens) + if is_show + show_TimerOutput_results(to) + end + println("") + return nothing +end + +function write_history(iThetas,Ts,llhs,tJNs,Ens) + try + mkdir("history") + catch + nothing + end + nr = length(Ts) + Ns = length(llhs[1]) + num_states = sum([tmp[2] for tmp in tJNs]) + + oupf="./history/PTlog_Theta_lowestT.dat" + io = open(oupf,"w") + write(io,nr) + write(io,Ns) + write(io,length(iThetas[1])) + for i=1:Ns + write(io,iThetas[i]) + end + close(io) + oupf="./history/PTlog_llhs.dat" + io = open(oupf,"w") + write(io,nr) + write(io,Ns) + write(io,Ts) + for i=1:nr + write(io,llhs[i]) + end + close(io) + + oupf="./history/PTlog_Ens.dat" + io = open(oupf,"w") + write(io,nr) + write(io,Ns) + for i=1:num_states + write(io,Ens[i]) + end + close(io) +end + +function intMCMC(itnum_MCMC,burnin,var_M,iThetas,Theta,c_Theta,Vint,Vopt, + nSPEs,nV1s,nV0s,SPEs,V1s,V0s, + idx_s_from_i,facs, + lnJ,tJNs,HNmats,AllTDs,Allinfos, + tMat,tildH,Erefs,errors,to) + llh = -1.e+10 + nllh = -1.e+10 + Accept = true + Acchit = 0 + llhs = Float64[ ] + tnum_states = sum([tJNs[i][2] for i =1:lnJ]) + Ens = [ Float64[] for i=1:tnum_states] + offs = [0] + for i = 2:lnJ + push!(offs,sum([tJNs[j][2] for j =1:i-1])) + end + for itM = 1:itnum_MCMC + #println("it=$itM") + @timeit to "int." begin + proposal_Vint!(nSPEs,nV1s,nV0s,SPEs,V1s,V0s,var_M) + evalVsnt!(c_Theta,nSPEs,nV1s,nV0s,idx_s_from_i,facs) + end + nllh = 0.0 + @timeit to "EC" @inbounds for jidx = 1:lnJ + @timeit to "Hmats" begin + num_ev_target = tJNs[jidx][2] + Hmat = HNmats[jidx][1] + Linv = HNmats[jidx][2] + TDs = AllTDs[jidx] + infos = Allinfos[jidx] + @inbounds @threads for idx = 1:length(TDs) + TD = TDs[idx] + i,j,Nij = infos[idx] + tmp = dot(TD,c_Theta) + Hmat[i,j] = tmp; Hmat[j,i] = tmp + end + end + BLAS.gemm!('N','N',1.0,Linv,Hmat,0.0,tMat) + BLAS.gemm!('T','N',1.0,Linv,tMat,0.0,tildH) + + #@timeit to "Arpack" vals = real.(Arpack.eigs(tildH,nev=num_ev_target,which=:SR,tol=1.e-8, maxiter=300))[1] + vals = real.(eigsolve(tildH,num_ev_target,:SR;tol=1.e-8, maxiter=300)[1])[1:num_ev_target] + + nllh += L2_llh(vals,Erefs[jidx],errors[jidx]) + if itM > burnin + @timeit to "push" for i=1:num_ev_target + push!(Ens[offs[jidx]+i],vals[i]) + end + end + end + if nllh > llh + Accept = true + elseif nllh-llh > log(rand()) + Accept = true + else + Accept = false + end + if Accept + llh = nllh + Theta .= c_Theta + SPEs .= nSPEs + V1s .= nV1s + V0s .= nV0s + Acchit += 1 + evalVint!(Vint,SPEs,V1s,V0s) + end + @timeit to "int." if itM > burnin + push!(iThetas,copy(Vint)) + push!(llhs,llh) + end + end + tx = "Accept rate" * @sprintf("%5.1f ",100*(Acchit-burnin)/(itnum_MCMC-burnin)) + return tx, llhs, Ens +end + +function singleH(TD,i,j,Nij,c_Theta, Hmat) + tmp = dot(TD,c_Theta) + Hmat[i,j]=tmp;Hmat[j,i]=tmp + return nothing +end + +function intMCMC_PT(itnum_MCMC,burnin,var_M,Ts, + iThetas,Thetas,c_Thetas,Vint, + RnSPEs,RnV1s,RnV0s,RSPEs,RV1s,RV0s, + oSPEs,oV1s,oV0s, + idx_s_from_i,facs, + lnJ,tJNs,HNmats,AllTDs,Allinfos, + tMat,tildH,Erefs,errors,to) + nr = length(Ts) # num_replica + llhs = [-1.e+30 for i =1:nr] + nllhs = [-1.e+10 for i =1:nr] + Accept = false; TAcchit=0 + Acchits = zeros(Int64,nr) + xi = copy(Thetas[1]) + xj = copy(Thetas[1]) + tnum_states = sum([tJNs[i][2] for i =1:lnJ]) + Ens = [ Float64[] for i=1:tnum_states] + Evals = deepcopy(Erefs) + illhs = [ Float64[] for i=1:nr] # llh history for lowest temperature + offs = [0];for i = 2:lnJ;push!(offs,sum([tJNs[j][2] for j =1:i-1]));end + + # for destructive_eigen my_eigvals + Dim = size(tMat)[1] + v0 = [randn() for i=1:Dim] + + vks = [ zeros(Float64,Dim) for i = 1:150 ] + for i=1:Dim; vks[1][i] = randn(); end + tmp = dot(vks[1],vks[1]) + vks[1] .*= 1.0/sqrt(tmp) + num_history = 2 + en_s = [ [ zeros(Float64,tmp[2]) for i =1:num_history] for tmp in tJNs ] + #### + + jobs = Vector{Int64}[ ] + for jidx = 1:lnJ + info = Allinfos[jidx] + for idx = 1:length(AllTDs[jidx]) + i,j,Nij = info[idx] + push!(jobs,[jidx,idx,i,j,Nij]) + end + end + + for itM = 1:itnum_MCMC + if itM % div(itnum_MCMC,10) ==0; println("it ",itM);end + @inbounds for ridx = 1:nr + nSPEs = RnSPEs[ridx]; nV1s = RnV1s[ridx]; nV0s = RnV0s[ridx] + SPEs = RSPEs[ridx] ; V1s = RV1s[ridx] ; V0s = RV0s[ridx] + Theta = Thetas[ridx]; c_Theta = c_Thetas[ridx] + proposal_Vint!(nSPEs,nV1s,nV0s,SPEs,V1s,V0s,var_M*(ridx^(1/3))) + evalVsnt!(c_Theta,nSPEs,nV1s,nV0s,idx_s_from_i,facs) + nllhs[ridx] = L2norm(nSPEs,nV1s,nV0s,oSPEs,oV1s,oV0s) + @timeit to "Hmats" @inbounds @threads for job in jobs + jidx,idx,mi,mj,Nij = job + Hmat = @views HNmats[jidx][1] + Linv = @views HNmats[jidx][2] + TD = @views AllTDs[jidx][idx] + info = @views Allinfos[jidx][idx] + singleH(TD,mi,mj,Nij,c_Theta,Hmat) + end + @timeit to "gem&diag" for jidx = 1:lnJ + num_ev_target = tJNs[jidx][2] + Hmat = @views HNmats[jidx][1] + Linv = @views HNmats[jidx][2] + Eval = @views Evals[jidx] + BLAS.gemm!('N','N',1.0,Linv,Hmat,0.0,tMat) + BLAS.gemm!('T','N',1.0,Linv,tMat,0.0,tildH) + + @timeit to "eigsolve" begin + #Eval .= Arpack.eigs(tildH,nev=num_ev_target,which=:SR,tol=1.e-6,maxiter=150)[1] + Eval .= real.(eigsolve(tildH,num_ev_target,:SR)[1])[1:num_ev_target] + end + nllhs[ridx] += L2_llh(Eval,Erefs[jidx],errors[jidx];T=Ts[ridx]) + if ridx == 1 + if itM > burnin + for i=1:num_ev_target + push!(Ens[offs[jidx]+i],Eval[i]) + end + end + end + end + Accept = false + if nllhs[ridx] > llhs[ridx] + Accept =true + elseif nllhs[ridx]-llhs[ridx] > log(rand()) + Accept = true + end + if Accept + llhs[ridx] = nllhs[ridx] + Theta .= c_Theta + SPEs .= nSPEs; V1s .= nV1s;V0s .= nV0s + evalVint!(Vint,SPEs,V1s,V0s) + Acchits[ridx] += 1 + end + if itM > burnin + if ridx == 1; push!(iThetas,copy(Vint)); end + push!(illhs[ridx],llhs[ridx]) + continue + end + if itM == burnin;Acchits .= 0;TAcchit=0;end + end + ### exchange + if itM <= burnin || nr == 1; continue ;end + ri = itM % 2 + 1 + for r1 in ri:2:nr-1 + r2 = r1 + 1 + Ti = Ts[r1]; Tj = Ts[r2] + xi .= Thetas[r1]; xj .= Thetas[r2] + llh_i = llhs[r1]; llh_j = llhs[r2] + logp = (Ti-Tj) * (llh_i/Tj - llh_j/Ti) + TAccept = false + if logp > 0.0 + TAccept=true + elseif log(rand()) < logp + TAccept=true + end + if TAccept && itM > burnin + Thetas[r1] .= xj; Thetas[r2] .= xi + illhs[r1][end] = llh_j + illhs[r2][end] = llh_i + TAcchit += 1 + if r1 ==1 + evalVint!(Vint,RSPEs[r2],RV1s[r2],RV1s[r2]) + iThetas[end] .= Vint + end + end + end + end + rates = (100/(itnum_MCMC-burnin)) .* Acchits + tx = "Accept rates " + for i = 1:nr + tx *= " "*@sprintf("%5.1f ", rates[i]) + end + tx *= "\t exchange rate " *@sprintf("%5.1f ", 100*TAcchit/((itnum_MCMC-burnin)*nr)) + return tx, illhs, Ens +end + +function my_eigvals!(A,T,Dim,vks,en,num_ev_target;TF=[false],tol=1.e-6) + T .= 0.0 + for i=1:length(en); en[i] .= 1.e+10;end + for i=1:Dim; vks[1][i] = randn(); end + tmp = 1.0 / sqrt(dot(vks[1],vks[1])) + vks[1] .*= tmp + + for it = 1:length(vks)-1 + v = vks[it]; Hv = vks[it+1] + mul!(Hv,A,v) + talpha = dot(v,Hv) + T[it,it] = talpha + #diagonalize_T!(it,num_ev_target,T,en,2,TF,tol) + #print_vec("HNLanczos @$it ",en[1]) + if TF[1];break;end + axpy!(-talpha,v,Hv) + svks = @views vks[1:it-1] + ReORTH(it,Hv,svks) + tbeta = sqrt(dot(Hv,Hv)) + tmp = 1.0/tbeta + Hv .*= tmp + T[it+1,it] = tbeta; T[it,it+1] = tbeta + end + return nothing +end + +function proposal_Theta!(c_Theta,Theta,ln,var_M) + @inbounds for i=1:ln + c_Theta[i] = Theta[i] + var_M .* randn() + end + return nothing +end +function L2norm(nSPEs,nV1s,nV0s,oSPEs,oV1s,oV0s;lam=1e1) + s = 0.0 + for i =1:length(nSPEs); s += (nSPEs[i]-oSPEs[i])^2; end + for i =1:length(nV1s); s += (nV1s[i]-oV1s[i])^2; end + for i =1:length(nV0s); s += (nV0s[i]-oV0s[i])^2; end + return - 0.5 * s * lam +end +function L2_llh(Eval,Eref,err;T=1.0,sigma_sm=0.126) + s=0.0 + @inbounds for i = 1:length(Eval) + #s -= (Eval[i] - Eref[i])^2 / (sigma_sm^2 + err[i]^2) #i-dependent + #s -= (Eval[i] - Eref[i])^2 # i-independent err=1.0 + s -= (Eval[i] - Eref[i])^2 / (1.0+err[i]^2) # i-independent err=1.0 + end + return 0.5 * s / (T*length(Eval)) +end + +function read_exact(target_nuc,targetJ,sntf,lines) + hit = 0 + for (i,line) in enumerate(lines) + if occursin("$target_nuc",line) + if parse(Int,split(lines[i+1],"targetJ")[end]) == targetJ + hit = 1 + continue + end + end + if hit ==0;continue;end + csntf = split(sntf,"/")[end] + if occursin(csntf,line) + oJs = split(split(split(lines[i+1],"[")[2],"]")[1],",") + oEs = split(lines[i+2])[2:end] + Es=[] + for (idx,J) in enumerate(oJs) + tJ = Int(2*parse(Float64,J)) + if tJ != targetJ; continue;end + push!(Es,parse(Float64,oEs[idx])) + end + return Es + end + end + return [] +end + +function plot_EC_scatter(target_nuc,Hs,sumV,tJNs,Dims,exlines) + js = [ tJNs[i][1] for i=1:length(tJNs)] + xs = [ Float64[] for i=1:length(js)] + ys = [ Float64[] for i=1:length(js)] + minmax=Float64[1.e+10,-1.e+10] + yerrs = Float64[] + for tmp in sumV + sntf,tJ,EnsEC = tmp + Ens = read_exact(target_nuc,tJ,sntf,exlines) + tl = Float64[] + idx = 0 + for k = 1:length(js) + if js[k]==tJ + idx = k;break + end + end + try + Ens[1]; EnsEC[1] + catch + continue + end + for i=1:1 # only the yrast state + #if EnsEC[i] < Ens[i] + # println("variational err ", + # " EC ",EnsEC[i], " Exact ",Ens) + #end + push!(xs[idx],Ens[i]) + push!(ys[idx],EnsEC[i]) + push!(tl,Ens[i]) + push!(tl,EnsEC[i]) + push!(yerrs,(EnsEC[i] - Ens[i])/abs(Ens[i])) + end + if minimum(tl) < minmax[1] + minmax[1] = minimum(tl) + end + if maximum(tl) > minmax[2] + minmax[2] = maximum(tl) + end + end + if length(Hs) > 1; println(" typical error ", mean(yerrs));end + ##scatter + plt=pyimport("matplotlib.pyplot") + cm=pyimport("matplotlib.cm") + patches=pyimport("matplotlib.patches") + fig = plt.figure(figsize=(6,4)) + ax = fig.add_subplot(111); #axB = fig.add_subplot(212) + ax.set_xlabel("Exact (MeV)") + ax.set_ylabel("EC estimate (MeV)") + cols = ["red","blue","darkgreen","orange","purple","magenta","cyan"] + tms = ["o","^",",","d","p","h","8","v","*","1"] + ax.plot([minmax[1]-5,minmax[2]+5],[minmax[1]-5,minmax[2]+5], + linestyle="dotted",color="gray",alpha=0.6) + for jidx = 1:length(js) + nth = 1 #yrast only + tl = nothing + if nth == 1 + tl = "J="*latexstring(string(J2_to_J(js[jidx]))*"_{"*string(nth)*"}") + end + ax.scatter(xs[jidx],ys[jidx],marker=tms[jidx],s=20,lw=0.8, + zorder=150,color=cols[jidx],facecolor="none",label=tl,alpha=0.5) + end + Anum = parse(Int64, match(reg,target_nuc).match) + nuc_latex = latexstring("{}^{$Anum}")*split(target_nuc,string(Anum))[1] + + tl = " $nuc_latex ("*latexstring("N_s")*"="*string(Dims[1])*")" + ax.text(0.05, 0.9, tl,fontsize=15, transform=ax.transAxes) + ax.legend(loc="lower right",fontsize=12) + plt.savefig("./pic/scatter_ECestimates_"*target_nuc*".pdf",bbox_iinches="tight",pad_inches=0) + plt.close() +end + + +function ngauss(x,mu,varV) + return 1.0/(sqrt(2*pi*varV)) * exp( - 0.5* (x-mu)^2 / varV ) +end + +struct Tridx + i::Int + j::Int + Nth::Int +end + +function calcOBTD(OBTDs,idxs,p_sps,n_sps, + msps_p,msps_n, + tdims, + jocc_p,jocc_n, + pbits,nbits, + wfs) + lps=length(p_sps);lns=length(n_sps);lblock = length(nbits) + @inbounds for k=1:length(idxs) + #for k=1:length(idxs) + tmp = idxs[k] + i = tmp.i; j=tmp.j; Nij=tmp.Nth + w_i=wfs[i]; w_j = wfs[j] + pOBTD = OBTDs[Nij][1] + nOBTD = OBTDs[Nij][2] + @inbounds for bi = 1:lblock + idim = tdims[bi] + jocc_p_bi = jocc_p[bi] + jocc_n_bi = jocc_n[bi] + pbit = pbits[bi] + nbit=nbits[bi] + l_Nn = length(nbit) + offset = idim - l_Nn + @inbounds for (pidx,pocc) in enumerate(jocc_p_bi) + tM = offset + pidx*l_Nn + @inbounds for (nidx,nocc) in enumerate(jocc_n_bi) + Mi = tM + nidx + wfprod = w_i[Mi] .* w_j[Mi] + @inbounds for i = 1:lps + pOBTD[i] += wfprod .* pocc[i] + end + @inbounds for i = 1:lns + nOBTD[i] += wfprod .* nocc[i] + end + #pOBTD .+= wfprod .* pocc;nOBTD .+= wfprod .* nocc + #axpy!(wfprod,pocc,pOBTD);axpy!(wfprod,nocc,nOBTD) + end + end + end + end + if false + println("OBTDs: ") + for idx in idxs + i = idx.i; j=idx.j; Nij=idx.Nth + if i == j + println("_[$i] Nij $Nij") + println(" proton: ",OBTDs[Nij][1]) + println("neutron: ",OBTDs[Nij][2]);println("") + end + end + end + return nothing +end + +struct T1ifc + nth::Int64 + bi::Int64 + i::Int64 + f::Int64 + fac::Float64 +end +struct T0ifc + i::Int64 + f::Int64 + fac::Float64 +end + +function prepTBTD(tJ,idxs,p_sps,n_sps,msps_p::Array{Array{Int64,1},1},msps_n::Array{Array{Int64,1},1}, + pbits,nbits,labels,TBMEs,wfs::Array{Array{Float64,1}},tdims,Mps) + lblock=length(pbits) + lp = length(msps_p); ln = length(msps_n) + mstates = [msps_p,msps_n] + sps = [p_sps,n_sps] + loffs = [ 0, length(p_sps)] + TBTD1 = [ T1ifc[ ] , T1ifc[] ] + n = nthreads() + #@threads + for vrank =1:2 #pp:1, nn:2 + loff = loffs[vrank] + T1info = TBTD1[vrank] + vecs= [ [ [ false for i = 1:lp] for j=1:2], + [ [ false for i = 1:ln] for j=1:2]] + @inbounds for (i,ME) in enumerate(TBMEs[vrank]) + a,b,c,d,totJ = labels[vrank][i] + J2 = 2*totJ + ja,ma_s,ma_idxs = possible_mz(sps[vrank][a-loff],mstates[vrank]) + jb,mb_s,mb_idxs = possible_mz(sps[vrank][b-loff],mstates[vrank]) + jc,mc_s,mc_idxs = possible_mz(sps[vrank][c-loff],mstates[vrank]) + jd,md_s,md_idxs = possible_mz(sps[vrank][d-loff],mstates[vrank]) + @inbounds for (ic,mc) in enumerate(mc_s) + @inbounds for (id,md) in enumerate(md_s) + if c == d && mc >= md; continue;end + if abs(mc + md) > J2; continue;end + M_ani = mc + md + initialize_tvec!(vecs[vrank][1]); vecs[vrank][1][mc_idxs[ic]] = true + bit_c = bitarr_to_int(vecs[vrank][1]) + initialize_tvec!(vecs[vrank][1]); vecs[vrank][1][md_idxs[id]] = true + bit_d = bitarr_to_int(vecs[vrank][1]) + @inbounds for (ia,ma) in enumerate(ma_s) + @inbounds for (ib,mb) in enumerate(mb_s) + if a == b && ma>=mb; continue;end + if ma + mb != M_ani;continue;end + initialize_tvec!(vecs[vrank][2]);vecs[vrank][2][ma_idxs[ia]] = true + bit_a = bitarr_to_int(vecs[vrank][2]) + initialize_tvec!(vecs[vrank][2]);vecs[vrank][2][mb_idxs[ib]] = true + bit_b = bitarr_to_int(vecs[vrank][2]) + CG1 = clebschgordan(Float64,ja//2,ma//2,jb//2,mb//2,J2//2,M_ani//2) + CG2 = clebschgordan(Float64,jc//2,mc//2,jd//2,md//2,J2//2,M_ani//2) + bitlist = bit2b(bit_a,bit_b,bit_c,bit_d) + coeff = sqrt( (1.0+deltaf(a,b)) *(1.0+deltaf(c,d)) ) * CG1 * CG2 + if vrank==1 # pp + @inbounds for bi = 1:lblock + TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1] + idim = tdims[bi] + l_Nn = length(nbits[bi]) + @inbounds for (Npi,Phi) in enumerate(pbits[bi]) + TF_connectable(Phi,bitlist,TF) + if TF[1]==false; continue;end + calc_phase!(Phi,bitlist,lp,ret) + bisearch!(pbits[bi],ret[2],ridx) + Npf = ridx[1] + fac = coeff .*ret[1] + push!(T1info,T1ifc(i,bi,Npi,Npf,fac)) + end + end + else ## nn + @inbounds for bi = 1:lblock + TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1] + idim = tdims[bi] + l_Nn = length(nbits[bi]) + @inbounds for (Nni,Phi) in enumerate(nbits[bi]) + TF_connectable(Phi,bitlist,TF) + if TF[1]==false; continue;end + calc_phase!(Phi,bitlist,ln,ret) + bisearch!(nbits[bi],ret[2],ridx) + Nnf = ridx[1] + fac = coeff .*ret[1] + push!(T1info,T1ifc(i,bi,Nni,Nnf,fac)) + end + end + end + end + end + end + end + end + end + + ## pn + vrank=3 + loff = loffs[2] # for neutron + delMs = Int64[] + for j = 1:lp + for i = 1:lp + push!(delMs,msps_p[i][5]-msps_p[j][5]) + end + end + unique!(delMs);sort!(delMs,rev=true) + maxDeltaM = maximum(delMs) + bfs = [ Int64[ ] for i = 1:lblock] + bis = [ Int64[ ] for i = 1:lblock] + for i = 1:lblock + for j = i:lblock + if Mps[j] - Mps[i] in delMs + push!(bfs[i],j); push!(bis[j],i) + end + end + end + + vec_ani_p = [false for i = 1:lp];vec_cre_p = [false for i = 1:lp] + vec_ani_n = [false for i = 1:ln];vec_cre_n = [false for i = 1:ln] + retM = [0,0,0] + bif_idxs= Vector{Int64}[ ] + for i = 1:lblock + for j=i:lblock + push!(bif_idxs,[i,j,idx_from_ij(i,j,lblock)]) + end + end + pjumps = [ [ [ [T0ifc(0,0,0.0)] ] ] for i=1:length(TBMEs[vrank]) ] + njumps = [ [ [ [T0ifc(0,0,0.0)] ] ] for i=1:length(TBMEs[vrank]) ] + for i = 1:length(TBMEs[vrank]) + deleteat!(pjumps[i],1) + deleteat!(njumps[i],1) + end + for (nth,ME) in enumerate(TBMEs[vrank]) + a,b,c,d,totJ,vidx = labels[vrank][nth] + J2 = 2*totJ + ja,ma_s,ma_idxs = possible_mz(p_sps[a],msps_p) + jc,mc_s,mc_idxs = possible_mz(p_sps[c],msps_p) + jb,mb_s,mb_idxs = possible_mz(n_sps[b-loff],msps_n) + jd,md_s,md_idxs = possible_mz(n_sps[d-loff],msps_n) + for (ic,mc) in enumerate(mc_s) + initialize_tvec!(vec_ani_p); vec_ani_p[mc_idxs[ic]] = true + bit_c = bitarr_to_int(vec_ani_p) + for (ia,ma) in enumerate(ma_s) + initialize_tvec!(vec_cre_p); vec_cre_p[ma_idxs[ia]] = true + bit_a = bitarr_to_int(vec_cre_p) + Mp = ma - mc + bisearch!(delMs,Mp,retM);idx = retM[1] + for (id,md) in enumerate(md_s) + if abs(mc + md) > J2; continue;end + initialize_tvec!(vec_ani_n); vec_ani_n[md_idxs[id]] = true + bit_d = bitarr_to_int(vec_ani_n) + for (ib,mb) in enumerate(mb_s) + if mb - md + Mp != 0; continue;end + initialize_tvec!(vec_cre_n); vec_cre_n[mb_idxs[ib]] = true + bit_b = bitarr_to_int(vec_cre_n) + CG1 = clebschgordan(Float64,ja//2,ma//2,jb//2,mb//2, + J2//2,(ma+mb)//2) + CG2 = clebschgordan(Float64,jc//2,mc//2,jd//2,md//2, + J2//2,(mc+md)//2) + fac = CG1*CG2 #* sqrt(tJ+1) /sqrt(J2+1.0) + pjump = [ T0ifc[] for i=1:length(bif_idxs)] + njump = [ T0ifc[] for i=1:length(bif_idxs)] + TF=[true]; ret_p=[0,0];ret_n=[0,0] + ridx=[0,0,0] + + for bi = 1:lblock + idim = tdims[bi] + l_Nn_i = length(nbits[bi]) + for bfidx = 1:length(bfs[bi]) + bf = bfs[bi][bfidx] + deltaM = Mps[bf] - Mps[bi] + tfac = fac * ifelse(bi==bf,0.5,1.0) + if (deltaM in delMs) == false;continue;end + fdim = tdims[bf] + l_Nn_f = length(nbits[bf]) + bif_idx = idx_from_ij(bi,bf,lblock) + for (Npi,pPhi) in enumerate(pbits[bi]) + TF_connectable_1(pPhi,bit_a,bit_c,TF) + if TF[1]==false; continue;end + calc_phase_1!(pPhi,bit_a,bit_c,ret_p) + bisearch!(pbits[bf],ret_p[2],ridx) + if ridx[1] == 0;continue;end + Npf = ridx[1] + push!(pjump[bif_idx], + T0ifc(Npi,Npf,tfac*ret_p[1])) + end + for (Nni,nPhi) in enumerate(nbits[bi]) + TF_connectable_1(nPhi,bit_b,bit_d,TF) + if TF[1]==false; continue;end + calc_phase_1!(nPhi,bit_b,bit_d,ret_n) + bisearch!(nbits[bf],ret_n[2],ridx) + if ridx[1] == 0;continue;end + Nnf = ridx[1] + push!(njump[bif_idx], + T0ifc(Nni,Nnf,1.0*ret_n[1])) + end + # if ma==mb==mc==md==-1 + # if a==1 && b==3 && c==1 && ; + # d==3 && totJ==1 && bi==bf + # println("V($a$b$c$d$totJ) hit ", + # labels[vrank][nth], + # " ms $ma $mb $mc $md ", + # " bi $bi bf $bf bifidx $bif_idx") + # println("pjump ", pjump[bif_idx], + # " njump ", njump[bif_idx]) + # println("") + # end + # end + end + end + push!(pjumps[nth],pjump) + push!(njumps[nth],njump) + end + end + end + end + end + return TBTD1, pjumps,njumps,bif_idxs +end + +function calcTBTD(TBTDs, + opTBTD1,pjumps,njumps, + pbits,nbits,tdims,wfs,idxs,bif_idxs, + olabels,oTBMEs,labels,to) + @inbounds @threads for k in eachindex(idxs) + tmp = idxs[k] + i = tmp.i; j=tmp.j; Nij=tmp.Nth + w_i=wfs[i]; w_j = wfs[j] + coeffs = [ zeros(Float64,length(oTBMEs)) for vrank=1:3] + #@timeit to "pp/nn" begin + ## pp/nn + vrank =1 + for tmp in opTBTD1[vrank] ## pp + nth = tmp.nth; bi = tmp.bi + nbit = nbits[bi] + pbit = pbits[bi] + lN = length(nbit) + vidx = labels[vrank][nth][end] + Npi = tmp.i;Npf = tmp.f;fac = tmp.fac + tMi = tdims[bi]+ (Npi-1)*length(nbit) + tMf = tdims[bi]+ (Npf-1)*length(nbit) + @inbounds for nidx = 1:length(nbit) + Mi = tMi+nidx; Mf = tMf+nidx + coeffs[vrank][vidx] += fac .* (w_i[Mf].*w_j[Mi]) + end + end + vrank =2 + for tmp in opTBTD1[vrank] ## nn + nth = tmp.nth; bi = tmp.bi + nbit = nbits[bi] + pbit = pbits[bi] + lN = length(nbit) + vidx = labels[vrank][nth][end] + Nni = tmp.i;Nnf = tmp.f;fac = tmp.fac + tMi = tdims[bi]+ Nni -lN + tMf = tdims[bi]+ Nnf -lN + @inbounds for pidx = 1:length(pbit) + Mi = tMi + pidx*lN; Mf = tMf + pidx*lN + coeffs[vrank][vidx] += fac .* (w_i[Mf].*w_j[Mi]) + end + end + #@timeit to "pn" begin + ### pn + vrank=3 + @inbounds for (nth,label) in enumerate(labels[vrank]) + vidx = label[end] + pjump = pjumps[nth] + njump = njumps[nth] + @inbounds for opidx in eachindex(pjump) + tmp_pj = pjump[opidx] + tmp_nj = njump[opidx] + @inbounds for bidx in eachindex(bif_idxs) + bi,bf,dummy = bif_idxs[bidx] + idim = tdims[bi]; lNi = length(nbits[bi]) + fdim = tdims[bf]; lNf = length(nbits[bf]) + pjs = tmp_pj[bidx] + njs = tmp_nj[bidx] + @inbounds for pj in pjs + Npi = pj.i;Npf = pj.f; fac_p = pj.fac + tMi = idim + (Npi-1)*lNi + tMf = fdim + (Npf-1)*lNf + @inbounds for nj in tmp_nj[bidx] + Nni = nj.i; Nnf=nj.f; fac_n = nj.fac + Mi = tMi + Nni; Mf = tMf + Nnf + coeffs[vrank][vidx] += fac_p .* fac_n .* (w_i[Mf].*w_j[Mi]) + coeffs[vrank][vidx] += fac_p .* fac_n .* (w_i[Mi].*w_j[Mf]) + end + end + end + end + end + for k = 1:3 + TBTDs[Nij] .+= coeffs[k] + end + end + return nothing +end diff --git a/src/ShellModel/lanczos_methods.jl b/src/ShellModel/lanczos_methods.jl index 31b855a7..68d14108 100644 --- a/src/ShellModel/lanczos_methods.jl +++ b/src/ShellModel/lanczos_methods.jl @@ -134,7 +134,6 @@ function Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj, ReORTH(it,vkp1,svs) beta = sqrt(dot(vkp1,vkp1)) - println(" betaJ $beta") #if beta < 1.e-4;eljit=it;TF[1]=true;JTF[1]=true;break;end vkp1 .*= 1.0/beta Jmat[it+1,it] = beta; Jmat[it,it+1] = beta diff --git a/src/ShellModel/lanczos_methods.jl~ b/src/ShellModel/lanczos_methods.jl~ new file mode 100644 index 00000000..31b855a7 --- /dev/null +++ b/src/ShellModel/lanczos_methods.jl~ @@ -0,0 +1,871 @@ +function TRL(vks,uks,Tmat,k,pbits,nbits,jocc_p,jocc_n, + SPEs,pp_2bjump,nn_2bjump,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,checkorth=true,debugmode=true) + + mdim = tdims[end] + TF=[false] + + lnJ = ls + Jvs = [zeros(Float64,mdim) for i=1:lnJ] + Jvret = [zeros(Float64,mdim)] + Jmat = zeros(Float64,lnJ,lnJ) + Jvs[1] .= vks[1] + Jtol = 1.e-7 + JTF = [false];beta_J = 1.0 + + ### + num_restart = 0 + Jvlog_before = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + Jvlog_after = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + Hvs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + HvRs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + ### + + if doubleLanczos + Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to) + vks[1] .= Jvs[1] + vks[1] .*= 1.0 / sqrt(dot(vks[1],vks[1])) + end + + elit = 1 + while TF[1]==false + for it = k:lm-1 + vk =vks[it]; vkp1 =vks[it+1] + @timeit to "operate H" begin + operate_H!(vk, vkp1, pbits,nbits,jocc_p,jocc_n,SPEs, + pp_2bjump,nn_2bjump,tdims,bis,bfs,block_tasks, + p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to) + end + if debugmode + Hvs[it] .= vkp1 + end + talpha = dot(vk,vkp1) + Tmat[it,it] = talpha + 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] + @timeit to "ReORTH" ReORTH(it,vkp1,svks) + if debugmode + HvRs[it] .= vkp1 + end + tbeta = sqrt(dot(vkp1,vkp1)) + vkp1 .*= 1.0/tbeta + #if checkorth; Check_Orthogonality(it,vks,en); end + if debugmode + print_vec("En TRL $it ",en[1]) + println(" beta $tbeta") + end + Tmat[it+1,it] = tbeta; Tmat[it,it+1] = tbeta + + if doubleLanczos + if tbeta < Jtol;TF[1]=true;elit=it;break;end + @timeit to "JJ lanczos" begin + for i = 2:length(Jvs);Jvs[i] .=0.0;end + Jmat .= 0.0;JTF[1] = false + Jvs[1] .= vkp1 + if debugmode; Jvlog_before[it] .= vkp1; end ## + Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to) + vkp1 .= Jvs[1] + if debugmode; Jvlog_after[it] .= vkp1; end ## + if TF[1];elit=it;break;end + end + #println(" beta $tbeta beta_inJ $beta_in_J") + sv = @view vks[1:it] + ReORTH(it,vkp1,sv) + vkp1 .*= 1.0 / sqrt(dot(vkp1,vkp1)) + end + end + if TF[1] == false + @timeit to "Restart" begin + num_restart += 1 + if debugmode; write_wf_hdf5(vks, Hvs, HvRs,Tmat,Jvlog_before, Jvlog_after, "JJ_TR_res$(num_restart).h5", false); end + ThickRestart(vks,uks,Tmat,lm,ls,debugmode) + end + end + k = ls+1 + end + return elit +end + +function write_wf_hdf5(vks,Hvs,HvRs,Tmat,Jv_bef, Jv_aft, fname::String, is_block) + # for debug + io = h5open(fname, "w") + mdim = length(vks[1]) + write(io,"len", length(vks)) + write(io,"mdim", mdim) + for i = 1:length(vks) + write(io, "Hvs_$i", Hvs[i][:]) + write(io,"HvRs_$i", HvRs[i][:]) + end + write(io, "Tmat",Tmat) + for i = 1:length(Jv_bef) + write(io, "Jv_bef_$i", Jv_bef[i]) + end + for i = 1:length(Jv_aft) + write(io, "Jv_aft_$i", Jv_aft[i]) + end + close(io) + return nothing +end + +function Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj, + Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to) + mdim = tdims[end] + lnJ = length(Jvs) + eljit=1; k = 1; inow=k + beta = 0.0 + while JTF[1] == false + for it = k:lnJ-1 + inow = it + vk = Jvs[it]; vkp1 = Jvs[it+1] + operate_J!(vk,vkp1,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd,beta_J) + axpy!(eval_jj,vk,vkp1)## vkp1 .+= eval_jj .* vk + Jmat[it,it] = dot(vk,vkp1) + teval = eigvals(@views Jmat[1:it,1:it])[1] + #println("JJ it ", @sprintf("%3i",it), " eval ", @sprintf("%12.8f",teval)) + if abs(teval-eval_jj) < Jtol + eljit=it;JTF[1] = true;break + end + + svs = @views Jvs[1:it] + ReORTH(it,vkp1,svs) + + beta = sqrt(dot(vkp1,vkp1)) + println(" betaJ $beta") + #if beta < 1.e-4;eljit=it;TF[1]=true;JTF[1]=true;break;end + vkp1 .*= 1.0/beta + Jmat[it+1,it] = beta; Jmat[it,it+1] = beta + eljit = it + end + if JTF[1]==false + ThickRestart_J(Jvs,Jvret,Jmat,eljit+1,1,eval_jj,Jtol) + k=2 + end + end + if inow > k + ThickRestart_J(Jvs,Jvret,Jmat,eljit+1,1,eval_jj,Jtol) + end + return nothing +end + +function operate_H!(wf,twf,pbits,nbits,jocc_p,jocc_n,SPEs,pp_2bjump,nn_2bjump,tdims,bis,bfs,block_tasks,p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to=nothing) + @inbounds @threads for bi in block_tasks + if bi==0; continue;end #empty job + ret = [0,0,0]; ret2 = [0,0,0] + idim = tdims[bi] + l_Np = length(pbits[bi]) + l_Nn = length(nbits[bi]) + offset = idim -l_Nn + #@timeit to "pn1" begin + for (bfidx,bf) in enumerate(bfs[bi]) + bisearch!(delMs,Mps[bf]-Mps[bi],ret) #!! bf = bfs[bi][j] + Vs = Vpn[ret[1]] + fdim = tdims[bf] + l_Nn_f = length(nbits[bf]) + p_NiNf = p_NiNfs[bi][bfidx] + n_NiNf = n_NiNfs[bi][bfidx] + off_f = fdim-l_Nn_f + for (nth,V) in enumerate(Vs) + Npifs = p_NiNf[nth]; Nnifs = n_NiNf[nth] + for Npif in Npifs + tMi = offset+ Npif.i * l_Nn + tMf = off_f + Npif.f * l_Nn_f + phase_p = Npif.phase + for Nnif in Nnifs + Mi = tMi + Nnif.i; Mf = tMf + Nnif.f + phase_n = Nnif.phase + twf[Mi] += ifelse(phase_p!=phase_n,-V,V) * wf[Mf] + end + end + end + end + #@timeit to "pn2" begin + for j = 1:length(bis[bi])-1 #### tbf = bi !!! + tbi = bis[bi][j] + bisearch!(delMs,Mps[bi]-Mps[tbi],ret) + bisearch_ord!(bfs[tbi],bi,ret2) + fdim=tdims[tbi] + l_Nn_i=length(nbits[tbi]) + p_NiNf = p_NiNfs[tbi][ret2[1]] + n_NiNf = n_NiNfs[tbi][ret2[1]] + off_f = fdim - l_Nn_i + for (nth,V) in enumerate(Vpn[ret[1]]) + Npifs = p_NiNf[nth] + Nnifs = n_NiNf[nth] + for Npif in Npifs + tMi = off_f + Npif.i *l_Nn_i #idim <-> fdim + tMf = offset + Npif.f*l_Nn + phase_p = Npif.phase + for Nnif in Nnifs + Mi = tMi + Nnif.i; Mf = tMf + Nnif.f + phase_n = Nnif.phase + twf[Mf] += ifelse(phase_p!=phase_n,-V,V) * wf[Mi] + end + end + end + end + #@timeit to "pp/nn" begin + ### pp/nn interaction + for (Npi,tinfo) in enumerate(pp_2bjump[bi]) + tMi = offset + Npi*l_Nn + for (jj,tmp) in enumerate(tinfo) + tMf = offset + l_Nn*tmp.f + fac = tmp.coef + for nidx = 1:l_Nn + Mi = tMi + nidx; Mf = tMf + nidx + twf[Mf] += fac * wf[Mi] + twf[Mi] += fac * wf[Mf] + end + end + end + for (Nni,tinfo) in enumerate(nn_2bjump[bi]) + tMi = offset + Nni + for (jj,tmp) in enumerate(tinfo) + tMf = offset + tmp.f + fac = tmp.coef + for pidx = 1:l_Np + Mi = tMi + pidx*l_Nn; Mf = tMf + pidx*l_Nn + twf[Mf] += fac * wf[Mi] + twf[Mi] += fac * wf[Mf] + end + end + end + + #@timeit to "1b" begin + ### one-body operator + for pidx = 1:l_Np + tMi = offset + pidx*l_Nn + for nidx =1:l_Nn + Mi = tMi + nidx + twf[Mi] += (dot(SPEs[1],jocc_p[bi][pidx])+ + dot(SPEs[2],jocc_n[bi][nidx])) * wf[Mi] + end + end + end + return nothing +end + +function diagonalize_T!(k::Int64,num_ev::Int64,Tmat,en::Array{Array{Float64,1}},num_history::Int64, TF::Array{Bool,1}, tol::Float64) + for ith = num_history:-1:2 + en[ith] .= en[ith-1] + end + n = minimum([num_ev,k]) + @views en[1][1:n] .= eigvals(@views Tmat[1:k,1:k] )[1:n] + if all( en[2]-en[1] .< tol) ; TF[1]=true; end + return nothing +end + +function ThickRestart(vks,uks,Tmat,lm,ls,debugmode=false,make_sure=true) + vals,vecs = eigen(@views Tmat[1:lm-1,1:lm-1]) + Tmat_before = ifelse(debugmode,copy(Tmat),[0.0]) + beta = Tmat[lm-1,lm] + Tmat .= 0.0 + @inbounds for k = 1:ls + Tmat[k,k] = vals[k] + end + @inbounds for (k,uk) in enumerate(uks) + tmp = beta .* vecs[lm-1,k] + Tmat[ls+1,k] = tmp; Tmat[k,ls+1] = tmp + uk .= 0.0 + for j=1:lm-1 + axpy!(vecs[j,k],vks[j],uk) + end + end + @inbounds for (k,uk) in enumerate(uks) + if make_sure && k > 1 + svs = @view uks[1:k-1] + ReORTH(k-1,uk,svs) + end + vks[k] .= uk .* (1.0 /sqrt(dot(uk,uk))) + end + vks[ls+1] .= vks[lm] + if make_sure + svs = @view vks[1:ls] + ReORTH(ls,vks[ls+1],svs) + end + for k = ls+2:lm + vks[k] .= 0.0 + end + # if debugmode + # io = h5open("restart.hdf5", "w") + # write(io,"lm", lm) + # write(io,"ls", ls) + # write(io, "beta", beta) + # write(io,"vals", vals) + # write(io,"vecs", vecs) + # for i = 1:ls + # write(io, "uks_$i", uks[i]) + # end + # write(io, "Tmat_before", Tmat_before) + # write(io, "Tmat", Tmat) + # close(io) + # end + return nothing +end + +""" +Thick-Restart Block Lanczos method +""" +function TRBL(q,vks,uks,Tmat,Beta_H,pbits,nbits,jocc_p,jocc_n,SPEs, + pp_2bjump,nn_2bjump,bis,bfs,block_tasks, + p_NiNfs,n_NiNfs,Mps,delMs,Vpn,tdims, + eval_jj,oPP,oNN,oPNu,oPNd,Jidxs, + num_ev,num_history,lm,ls_sub,en,tol,to; + doubleLanczos=false,verbose=false,debugmode=false) + ls = q*ls_sub + mdim = tdims[end] + if doubleLanczos && mdim < 100 + @warn "Block Lanczos with J projection is not stable for small systems.\nUse TRL (is_block=false) instead." + end + + TF=[false] + elit = 1 + inow = 1; itmin = 1; itmax = div(lm,q)-1 + rescount = 0 + Vt = zeros(Float64,q,mdim) + R = zeros(Float64,q,q) + Beta_J = zeros(Float64,q,q) + + lnJ = 20 + Jvs = [zeros(Float64,q,mdim) for i=1:lnJ] + Jvret = [zeros(Float64,q,mdim) ] + Jmat = zeros(Float64,lnJ*q,lnJ*q) + JTF = [false] + Jtol = 1.e-7 + U = zeros(Float64,mdim,lnJ*q) + Mat = zeros(Float64,mdim,q) + + Jvlog_before = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + Jvlog_after = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + Hvs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + HvRs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] + + if doubleLanczos + @timeit to "JJ lanczos" begin + Jvs[1] .= vks[1] + bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,to,Beta_H,U,Mat) + vks[1] .= Jvs[1];Beta_J .= 0.0;JTF[1]=false + for i=1:lnJ;Jvs[i] .= 0.0;end + V = vks[1] + end + end + num_restart = 0 + while TF[1]==false + for it = itmin:itmax + inow = it + V = vks[it] + HV = vks[it+1] + @timeit to "bl_operateH" begin + bl_operate_H!(q,V,HV,pbits,nbits, + jocc_p,jocc_n,SPEs,pp_2bjump,nn_2bjump, + tdims,bis,bfs,block_tasks, + p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to) + end + if debugmode + Hvs[it] .= HV[1,:] + end + BLAS.gemm!('N','T',1.0,V,HV,0.0,R) + if issymmetric(R) == false + @inbounds for i=1:q;for j=i:q + R[i,j] = R[j,i] + end;end + end + + @views Tmat[q*it-q+1:q*it,q*it-q+1:q*it] .= R + diagonalize_T!(it*q,num_ev,Tmat,en,num_history,TF,tol) + if debugmode + print_vec("En TRBL $it ",en[1]) + end + BLAS.gemm!('N','N',-1.0,R,V,1.0,HV)#mul!(HV,R,V,-1.0,1.0) + s_vks = @views vks[1:it-1] + @timeit to "ReORTH" bl_ReORTH(q,it,mdim,s_vks,HV,Vt,R) + if debugmode + HvRs[it] .= HV[1,:] + end + bl_QR!(HV',Beta_H,mdim,q)#myQR!(HV',Beta_H,mdim,q) + + if doubleLanczos + tnorm = norm(Diagonal(Beta_H),Inf) + if tnorm < Jtol + if verbose; println("Hbn norm $tnorm");end + TF[1]=true;elit=it;break + end + @timeit to "JJ lanczos" begin + for i=2:lnJ; Jvs[i] .= 0.0; end + JTF[1] = false; Jmat .= 0.0;Jvs[1] .= HV + Vt .= 0.0; R .= 0.0 + if debugmode; Jvlog_before[it] .= Jvs[1][1,:]; end + bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,to,Beta_H,U,Mat) + HV .= Jvs[1]; Beta_J .=0.0 + if debugmode; Jvlog_after[it] .= Jvs[1][1,:]; end + end + end + if debugmode + println(" R $R betaH $Beta_H diff $(R-Beta_H)") + end + add_bl_T!(q,it,Tmat,Beta_H) + if TF[1];elit = it;break;end + end + if TF[1] == false + @timeit to "bl_Restart" begin + num_restart += 1 + if debugmode; write_wf_hdf5(vks,Hvs,HvRs,Tmat,Jvlog_before, Jvlog_after,"JJ_TRBL_res$(num_restart).h5",debugmode); end + bl_ThickRestart(q,vks,uks,Beta_H,Tmat,inow,ls_sub,mdim,Vt,debugmode) + end + itmin = ls_sub + 1 + end + end + return elit +end + +function bl_ReORTH(q,it,mdim,vks,HV,Vt,R) + #LinearAlgebra.jl: mul!(C, A, B, α, β) -> C: ABalpha + Cbeta + #BLAS.gemm!(tA, tB, alpha, A, B,beta,C) C := alpha*A*B + beta*C + # :tA,tB 'N': normal 'T': transpose + @inbounds for l = it-1:-1:1 + Vt .= vks[l] + BLAS.gemm!('N','T', 1.0,HV,Vt,0.0,R) #mul!(R,HV,Vt') + BLAS.gemm!('N','N',-1.0,R,Vt,1.0,HV) #mul!(HV,R,Vt,-1.0,1.0) + end + R .= 0.0 + return nothing +end + +# wf:V, twf:HV +function bl_operate_H!(q,wf,twf,pbits,nbits,jocc_p,jocc_n, + SPEs,pp_2bjump,nn_2bjump,tdims,bis,bfs,block_tasks, + p_NiNfs,n_NiNfs,Vpn, Mps,delMs,to=nothing) + #lblock = length(pbits) + @inbounds @threads for bi in block_tasks + if bi==0; continue;end #empty job + ret = [0,0,0]; ret2 = [0,0,0] + idim = tdims[bi] + l_Np = length(pbits[bi]) + l_Nn = length(nbits[bi]) + offset = idim -l_Nn + #@timeit to "pn1" begin + @inbounds for (bfidx,bf) in enumerate(bfs[bi]) + bisearch!(delMs,Mps[bf]-Mps[bi],ret) #!! bf = bfs[bi][j] + Vs = Vpn[ret[1]] + fdim = tdims[bf] + l_Nn_f = length(nbits[bf]) + p_NiNf = p_NiNfs[bi][bfidx] + n_NiNf = n_NiNfs[bi][bfidx] + off_f = fdim-l_Nn_f + @inbounds for (nth,V) in enumerate(Vs) + Npifs = p_NiNf[nth]; Nnifs = n_NiNf[nth] + @inbounds for Npif in Npifs + tMi = offset+ Npif.i * l_Nn + tMf = off_f + Npif.f * l_Nn_f + phase_p = Npif.phase + @inbounds for Nnif in Nnifs + Mi = tMi + Nnif.i; Mf = tMf + Nnif.f + phase_n = Nnif.phase + coeff = ifelse(phase_p!=phase_n,-V,V) + w_f = @views twf[:,Mi] + w_i = @views wf[:,Mf] + @inbounds for b=1:q + w_f[b] += coeff * w_i[b] + end + end + end + end + end + #end + #@timeit to "pn2" begin + @inbounds for j = 1:length(bis[bi])-1 #### tbf = bi !!! + tbi = bis[bi][j] + bisearch!(delMs,Mps[bi]-Mps[tbi],ret) + bisearch_ord!(bfs[tbi],bi,ret2) + fdim=tdims[tbi] + l_Nn_i=length(nbits[tbi]) + p_NiNf = p_NiNfs[tbi][ret2[1]] + n_NiNf = n_NiNfs[tbi][ret2[1]] + off_f = fdim - l_Nn_i + @inbounds for (nth,V) in enumerate(Vpn[ret[1]]) + Npifs = p_NiNf[nth] + Nnifs = n_NiNf[nth] + @inbounds for Npif in Npifs + tMi = off_f + Npif.i *l_Nn_i #idim <-> fdim + tMf = offset + Npif.f*l_Nn + phase_p = Npif.phase + @inbounds for Nnif in Nnifs + Mi = tMi + Nnif.i; Mf = tMf + Nnif.f + phase_n = Nnif.phase + coeff = ifelse(phase_p!=phase_n,-V,V) + w_f = @views twf[:,Mf] + w_i = @views wf[:,Mi] + @inbounds for b=1:q + w_f[b] += coeff * w_i[b] + end + end + end + end + end + #end + #@timeit to "pp/nn" begin + ### pp/nn interaction + @inbounds for (Npi,tinfo) in enumerate(pp_2bjump[bi]) + tMi = offset + Npi*l_Nn + @inbounds for (jj,tmp) in enumerate(tinfo) + tMf = offset + l_Nn*tmp.f + fac = tmp.coef + @inbounds for nidx = 1:l_Nn + Mi = tMi + nidx; Mf = tMf + nidx + w_f1 = @views twf[:,Mf] + w_i1 = @views wf[:,Mi] + w_f2 = @views twf[:,Mi] + w_i2 = @views wf[:,Mf] + @inbounds for b=1:q + w_f1[b] += fac * w_i1[b] + w_f2[b] += fac * w_i2[b] + end + end + end + end + @inbounds for (Nni,tinfo) in enumerate(nn_2bjump[bi]) + tMi = offset + Nni + @inbounds for (jj,tmp) in enumerate(tinfo) + tMf = offset + tmp.f + fac = tmp.coef + @inbounds for pidx = 1:l_Np + Mi = tMi + pidx*l_Nn; Mf = tMf + pidx*l_Nn + w_f1 = @views twf[:,Mf] + w_i1 = @views wf[:,Mi] + w_f2 = @views twf[:,Mi] + w_i2 = @views wf[:,Mf] + @inbounds for b=1:q + w_f1[b] += fac * w_i1[b] + w_f2[b] += fac * w_i2[b] + end + end + end + end + #end + ### one-body operator + @inbounds for pidx = 1:l_Np + tMi = offset + pidx*l_Nn + @inbounds for nidx =1:l_Nn + Mi = tMi + nidx + coeff = (dot(SPEs[1],jocc_p[bi][pidx])+ + dot(SPEs[2],jocc_n[bi][nidx])) + w_f = @views twf[:,Mi] + w_i = @views wf[:,Mi] + @inbounds for b=1:q + w_f[b] += coeff * w_i[b] + end + end + end + end + return nothing +end + +function Jcompress(q,vks,Tmat,inow,ls_sub,mdim,R,bnout,U,Mat,Beta_J;use_LAPACK=false) + lm = q*inow; ls = q*ls_sub + vals = [0.0] + vecs = [0.0;0.0] + if use_LAPACK # used for only for debug (you need "l_diag.F90" by S.Yoshida) + M = Tmat[1:lm,1:lm] + ccall((:diagonalize_double_,"l_diag.so"),Nothing, + (Ref{Int64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Int64}), + lm,M,vals,vecs,ls) + else + M = @views Tmat[1:lm,1:lm] + vals,vecs = eigen(M) # eigen(Symmetric(M,:L)) + end + + tv = @views vecs[1:ls,1:ls] + BLAS.gemm!('T','N',1.0,tv,bnout,0.0,R)# mul!(R,tv,bnout) + bnout .= R + + Tmat .= 0.0; U .= 0.0 + for i=1:inow + bv = vks[i] + for b=1:q + j = q*(i-1)+b + u = @views U[:,j] + u .= @views bv[b,:] + end + end + tU = @views U[:,1:lm] + tv = @views vecs[1:lm,1:ls] + BLAS.gemm!('T','T',1.0,tv,tU,0.0,vks[1]) + Beta_J .= @views vecs[1:ls,1:ls] + + return nothing +end + +function like_bl_ThickRestart(vks,uks,Tmat,lm,ls,debugmode) + q = 1 + ls_sub = div(ls,q) + mdim = size(vks[1]) + R = beta = Tmat[lm-1,lm] + vals,vecs = eigen(@views Tmat[1:lm-1,1:lm-1]) + Tmat .= 0.0 + + for k = 1:ls; Tmat[k,k] = vals[k]; end + tv = @views vecs[lm-1,1:ls] + r = R .* tv + + Tmat[ls+1,1:ls] .= r + Tmat[1:ls,ls+1] .= r + + for k = 1:ls + uk = uks[k] #.= 0.0 + uk .= 0.0 + for j=1:lm-1 + vk = vks[j] + fac = vecs[j,k] + uk .+= fac .* vk + end + end + for i = 1:ls + vks[i] .= uks[i] + println("norm vks[$i] $(norm(vks[i]))") + end + vks[ls+1] .= vks[lm] + for i= ls_sub+2:length(vks) + vks[i] .= 0.0 + end + for i=2:ls+1 + v = vks[i] + for j = i-1:-1:1 + fac = dot(vks[j],v) + v .-= fac .* vks[j] + end + w = sqrt(dot(v,v)) + v .*= 1.0/w + println(" => norm vks[$i] $(norm(vks[i]))") + end + return nothing +end + +function bl_ThickRestart(q,vks,uks,R,Tmat,inow,ls_sub,mdim,Vt,debug_mode) + lm = q*inow + ls = q*ls_sub + r = zeros(Float64,q,ls) + vals,vecs = eigen(@views Tmat[1:lm,1:lm]) + Tmat_before = copy(Tmat) + + println("βm1 ",Tmat[lm-1,lm], " ", Tmat[lm-2,lm-1]) + + Tmat .= 0.0 + for k = 1:ls; Tmat[k,k] = vals[k]; end + tv = @views vecs[lm-q+1:lm,1:ls] + mul!(r,R,tv) #BLAS.gemm!('N','N',1.0,R,tv,0.0,r) + + println("R = beta? $R") + println("r in blockRestart $r") + println("tv $tv") + + Tmat[ls+1:ls+q,1:ls] .= r + Tmat[1:ls,ls+1:ls+q] .= r' + for bi = 1:ls_sub + uk = uks[bi] #.= 0.0 + uk .= 0.0 + for b=1:q + k = q*(bi-1) + b + for j=1:inow + vk = vks[j] + for bj = 1:q + v = @views vk[bj,:] + idx = q*(j-1) +bj + fac = vecs[idx,k] + @views u_vv = @views uk[b,:] + axpy!(fac, v, u_vv) + end + end + end + end + for i = 1:ls_sub + vks[i] .= uks[i] + end + vks[ls_sub+1] .= vks[inow+1] + for i= ls_sub+2:length(vks) + vks[i] .= 0.0 + end + for i=1:ls_sub+1 + V = vks[i] + if i>1 + for j = i-1:-1:1 + Vt .= vks[j] + mul!(R,V,Vt') + mul!(V,R,Vt,-1.0,1.0) + end + bl_QR!(V',R,mdim,q) + end + end + # if debug_mode + # io = h5open("restart_blocck.hdf5", "w") + # write(io,"lm", lm) + # write(io,"ls", ls) + # write(io,"vals", vals) + # write(io,"vecs", vecs) + # for i = 1:ls + # write(io, "uks_$i", uks[i][1,:]) + # end + # write(io, "Tmat_before", Tmat_before) + # write(io, "Tmat", Tmat) + # close(io) + # end + return nothing +end + +function bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj, + Jidxs,oPP,oNN,oPNu,oPNd,to,bnout,U,Mat;verbose=false) + mdim = tdims[end] + lnJ = length(Jvs) + itmin = 1;itmax = lnJ-1; inow = 0 + rescount = 0 + while JTF[1] == false + for it = itmin:itmax + inow = it + V = Jvs[it]; JV = Jvs[it+1] + bl_operate_J!(q,V,JV,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd) + axpy!(eval_jj,V,JV)##JV .+= eval_jj .* V + BLAS.gemm!('N','T',1.0,V,JV,0.0,R) + @inbounds for j=1:q;for i=1:j + R[j,i] = R[i,j] + end;end + @views Jmat[q*it-q+1:q*it,q*it-q+1:q*it] .= R + tJmat = @views Jmat[1:it*q,1:it*q] + + jeval = eigvals(tJmat)[1:q] + if verbose; print_vec("$it jeval=> ",jeval;long=true);end + if jeval[1] - eval_jj < -Jtol + if verbose;println("neg JJ @$it jeval");end + JTF[1] = true;break + end + if all( abs.(jeval.-eval_jj) .< Jtol) + if verbose;println("J converged @$it");end + JTF[1] = true;break + end + BLAS.gemm!('N','N',-1.0,R,V,1.0,JV) #mul!(JV,R,V,-1.0,1.0) + s_vs = @views Jvs[1:it-1] + bl_ReORTH(q,it,mdim,s_vs,JV,Vt,R) + bl_QR!(JV',Beta_J,mdim,q) + + tnorm = norm(Diagonal(Beta_J),Inf) + add_bl_T!(q,it,Jmat,Beta_J) + if tnorm < Jtol + println("Jbn norm $tnorm") + JTF[1]=true;break + end + end + if rescount == 20;println("JJ not converged");return nothing;end + if JTF[1]==false + tJmat = @views Jmat[1:inow*q,1:inow*q] + bl_ThickRestart(q,Jvs,Jvret,Beta_J,Jmat,inow,1,mdim,Vt,false) + itmin = 2;rescount += 1 + end + end + if inow > itmin + Jcompress(q,Jvs,Jmat,inow,1,mdim,R,bnout,U,Mat,Beta_J) + end + return nothing +end + +function bl_operate_J!(q,Rvec,Jv,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd,beta_J=1.0) + lblock=length(pbits) + @inbounds @threads for bi in Jidxs + if bi==0;continue;end + idim = tdims[bi] + lNn = length(nbits[bi]) + lNp = length(pbits[bi]) + opPP = oPP[bi] + opNN = oNN[bi] + offset = idim-lNn + @inbounds for tmp in opPP + Npi =tmp.Mi; Npf=tmp.Mf; fac=tmp.fac * beta_J + tMi = offset + Npi*lNn + tMf = offset + Npf*lNn + @inbounds for nidx = 1:lNn + Mi = tMi+nidx; Mf = tMf+nidx + w_f1 = @views Jv[:,Mf] + w_i1 = @views Rvec[:,Mi] + w_f2 = @views Jv[:,Mi] + w_i2 = @views Rvec[:,Mf] + @inbounds for b=1:q + w_f1[b] += fac * w_i1[b] + w_f2[b] += fac * w_i2[b] + end + end + end + @inbounds for tmp in opNN #nn + Nni =tmp.Mi; Nnf=tmp.Mf; fac=tmp.fac * beta_J + tMi = offset + Nni + tMf = offset + Nnf + @inbounds for pidx = 1:lNp + Mi = tMi+pidx*lNn; Mf = tMf+pidx*lNn + w_f1 = @views Jv[:,Mf] + w_i1 = @views Rvec[:,Mi] + w_f2 = @views Jv[:,Mi] + w_i2 = @views Rvec[:,Mf] + @inbounds for b=1:q + w_f1[b] += fac * w_i1[b] + w_f2[b] += fac * w_i2[b] + end + end + end + end + @inbounds @threads for bi = 2:lblock + operator = oPNd[bi] + bf = bi-1 + l_Nn_i = length(nbits[bi]) + l_Nn_f = length(nbits[bf]) + off_i = tdims[bi] - l_Nn_i + off_f = tdims[bf] - l_Nn_f + @inbounds for top in operator + pj = top.pjump + nj = top.njump + fac =top.fac * beta_J + @inbounds for tmp_p in pj + phase_p=tmp_p.phase + tMi = off_i + tmp_p.i * l_Nn_i + tMf = off_f + tmp_p.f * l_Nn_f + @inbounds for tmp_n in nj + phase_n=tmp_n.phase + Mi = tMi + tmp_n.i + Mf = tMf + tmp_n.f + coeff = ifelse(phase_p!=phase_n,-fac,fac) + w_f = @views Jv[:,Mf]; w_i = @views Rvec[:,Mi] + @inbounds for b=1:q + w_f[b] += coeff * w_i[b] + end + end + end + end + end + @inbounds @threads for bi = 1:lblock-1 + operator = oPNu[bi] + bf = bi+1 + l_Nn_i = length(nbits[bi]) + l_Nn_f = length(nbits[bf]) + off_i = tdims[bi] - l_Nn_i + off_f = tdims[bf] - l_Nn_f + @inbounds for top in operator + pj = top.pjump + nj = top.njump + fac = top.fac * beta_J + @inbounds for tmp_p in pj + phase_p=tmp_p.phase + tMi = off_i + tmp_p.i * l_Nn_i + tMf = off_f + tmp_p.f * l_Nn_f + @inbounds for tmp_n in nj + Nni = tmp_n.i; Nnf = tmp_n.f; phase_n=tmp_n.phase + Mi = tMi + tmp_n.i + Mf = tMf + tmp_n.f + coeff = ifelse(phase_p!=phase_n,-fac,fac) + w_f = @views Jv[:,Mf]; w_i = @views Rvec[:,Mi] + @inbounds for b=1:q + w_f[b] += coeff * w_i[b] + end + end + end + end + end + return nothing +end diff --git a/test/HFMBPT_IMSRG_test.jl b/test/HFMBPT_IMSRG_test.jl index 49381ea8..7f5bb740 100644 --- a/test/HFMBPT_IMSRG_test.jl +++ b/test/HFMBPT_IMSRG_test.jl @@ -51,7 +51,7 @@ @testset "testing DMD" begin nuc = "He4"; emax = 2; hw = 20 s_pred = [30.0, 50.0] - ds = 0.25; smin = 15.0 + ds = 0.5; smin = 15.0 # generating Omega sntf = "tbme_emn500n4lo_2n3n_srg10.0hw20emax2.snt.bin" pid = getpid() @@ -73,4 +73,4 @@ end end -end \ No newline at end of file +end diff --git a/test/HFMBPT_IMSRG_test.jl~ b/test/HFMBPT_IMSRG_test.jl~ new file mode 100644 index 00000000..49381ea8 --- /dev/null +++ b/test/HFMBPT_IMSRG_test.jl~ @@ -0,0 +1,76 @@ +@testset "HFMBPT & VS-IMSRG calculations" begin + hw = 20; emax=2 + nuc = "He4"; core = "He4"; vspace="p-shell" + nucs = ["He4"] + sntf = "tbme_em500n3lo_barehw20emax2.snt.bin" + ## HF-MBPT from snt/snt.bin + @testset "HFMBPT results under bare EM500,hw20,e2,nmesh50" begin + Eref = [1.493, -5.805, 0.395] + HFobj1 = hf_main(nucs,sntf,hw,emax;return_obj=true,verbose=true) + Es1 = [HFobj1.E0, HFobj1.EMP2, HFobj1.EMP3] + println("Eref $Eref") + println("Es1 $Es1") + @testset "HF" begin + @test (Es1[1] - Eref[1])^2 < 1.e-4 + end + @testset "EMP2" begin + @test (Es1[2] - Eref[2])^2 < 1.e-4 + end + @testset "EMP3" begin + @test (Es1[3] - Eref[3])^2 < 1.e-4 + end + @testset "snt & snt.bin must give identical results" begin + tsntf = replace(sntf,".snt.bin" => ".snt") + HFobj2 = hf_main(nucs,tsntf,hw,emax;return_obj=true,fn_params="parameters/optional_parameters.jl") + Es2 = [HFobj2.E0, HFobj2.EMP2, HFobj2.EMP3] + @test ((HFobj1.E0-HFobj2.E0)^2 + (HFobj1.EMP2-HFobj2.EMP2)^2 + (HFobj1.EMP3-HFobj2.EMP3)^2) < 1.e-6 + end + end + Eref = -4.05225276 + @testset "IMSRG results under bare EM500,hw20,e2,nmesh50" begin + IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,return_obj=true,fn_params="parameters/optional_parameters.jl") + Es = IMSRGobj.H.zerobody[1] + @test abs(Eref-Es[1]) < 1.e-6 + end + @testset "VSIMSRG results under bare EM500,hw20,e2,nmesh50" begin + IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,corenuc=core,ref="nuc",valencespace=vspace,return_obj=true,fn_params="parameters/optional_parameters.jl") + Es = IMSRGobj.H.zerobody[1] + @test abs(Eref - Es[1]) < 1.e-6 + end + + @testset "testing HFMBPT/IMSRG for Rp2 with BetaCM=1.0" begin + nucs = ["He4"] + sntf = "tbme_emn500n4lo_2n3n_srg10.0hw20emax2.snt.bin" + hw = 20 + emax = 2 + IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,Operators=["Rp2"],return_obj=true,fn_params="parameters/optional_parameters.jl") + Rp2_ref = 1.559825^2 + @test abs(IMSRGobj.ExpectationValues["Rp2"] - Rp2_ref)^2 < 1.e-6 + end + + @testset "testing DMD" begin + nuc = "He4"; emax = 2; hw = 20 + s_pred = [30.0, 50.0] + ds = 0.25; smin = 15.0 + # generating Omega + sntf = "tbme_emn500n4lo_2n3n_srg10.0hw20emax2.snt.bin" + pid = getpid() + hf_main([nuc],sntf,hw,emax;doIMSRG=true,Hsample=1,fn_params="parameters/optional_parameters_forDMD.jl") + fn_exact = [ "flowOmega/omega_vec_$(pid)$(nuc)_s$(strip(@sprintf("%6.2f",s))).h5" for s in s_pred] + fns = [ "flowOmega/omega_vec_$(pid)$(nuc)_s$(strip(@sprintf("%6.2f",s))).h5" for s = 15.0:ds:20.0] + # generating DMD vectors + trank = 10 + dmd_main(emax, nuc, fns, trank, smin, ds; s_pred=s_pred, fn_exact=fn_exact) + # restart from DMD + for s in s_pred + s_str = strip(@sprintf("%6.2f",s)) + println("s = $s") + println("DMD:") + fns_conv = ["flowOmega/omega_dmdvec_e$(emax)_$(nuc)_s$(s_str).h5"] + hf_main([nuc],sntf,hw,emax;doIMSRG=true, + restart_from_files=[fns_conv,[]], + fn_params="parameters/optional_parameters_forDMD.jl") + end + end + +end \ No newline at end of file diff --git a/test/ShellModel_test.jl b/test/ShellModel_test.jl index 1aedaf0f..0ca3435b 100644 --- a/test/ShellModel_test.jl +++ b/test/ShellModel_test.jl @@ -57,7 +57,7 @@ end @test prepEC(Hs,target_nuc,num_ev,num_ECsample,targetJ,mode;path_to_samplewav=spath,save_wav=true) == nothing println("make TDmat...") - num_ECsamples = length(Hs) + num_ECsample = length(Hs) mode = "TD" @test prepEC(Hs,target_nuc,num_ev,num_ECsample,targetJ,mode;path_to_samplewav=spath) == nothing diff --git a/test/parameters/optional_parameters_forDMD.jl b/test/parameters/optional_parameters_forDMD.jl index 96b108ac..d5193102 100644 --- a/test/parameters/optional_parameters_forDMD.jl +++ b/test/parameters/optional_parameters_forDMD.jl @@ -10,7 +10,7 @@ chi_order = 3; pottype="em500n3lo" ### --- IMSRG --- smax = 500.0 -dsmax = 0.25 +dsmax = 0.5 denominatorDelta=0.0 BetaCM = 0.0 -magnusmethod="NS" \ No newline at end of file +magnusmethod="NS" diff --git a/test/parameters/optional_parameters_forDMD.jl~ b/test/parameters/optional_parameters_forDMD.jl~ new file mode 100644 index 00000000..96b108ac --- /dev/null +++ b/test/parameters/optional_parameters_forDMD.jl~ @@ -0,0 +1,16 @@ +###--- ChiralEFT --- +n_mesh = 50 +emax = 2 +calc_NN = true +calc_3N = false +hw = 20.0 +srg = false +tbme_fmt = "snt.bin" +chi_order = 3; pottype="em500n3lo" + +### --- IMSRG --- +smax = 500.0 +dsmax = 0.25 +denominatorDelta=0.0 +BetaCM = 0.0 +magnusmethod="NS" \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 274eb6a4..560e8c1b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,9 +4,6 @@ using Test using Printf @testset "NuclearToolkit.jl" begin - println("pwd: ", pwd()) - println("readdir() ",readdir("interaction_file")) - include("chiEFTint_test.jl") include("HFMBPT_IMSRG_test.jl") include("ShellModel_test.jl") From aa29e27e7e87a341563b14c74bd29925afe232d0 Mon Sep 17 00:00:00 2001 From: SotaYoshida Date: Tue, 5 Mar 2024 13:24:17 +0900 Subject: [PATCH 2/4] redundant stdout removed --- src/ShellModel/eigenvector_continuation.jl~ | 1357 ------------------- src/ShellModel/lanczos_methods.jl~ | 871 ------------ test/HFMBPT_IMSRG_test.jl~ | 76 -- test/ShellModel_test.jl | 2 +- 4 files changed, 1 insertion(+), 2305 deletions(-) delete mode 100644 src/ShellModel/eigenvector_continuation.jl~ delete mode 100644 src/ShellModel/lanczos_methods.jl~ delete mode 100644 test/HFMBPT_IMSRG_test.jl~ diff --git a/src/ShellModel/eigenvector_continuation.jl~ b/src/ShellModel/eigenvector_continuation.jl~ deleted file mode 100644 index 8b10ddf2..00000000 --- a/src/ShellModel/eigenvector_continuation.jl~ +++ /dev/null @@ -1,1357 +0,0 @@ -function J2_to_J(J2) - J="" - if Int(J2) % 2 == 0 - J=string(Int(div(J2,2))) - else - J=string(Int(J2))*"/2" - end - return J -end - -function myCholesky!(tmpA,ln::Int64,cLL::LowerTriangular{Float64,Array{Float64,2}}) - l11 = sqrt(tmpA[1,1]) - cLL[1,1] = l11 - cLL[2,1] = tmpA[2,1]/l11; cLL[2,2] = sqrt( tmpA[2,2]-cLL[2,1]^2) - for i=3:ln - for j=1:i-1 - cLL[i,j] = tmpA[i,j] - for k = 1:j-1 - cLL[i,j] += - cLL[i,k]*cLL[j,k] - end - cLL[i,j] = cLL[i,j] / cLL[j,j] - end - cLL[i,i] = tmpA[i,i] - for j=1:i-1 - cLL[i,i] += -cLL[i,j]^2 - end - cLL[i,i] = sqrt(cLL[i,i]) - end - nothing -end - -function read_wfinfos(inpf) - f = open(inpf,"r");tlines = readlines(f);close(f) - lines = rm_comment(tlines) - ret = [] - for line in lines - tl = rm_nan(split(line," ")) - push!(ret,[tl[1],parse(Int64,tl[3]),parse(Int64,tl[5])]) - end - return ret -end - -function read_tdmat(inpf) - io=open(inpf,"r") - Dim=read(io,Int) - ln = read(io,Int) - lTD = read(io,Int) - TDs = [ Float64[] for i=1:ln] - infos = [ Int64[] for i=1:ln] - Nmat = zeros(Float64,Dim,Dim) - for nth = 1:ln - i = read(io,Int); j = read(io,Int); Nij = read(io,Int) - infos[nth] = [i,j,Nij] - tN = read(io,Float64) - Nmat[i,j] = tN - Nmat[j,i] = tN - TDs[nth] = [ read(io,Float64) for k=1:lTD] - end - close(io) - return Dim,ln,lTD,Nmat,infos,TDs -end - -function read_wf_from_info(wfinfos,mdim,vs,wpath) - ## e.g., ./wavsamples_sigma3/Mg24_sample_1.wav nth= 1 TDmatidx 1 - for tmp in wfinfos - inpf,nth,TDidx = tmp - if wpath != "./" && wpath != "" - tn = split(inpf,"/") - inpf = wpath* ifelse(wpath[end]=="/","","/") * tn[end] - end - io = open(inpf,"r") - num_ev = read(io,Int32) - mtot = read(io,Int32) - vals = [read(io,Float64) for i = 1:num_ev] - ojs = [read(io,Int32) for i=1:num_ev] - Css = [ [read(io,Float64) for i=1:mdim] for j=1:num_ev] - vs[TDidx] .= Css[nth] - close(io) - end - return nothing -end - -function readRvecs!(mdim,inpf,tJ,vecs,wfinfos;num=100) - io = open(inpf,"r") - num_ev = read(io,Int32) - mtot = read(io,Int32) - vals = [read(io,Float64) for i = 1:num_ev] - ojs = [read(io,Int32) for i=1:num_ev] - Css = [ [read(io,Float64) for i=1:mdim] for j=1:num_ev] - n = minimum([num,num_ev]) - for i=1:n - j = Int64(ojs[i]) - if j == tJ - push!(vecs,Css[i]) - push!(wfinfos,[inpf,i,length(vecs)]) - end - end - close(io) - return nothing -end - -function prepEC(Hs,target_nuc,num_ev,num_ECsample,tJ,mode; - num_history=3,lm=100,ls=15,tol=1.e-9,q=1, - path_to_samplewav="",calc_moment=true, - save_wav = false,tdmatdir="./tdmat/", - gfactors = [1.0,0.0,5.586,-3.826], - effcharge=[1.5,0.5],debugmode=false,doublelanczos=true) - to = TimerOutput() - sntf = Hs[1] - 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) - hw, bpar = init_ho_by_mass(Anum,1) # mass formula - if 16<= Anum <= 40 - hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell - end - Mtot = 0;if Anum % 2 != 0; Mtot = 1;end - if typeof(tJ) != Int64;println("targetJ must be integer!!");exit();end - Mtot = tJ - eval_jj = 0.5*tJ*(tJ/2+1) - if Anum % 2 != Mtot % 2; println("invalid target J");exit();end - - target_el = replace.(target_nuc, string(Anum)=>"") - Z,N,vp,vn = getZNA(target_el,Anum,cp,cn) - msps_p,msps_n,mz_p,mz_n = construct_msps(p_sps,n_sps) - pbits,nbits,jocc_p,jocc_n,Mps,Mns,tdims = occ(p_sps,msps_p,mz_p,vp, - n_sps,msps_n,mz_n,vn,Mtot) - oPP,oNN,oPNu,oPNd = prep_J(tdims,p_sps,n_sps,msps_p,msps_n, - pbits,nbits) - mdim = tdims[end] - lblock=length(pbits) - mdim_print(target_nuc,Z,N,cp,cn,vp,vn,mdim,tJ) - if mode == "sample" - for (iter,sntf) in enumerate(Hs) - @timeit to "make sample" begin - println("sntf: $sntf") - lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsmsnt(sntf,Anum) - pp_2bjump,nn_2bjump = prep_Hamil_T1(p_sps,n_sps,msps_p,msps_n,pbits,nbits,labels,TBMEs) - bVpn,Vpn,delMs = prep_Hamil_pn(p_sps,n_sps,msps_p,msps_n,labels[3],TBMEs[3]) - l_pbit = length(msps_p);l_nbit = length(msps_n) - bis,bfs,p_NiNfs,n_NiNfs,num_task = prep_pn(lblock,pbits,nbits,Mps,delMs,bVpn,Vpn) - bVpn=nothing - block_tasks = make_distribute(num_task) - - Js = [ 0.5*Mtot*(0.5*Mtot+1) for i = 1:num_ev] - Jtasks = zeros(Int64,lblock) - for i = 1:lblock - Jtasks[i] = length(pbits[i])*length(nbits[i]) - end - Jidxs = make_distribute_J(Jtasks) - - 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) - itmin = 1; elit=1 - vks = [ zeros(Float64,mdim) for i=1:lm] - uks = [ zeros(Float64,mdim) for i=1:ls] - initialize_wf(vks[1],"rand",tJ,mdim) - - @timeit to "Lanczos" elit = TRL(vks,uks,Tmat,itmin, - pbits,nbits,jocc_p,jocc_n,SPEs, - pp_2bjump,nn_2bjump,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=doublelanczos,debugmode=debugmode) - Rvecs = @views uks[1:num_ev] - vals,vecs = eigen(@views Tmat[1:elit*q,1:elit*q]) - @inbounds for (nth,Rvec) in enumerate(Rvecs) - Rvec .= 0.0 - @inbounds for k=1:length(vals) - Rvec .+= vecs[k,nth] .* vks[k] - end - Rvec .*= 1.0/sqrt(dot(Rvec,Rvec)) - end - Jv = zeros(Float64,mdim) - for (nth,Rv) in enumerate(Rvecs) - Jv .= 0.0 - operate_J!(Rv,Jv,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd) - Js[nth] += dot(Rv,Jv) - end - totJs = J_from_JJ1.(Js) - println("J $totJs")#;println("Js $Js") - print("En. ");map(x -> @printf("%24.15f ",x), en[1]);print("\n") - csnt = split(split(sntf,"/")[end],".")[1] - oupf= path_to_samplewav*target_nuc*"_"*csnt*".wav" - if tJ != -1;oupf=path_to_samplewav*target_nuc*"_"*csnt*"_j"*string(tJ)*".wav";end - if save_wav - writeRitzvecs(mdim,Mtot,en[1],totJs,Rvecs,oupf) - end - if calc_moment - @timeit to "moment" tx_mom = eval_moment(Mtot,Rvecs,totJs,p_sps,n_sps, - msps_p,msps_n,tdims, - jocc_p,jocc_n,pbits,nbits,bpar, - gfactors,effcharge) - if tx_mom != "";print(tx_mom);end - end - println("") - end - end - elseif mode =="TD" - try - mkdir(tdmatdir) - catch - nothing - end - #@timeit to "read sample wavs." begin - svecs = Vector{Float64}[ ] - wfinfos = [ ] - for fidx = 0:num_ECsample-1 - inpf = path_to_samplewav*target_nuc*"_tmp_$fidx"*"_j"*string(tJ)*".wav" - readRvecs!(mdim,inpf,tJ,svecs,wfinfos;num=num_ev) - end - - io=open(tdmatdir*"wfinfos_"*target_nuc*"_J"*string(tJ)*".dat","w") - for tmp in wfinfos - println(io,split(tmp[1],"/")[end], " nth= ",tmp[2], " TDmatidx ",tmp[3]) - end - close(io) - @timeit to "misc" begin - Dim = length(svecs) - Hmat = zeros(Float64,Dim,Dim) - Nmat = zeros(Float64,Dim,Dim) - println("Dim $Dim") - end - @timeit to "Nmat calc." begin - @threads for j = 1:Dim - v = svecs[j] - for i = 1:j - tN = dot(svecs[i],v) - Nmat[i,j] = tN; Nmat[j,i]=tN[1] - end - end - end - @timeit to "OBTD" begin - numv = length(svecs) - idxs = Tridx[] - for i = 1:length(svecs) - for j=i:length(svecs) - push!(idxs, Tridx(i,j,idx_from_ij(i,j,numv))) - end - end - OBTDs = [ [ zeros(Float64,length(p_sps)), - zeros(Float64,length(n_sps)) - ] for Nij = 1:length(idxs) ] - calcOBTD(OBTDs,idxs,p_sps,n_sps, - msps_p,msps_n,tdims, - jocc_p,jocc_n,pbits,nbits,svecs) - end - - @timeit to "prepTBTD" begin - opTBTD1,pjumps,njumps,bif_idxs = prepTBTD(tJ,idxs,p_sps,n_sps, - msps_p,msps_n, - pbits,nbits, - labels,TBMEs,svecs, - tdims,Mps) - TBTDs= [zeros(Float64,length(oTBMEs)) for Nij = 1:length(idxs)] - end - @timeit to "calcTBTD" begin - calcTBTD(TBTDs,opTBTD1,pjumps,njumps, - pbits,nbits,tdims,svecs, - idxs,bif_idxs, - olabels,oTBMEs,labels,to) - end - @timeit to "write TDs" begin - oupf=tdmatdir*target_nuc*"_J"*string(tJ)*".dat" - lTD = length(p_sps)+length(n_sps)+length(oTBMEs) - io=open(oupf,"w") - write(io,Dim) - write(io,length(idxs)) - write(io,lTD) - for i = 1:length(svecs) - for j=i:length(svecs) - Nij = idx_from_ij(i,j,length(svecs)) - write(io,i);write(io,j) - write(io,Nij) - write(io,Nmat[i,j]) - TDs = Float64[] - for nth=1:length(OBTDs[Nij][1]) - push!(TDs,OBTDs[Nij][1][nth]) - end - for nth=1:length(OBTDs[Nij][2]) - push!(TDs,OBTDs[Nij][2][nth]) - end - for nth=1:length(TBTDs[Nij]) - push!(TDs,TBTDs[Nij][nth]) - end - for TD in TDs - write(io,TD) - end - end - end - close(io) - end - end - if is_show - show_TimerOutput_results(to) - end - return nothing -end - -""" - 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. -```math -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 - -!!! note - All the effective interactions must be in "the same order" and must be consistent with interaction file from which the transition density matrices were made. - -""" -function solveEC(Hs,target_nuc,tJNs; - write_appwav=false,verbose=false, - calc_moment=true,wpath="./",tdmatpath="./", - is_show=false, - gfactors = [1.0,0.0,5.586,-3.826], - effcharge=[1.5,0.5], - exact_logf="") - - if length(tJNs)==0;println("specify targetJ and # of states");exit();end - to = TimerOutput() - sntf = Hs[1] - 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) - hw, bpar = init_ho_by_mass(Anum,1) # mass formula - if 16<= Anum <= 40 - hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell - end - target_el = replace.(target_nuc, string(Anum)=>"") - Z,N,vp,vn = getZNA(target_el,Anum,cp,cn) - msps_p,msps_n,mz_p,mz_n = construct_msps(p_sps,n_sps) - #println("nuc: $target_el") - exlines = "" - try - f = open(exact_logf,"r");exlines = readlines(f);close(f) - catch - nothing - #println("logfile of exact calc. is missing ") - end - - Dims = Int64[] - sumV=[] - wfinfos=[] - lTD = 0 - - for (jidx,tmp) in enumerate(tJNs) - tJ,num_ev_target = tmp - vals = zeros(Float64,num_ev_target) - Mtot = tJ - pbits,nbits,jocc_p,jocc_n,Mps,Mns,tdims = occ(p_sps,msps_p,mz_p,vp,n_sps,msps_n,mz_n,vn,Mtot) - mdim = tdims[end] - @timeit to "read TDmat" begin - inpf = tdmatpath*target_nuc*"_J"*string(tJ)*".dat" - Dim,ln,lTD,Nmat,infos,TDs = read_tdmat(inpf) - push!(Dims,Dim) - if Dim+1 <= num_ev_target; - println("Error!: sample number must be larger than $num_ev_target") - exit() - end - end - #println("for J=$tJ/2 states Dim.=$Dim") - @timeit to "Cholesky" begin - tMat = zeros(Float64,Dim,Dim) - tildH = zeros(Float64,Dim,Dim) - L = LowerTriangular(zeros(Float64,Dim,Dim)) - try - myCholesky!(Nmat,Dim,L) - catch - u = 1.1 * 1.e-16 - shift = tr(Nmat) * u * ( (Dim+2)/(1.0-((Dim+1)*(Dim+2))*u)) - #println("#epsilon prescription is used!! shift: $shift" ) - Nmat = Nmat + shift .* Matrix(1.0I, Dim,Dim) - myCholesky!(Nmat,Dim,L) - end - Linv = inv(L) - end - Hmat = zeros(Float64,Dim,Dim) - - for sntidx = 1:length(Hs) - sntf = Hs[sntidx] - tMat .= 0.0 - @timeit to "Hmat" begin - lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsmsnt(sntf,Anum) - MEs = [SPEs[1],SPEs[2],oTBMEs] - MEs = [(MEs...)...] - for (idx,TD) in enumerate(TDs) - i,j,Nij = infos[idx] - tmp = dot(TD,MEs) - Hmat[i,j] = tmp; Hmat[j,i] = tmp - end - @timeit to "Gen. Eigen" begin - mul!(tMat,Linv,Hmat) - mul!(tildH,Linv',tMat) - vals,vecs = Arpack.eigs(tildH,nev=num_ev_target,which=:SR,tol=1.e-8, maxiter=300) - vals = real.(vals) - println("vals $vals ln $(length(vals)) num_ev_target $num_ev_target") - if verbose - print_vec("$target_nuc 2J=$tJ En(EC)",vals) - end - end - push!(sumV,[sntf,tJ,vals]) - end - ## write out approx. wav. for only the first Hamiltonian - if write_appwav && sntidx == 1 - wfinfos = read_wfinfos("./tdmat/wfinfos_"*target_nuc*"_J"*string(tJ)*".dat") - @timeit to "read sample wavs." begin - vs = [zeros(Float64,mdim) for i=1:Dim] - read_wf_from_info(wfinfos,mdim,vs,wpath) - end - @timeit to "make .appwav" begin - Rvecs = [ zeros(Float64,mdim) for i=1:length(vals)] - for nth = 1:length(vals) - for k = 1:length(vs) - @. Rvecs[nth] += vecs[k,nth] * vs[k] - end - Rvecs[nth] ./= sqrt(sum(Rvecs[nth].^2)) - end - end - @timeit to "write .appwav" begin - csnt = split(split(sntf, "/")[end],".")[1] - wavname = "./appwavs/"*target_nuc*"_"*csnt*"_j"*string(tJ)*".appwav" - Jlist = [tJ for i=1:length(vals)] - writeRitzvecs(mdim,Mtot,vals,Jlist,Rvecs,wavname) - end - end - end - end - if length(Hs) > 1 && exact_logf != "" - @timeit to "scatter plot" plot_EC_scatter(target_nuc,Hs,sumV,tJNs,Dims,exlines) - end - if is_show - show_TimerOutput_results(to) - end - return nothing -end - - -function solveEC_UQ(Hs,target_nuc,tJNs,Erefs,errors; - itnum_MCMC = 5000,burnin=500,var_M=2.e-2, - calc_moment=true,wpath="./", - is_show=false, - gfactors = [1.0,0.0,5.586,-3.826], - effcharge=[1.5,0.5], - exact_logf="", - intf="", - num_replica=10) - UQ=true; write_appwav=false - if length(tJNs)==0;println("specify targetJ and # of states");exit();end - to = TimerOutput() - - sntf = Hs[1] - Anum = parse(Int64, match(reg,target_nuc).match) - lp,ln,cp,cn,massop,Aref,pow,p_sps,n_sps,SPEs,olabels,oTBMEs,labels,TBMEs = readsnt(sntf,Anum) - hw, bpar = init_ho_by_mass(Anum,1) # mass formula - if 16<= Anum <= 40 - hw, bpar = init_ho_by_mass(Anum,2) # 2: mass formula for sd-shell - end - target_el = replace.(target_nuc, string(Anum)=>"") - Z,N,vp,vn = getZNA(target_el,Anum,cp,cn) - msps_p,msps_n,mz_p,mz_n = construct_msps(p_sps,n_sps) - println("nuc: $target_nuc") - - Dims = Int64[] - lTD = 0 - lnJ = length(tJNs) - HNmats = [ [ zeros(Float64,2,2) ] for i =1:lnJ] - for i=1:lnJ;deleteat!(HNmats[i],1);end - AllTDs = [ [ zeros(Float64,2) ] ]; deleteat!(AllTDs,1) - Allinfos = [ [ [0,0] ] ]; deleteat!(Allinfos,1) - - for (jidx,tmp) in enumerate(tJNs) - tJ,num_ev_target = tmp - Mtot = tJ - @timeit to "read TDmat" begin - inpf = "./tdmat/"*target_nuc*"_J"*string(tJ)*".dat" - Dim,ln,lTD,Nmat,infos,TDs = read_tdmat(inpf) - push!(Dims,Dim) - if Dim+1 <= num_ev_target; - println("Error!: sample number must be larger than $num_ev_target") - exit() - end - end - @timeit to "Cholesky" begin - tMat = zeros(Float64,Dim,Dim) - tildH = zeros(Float64,Dim,Dim) - L = LowerTriangular(zeros(Float64,Dim,Dim)) - try - myCholesky!(Nmat,Dim,L) - catch - u = 1.1 * 1.e-16 - shift = tr(Nmat) * u * ( (Dim+2)/(1.0-((Dim+1)*(Dim+2))*u)) - println("#epsilon prescription is used!! shift: $shift" ) - Nmat = Nmat + shift .* Matrix(1.0I, Dim,Dim) - myCholesky!(Nmat,Dim,L) - end - Linv = inv(L) - end - Hmat = zeros(Float64,Dim,Dim) - push!(AllTDs,copy(TDs)) - push!(Allinfos,copy(infos)) - push!(HNmats[jidx],Hmat) - push!(HNmats[jidx],Linv) - end - - ## MCMC - in_snt =false - Theta = zeros(Float64,length(SPEs[1])+length(SPEs[2])+length(oTBMEs)) - lME = length(Theta) - tDim = Dims[1] - tMat = zeros(Float64,tDim,tDim) - tildH = zeros(Float64,tDim,tDim) - label_T1,label_T0 = make_int(p_sps,n_sps) - idx_s_from_i,facs = int2snt(p_sps,n_sps,label_T1,label_T0,olabels) - - ### random walk in .snt space (not used now) - #MEs = [SPEs[1],SPEs[2],oTBMEs] - #Theta .= [(MEs...)...] - #Thetas = [ zeros(Float64,lME) ]; deleteat!(Thetas,1) - - ### random walk in .int space - ### normal Metropolis-Hastings - #c_Theta = zeros(Float64,lME); c_Theta .= Theta - # iSPEs,V1s,V0s = readVint(intf,label_T1,label_T0) - # ln_int = length(iSPEs)+length(V1s)+length(V0s) - # Vint = zeros(Float64,ln_int) - # iThetas = [ zeros(Float64,ln_int)];deleteat!(iThetas,1) - # nSPEs = copy(iSPEs); nV1s = copy(V1s); nV0s = copy(V0s) - # evalVsnt!(Theta,iSPEs,V1s,V0s,idx_s_from_i,facs) - # tx, llhs, Ens = intMCMC(itnum_MCMC,burnin,var_M,iThetas,Theta,c_Theta,Vint, - # nSPEs,nV1s,nV0s,iSPEs,V1s,V0s, - # idx_s_from_i,facs, - # lnJ,tJNs,HNmats,AllTDs,Allinfos, - # tMat,tildH,Erefs,errors,to) - # plot_MCMC(iThetas,llhs,tJNs,Ens,Erefs) - - ### parallel tempering - oSPEs,oV1s,oV0s = readVint(intf,label_T1,label_T0) - iSPEs = [ copy(oSPEs) for i = 1:num_replica] - V1s = [ copy(oV1s) for i = 1:num_replica] - V0s = [ copy(oV0s) for i = 1:num_replica] - Ts = [ 1.0 + 0.25*(i-1) for i=1:num_replica] - #Ts = [ 1.0 *i for i=1:num_replica] - - ln_int = length(iSPEs[1])+length(V1s[1])+length(V0s[1]) - Vint = zeros(Float64,ln_int) - iThetas = Vector{Float64}[ ] - #AllThetas = [ [ zeros(Float64,ln_int)] for i = 1:num_replica] - #for i=1:num_replica; deleteat!(AllThetas[i],1); end - nSPEs = deepcopy(iSPEs); nV1s = deepcopy(V1s); nV0s = deepcopy(V0s) - evalVsnt!(Theta,iSPEs[1],V1s[1],V0s[1],idx_s_from_i,facs) - Thetas = [ copy(Theta) for i=1:num_replica] - c_Thetas = deepcopy(Thetas) - tx, llhs, Ens = intMCMC_PT(itnum_MCMC,burnin,var_M,Ts, - iThetas,Thetas,c_Thetas,Vint, - nSPEs,nV1s,nV0s,iSPEs,V1s,V0s, - oSPEs,oV1s,oV0s, - idx_s_from_i,facs, - lnJ,tJNs,HNmats,AllTDs,Allinfos, - tMat,tildH,Erefs,errors,to) - - Ns = length(iThetas) - #for n =1:ln_int - # println("i= ", @sprintf("%4i", n), - # "mean: ",@sprintf("%10.4f ", mean([iThetas[i][n] for i =1:Ns])), - # "std: ", @sprintf("%10.4f ", std([iThetas[i][n] for i =1:Ns]))) - #end - println(tx) - write_history(iThetas,Ts,llhs,tJNs,Ens) - if is_show - show_TimerOutput_results(to) - end - println("") - return nothing -end - -function write_history(iThetas,Ts,llhs,tJNs,Ens) - try - mkdir("history") - catch - nothing - end - nr = length(Ts) - Ns = length(llhs[1]) - num_states = sum([tmp[2] for tmp in tJNs]) - - oupf="./history/PTlog_Theta_lowestT.dat" - io = open(oupf,"w") - write(io,nr) - write(io,Ns) - write(io,length(iThetas[1])) - for i=1:Ns - write(io,iThetas[i]) - end - close(io) - oupf="./history/PTlog_llhs.dat" - io = open(oupf,"w") - write(io,nr) - write(io,Ns) - write(io,Ts) - for i=1:nr - write(io,llhs[i]) - end - close(io) - - oupf="./history/PTlog_Ens.dat" - io = open(oupf,"w") - write(io,nr) - write(io,Ns) - for i=1:num_states - write(io,Ens[i]) - end - close(io) -end - -function intMCMC(itnum_MCMC,burnin,var_M,iThetas,Theta,c_Theta,Vint,Vopt, - nSPEs,nV1s,nV0s,SPEs,V1s,V0s, - idx_s_from_i,facs, - lnJ,tJNs,HNmats,AllTDs,Allinfos, - tMat,tildH,Erefs,errors,to) - llh = -1.e+10 - nllh = -1.e+10 - Accept = true - Acchit = 0 - llhs = Float64[ ] - tnum_states = sum([tJNs[i][2] for i =1:lnJ]) - Ens = [ Float64[] for i=1:tnum_states] - offs = [0] - for i = 2:lnJ - push!(offs,sum([tJNs[j][2] for j =1:i-1])) - end - for itM = 1:itnum_MCMC - #println("it=$itM") - @timeit to "int." begin - proposal_Vint!(nSPEs,nV1s,nV0s,SPEs,V1s,V0s,var_M) - evalVsnt!(c_Theta,nSPEs,nV1s,nV0s,idx_s_from_i,facs) - end - nllh = 0.0 - @timeit to "EC" @inbounds for jidx = 1:lnJ - @timeit to "Hmats" begin - num_ev_target = tJNs[jidx][2] - Hmat = HNmats[jidx][1] - Linv = HNmats[jidx][2] - TDs = AllTDs[jidx] - infos = Allinfos[jidx] - @inbounds @threads for idx = 1:length(TDs) - TD = TDs[idx] - i,j,Nij = infos[idx] - tmp = dot(TD,c_Theta) - Hmat[i,j] = tmp; Hmat[j,i] = tmp - end - end - BLAS.gemm!('N','N',1.0,Linv,Hmat,0.0,tMat) - BLAS.gemm!('T','N',1.0,Linv,tMat,0.0,tildH) - - #@timeit to "Arpack" vals = real.(Arpack.eigs(tildH,nev=num_ev_target,which=:SR,tol=1.e-8, maxiter=300))[1] - vals = real.(eigsolve(tildH,num_ev_target,:SR;tol=1.e-8, maxiter=300)[1])[1:num_ev_target] - - nllh += L2_llh(vals,Erefs[jidx],errors[jidx]) - if itM > burnin - @timeit to "push" for i=1:num_ev_target - push!(Ens[offs[jidx]+i],vals[i]) - end - end - end - if nllh > llh - Accept = true - elseif nllh-llh > log(rand()) - Accept = true - else - Accept = false - end - if Accept - llh = nllh - Theta .= c_Theta - SPEs .= nSPEs - V1s .= nV1s - V0s .= nV0s - Acchit += 1 - evalVint!(Vint,SPEs,V1s,V0s) - end - @timeit to "int." if itM > burnin - push!(iThetas,copy(Vint)) - push!(llhs,llh) - end - end - tx = "Accept rate" * @sprintf("%5.1f ",100*(Acchit-burnin)/(itnum_MCMC-burnin)) - return tx, llhs, Ens -end - -function singleH(TD,i,j,Nij,c_Theta, Hmat) - tmp = dot(TD,c_Theta) - Hmat[i,j]=tmp;Hmat[j,i]=tmp - return nothing -end - -function intMCMC_PT(itnum_MCMC,burnin,var_M,Ts, - iThetas,Thetas,c_Thetas,Vint, - RnSPEs,RnV1s,RnV0s,RSPEs,RV1s,RV0s, - oSPEs,oV1s,oV0s, - idx_s_from_i,facs, - lnJ,tJNs,HNmats,AllTDs,Allinfos, - tMat,tildH,Erefs,errors,to) - nr = length(Ts) # num_replica - llhs = [-1.e+30 for i =1:nr] - nllhs = [-1.e+10 for i =1:nr] - Accept = false; TAcchit=0 - Acchits = zeros(Int64,nr) - xi = copy(Thetas[1]) - xj = copy(Thetas[1]) - tnum_states = sum([tJNs[i][2] for i =1:lnJ]) - Ens = [ Float64[] for i=1:tnum_states] - Evals = deepcopy(Erefs) - illhs = [ Float64[] for i=1:nr] # llh history for lowest temperature - offs = [0];for i = 2:lnJ;push!(offs,sum([tJNs[j][2] for j =1:i-1]));end - - # for destructive_eigen my_eigvals - Dim = size(tMat)[1] - v0 = [randn() for i=1:Dim] - - vks = [ zeros(Float64,Dim) for i = 1:150 ] - for i=1:Dim; vks[1][i] = randn(); end - tmp = dot(vks[1],vks[1]) - vks[1] .*= 1.0/sqrt(tmp) - num_history = 2 - en_s = [ [ zeros(Float64,tmp[2]) for i =1:num_history] for tmp in tJNs ] - #### - - jobs = Vector{Int64}[ ] - for jidx = 1:lnJ - info = Allinfos[jidx] - for idx = 1:length(AllTDs[jidx]) - i,j,Nij = info[idx] - push!(jobs,[jidx,idx,i,j,Nij]) - end - end - - for itM = 1:itnum_MCMC - if itM % div(itnum_MCMC,10) ==0; println("it ",itM);end - @inbounds for ridx = 1:nr - nSPEs = RnSPEs[ridx]; nV1s = RnV1s[ridx]; nV0s = RnV0s[ridx] - SPEs = RSPEs[ridx] ; V1s = RV1s[ridx] ; V0s = RV0s[ridx] - Theta = Thetas[ridx]; c_Theta = c_Thetas[ridx] - proposal_Vint!(nSPEs,nV1s,nV0s,SPEs,V1s,V0s,var_M*(ridx^(1/3))) - evalVsnt!(c_Theta,nSPEs,nV1s,nV0s,idx_s_from_i,facs) - nllhs[ridx] = L2norm(nSPEs,nV1s,nV0s,oSPEs,oV1s,oV0s) - @timeit to "Hmats" @inbounds @threads for job in jobs - jidx,idx,mi,mj,Nij = job - Hmat = @views HNmats[jidx][1] - Linv = @views HNmats[jidx][2] - TD = @views AllTDs[jidx][idx] - info = @views Allinfos[jidx][idx] - singleH(TD,mi,mj,Nij,c_Theta,Hmat) - end - @timeit to "gem&diag" for jidx = 1:lnJ - num_ev_target = tJNs[jidx][2] - Hmat = @views HNmats[jidx][1] - Linv = @views HNmats[jidx][2] - Eval = @views Evals[jidx] - BLAS.gemm!('N','N',1.0,Linv,Hmat,0.0,tMat) - BLAS.gemm!('T','N',1.0,Linv,tMat,0.0,tildH) - - @timeit to "eigsolve" begin - #Eval .= Arpack.eigs(tildH,nev=num_ev_target,which=:SR,tol=1.e-6,maxiter=150)[1] - Eval .= real.(eigsolve(tildH,num_ev_target,:SR)[1])[1:num_ev_target] - end - nllhs[ridx] += L2_llh(Eval,Erefs[jidx],errors[jidx];T=Ts[ridx]) - if ridx == 1 - if itM > burnin - for i=1:num_ev_target - push!(Ens[offs[jidx]+i],Eval[i]) - end - end - end - end - Accept = false - if nllhs[ridx] > llhs[ridx] - Accept =true - elseif nllhs[ridx]-llhs[ridx] > log(rand()) - Accept = true - end - if Accept - llhs[ridx] = nllhs[ridx] - Theta .= c_Theta - SPEs .= nSPEs; V1s .= nV1s;V0s .= nV0s - evalVint!(Vint,SPEs,V1s,V0s) - Acchits[ridx] += 1 - end - if itM > burnin - if ridx == 1; push!(iThetas,copy(Vint)); end - push!(illhs[ridx],llhs[ridx]) - continue - end - if itM == burnin;Acchits .= 0;TAcchit=0;end - end - ### exchange - if itM <= burnin || nr == 1; continue ;end - ri = itM % 2 + 1 - for r1 in ri:2:nr-1 - r2 = r1 + 1 - Ti = Ts[r1]; Tj = Ts[r2] - xi .= Thetas[r1]; xj .= Thetas[r2] - llh_i = llhs[r1]; llh_j = llhs[r2] - logp = (Ti-Tj) * (llh_i/Tj - llh_j/Ti) - TAccept = false - if logp > 0.0 - TAccept=true - elseif log(rand()) < logp - TAccept=true - end - if TAccept && itM > burnin - Thetas[r1] .= xj; Thetas[r2] .= xi - illhs[r1][end] = llh_j - illhs[r2][end] = llh_i - TAcchit += 1 - if r1 ==1 - evalVint!(Vint,RSPEs[r2],RV1s[r2],RV1s[r2]) - iThetas[end] .= Vint - end - end - end - end - rates = (100/(itnum_MCMC-burnin)) .* Acchits - tx = "Accept rates " - for i = 1:nr - tx *= " "*@sprintf("%5.1f ", rates[i]) - end - tx *= "\t exchange rate " *@sprintf("%5.1f ", 100*TAcchit/((itnum_MCMC-burnin)*nr)) - return tx, illhs, Ens -end - -function my_eigvals!(A,T,Dim,vks,en,num_ev_target;TF=[false],tol=1.e-6) - T .= 0.0 - for i=1:length(en); en[i] .= 1.e+10;end - for i=1:Dim; vks[1][i] = randn(); end - tmp = 1.0 / sqrt(dot(vks[1],vks[1])) - vks[1] .*= tmp - - for it = 1:length(vks)-1 - v = vks[it]; Hv = vks[it+1] - mul!(Hv,A,v) - talpha = dot(v,Hv) - T[it,it] = talpha - #diagonalize_T!(it,num_ev_target,T,en,2,TF,tol) - #print_vec("HNLanczos @$it ",en[1]) - if TF[1];break;end - axpy!(-talpha,v,Hv) - svks = @views vks[1:it-1] - ReORTH(it,Hv,svks) - tbeta = sqrt(dot(Hv,Hv)) - tmp = 1.0/tbeta - Hv .*= tmp - T[it+1,it] = tbeta; T[it,it+1] = tbeta - end - return nothing -end - -function proposal_Theta!(c_Theta,Theta,ln,var_M) - @inbounds for i=1:ln - c_Theta[i] = Theta[i] + var_M .* randn() - end - return nothing -end -function L2norm(nSPEs,nV1s,nV0s,oSPEs,oV1s,oV0s;lam=1e1) - s = 0.0 - for i =1:length(nSPEs); s += (nSPEs[i]-oSPEs[i])^2; end - for i =1:length(nV1s); s += (nV1s[i]-oV1s[i])^2; end - for i =1:length(nV0s); s += (nV0s[i]-oV0s[i])^2; end - return - 0.5 * s * lam -end -function L2_llh(Eval,Eref,err;T=1.0,sigma_sm=0.126) - s=0.0 - @inbounds for i = 1:length(Eval) - #s -= (Eval[i] - Eref[i])^2 / (sigma_sm^2 + err[i]^2) #i-dependent - #s -= (Eval[i] - Eref[i])^2 # i-independent err=1.0 - s -= (Eval[i] - Eref[i])^2 / (1.0+err[i]^2) # i-independent err=1.0 - end - return 0.5 * s / (T*length(Eval)) -end - -function read_exact(target_nuc,targetJ,sntf,lines) - hit = 0 - for (i,line) in enumerate(lines) - if occursin("$target_nuc",line) - if parse(Int,split(lines[i+1],"targetJ")[end]) == targetJ - hit = 1 - continue - end - end - if hit ==0;continue;end - csntf = split(sntf,"/")[end] - if occursin(csntf,line) - oJs = split(split(split(lines[i+1],"[")[2],"]")[1],",") - oEs = split(lines[i+2])[2:end] - Es=[] - for (idx,J) in enumerate(oJs) - tJ = Int(2*parse(Float64,J)) - if tJ != targetJ; continue;end - push!(Es,parse(Float64,oEs[idx])) - end - return Es - end - end - return [] -end - -function plot_EC_scatter(target_nuc,Hs,sumV,tJNs,Dims,exlines) - js = [ tJNs[i][1] for i=1:length(tJNs)] - xs = [ Float64[] for i=1:length(js)] - ys = [ Float64[] for i=1:length(js)] - minmax=Float64[1.e+10,-1.e+10] - yerrs = Float64[] - for tmp in sumV - sntf,tJ,EnsEC = tmp - Ens = read_exact(target_nuc,tJ,sntf,exlines) - tl = Float64[] - idx = 0 - for k = 1:length(js) - if js[k]==tJ - idx = k;break - end - end - try - Ens[1]; EnsEC[1] - catch - continue - end - for i=1:1 # only the yrast state - #if EnsEC[i] < Ens[i] - # println("variational err ", - # " EC ",EnsEC[i], " Exact ",Ens) - #end - push!(xs[idx],Ens[i]) - push!(ys[idx],EnsEC[i]) - push!(tl,Ens[i]) - push!(tl,EnsEC[i]) - push!(yerrs,(EnsEC[i] - Ens[i])/abs(Ens[i])) - end - if minimum(tl) < minmax[1] - minmax[1] = minimum(tl) - end - if maximum(tl) > minmax[2] - minmax[2] = maximum(tl) - end - end - if length(Hs) > 1; println(" typical error ", mean(yerrs));end - ##scatter - plt=pyimport("matplotlib.pyplot") - cm=pyimport("matplotlib.cm") - patches=pyimport("matplotlib.patches") - fig = plt.figure(figsize=(6,4)) - ax = fig.add_subplot(111); #axB = fig.add_subplot(212) - ax.set_xlabel("Exact (MeV)") - ax.set_ylabel("EC estimate (MeV)") - cols = ["red","blue","darkgreen","orange","purple","magenta","cyan"] - tms = ["o","^",",","d","p","h","8","v","*","1"] - ax.plot([minmax[1]-5,minmax[2]+5],[minmax[1]-5,minmax[2]+5], - linestyle="dotted",color="gray",alpha=0.6) - for jidx = 1:length(js) - nth = 1 #yrast only - tl = nothing - if nth == 1 - tl = "J="*latexstring(string(J2_to_J(js[jidx]))*"_{"*string(nth)*"}") - end - ax.scatter(xs[jidx],ys[jidx],marker=tms[jidx],s=20,lw=0.8, - zorder=150,color=cols[jidx],facecolor="none",label=tl,alpha=0.5) - end - Anum = parse(Int64, match(reg,target_nuc).match) - nuc_latex = latexstring("{}^{$Anum}")*split(target_nuc,string(Anum))[1] - - tl = " $nuc_latex ("*latexstring("N_s")*"="*string(Dims[1])*")" - ax.text(0.05, 0.9, tl,fontsize=15, transform=ax.transAxes) - ax.legend(loc="lower right",fontsize=12) - plt.savefig("./pic/scatter_ECestimates_"*target_nuc*".pdf",bbox_iinches="tight",pad_inches=0) - plt.close() -end - - -function ngauss(x,mu,varV) - return 1.0/(sqrt(2*pi*varV)) * exp( - 0.5* (x-mu)^2 / varV ) -end - -struct Tridx - i::Int - j::Int - Nth::Int -end - -function calcOBTD(OBTDs,idxs,p_sps,n_sps, - msps_p,msps_n, - tdims, - jocc_p,jocc_n, - pbits,nbits, - wfs) - lps=length(p_sps);lns=length(n_sps);lblock = length(nbits) - @inbounds for k=1:length(idxs) - #for k=1:length(idxs) - tmp = idxs[k] - i = tmp.i; j=tmp.j; Nij=tmp.Nth - w_i=wfs[i]; w_j = wfs[j] - pOBTD = OBTDs[Nij][1] - nOBTD = OBTDs[Nij][2] - @inbounds for bi = 1:lblock - idim = tdims[bi] - jocc_p_bi = jocc_p[bi] - jocc_n_bi = jocc_n[bi] - pbit = pbits[bi] - nbit=nbits[bi] - l_Nn = length(nbit) - offset = idim - l_Nn - @inbounds for (pidx,pocc) in enumerate(jocc_p_bi) - tM = offset + pidx*l_Nn - @inbounds for (nidx,nocc) in enumerate(jocc_n_bi) - Mi = tM + nidx - wfprod = w_i[Mi] .* w_j[Mi] - @inbounds for i = 1:lps - pOBTD[i] += wfprod .* pocc[i] - end - @inbounds for i = 1:lns - nOBTD[i] += wfprod .* nocc[i] - end - #pOBTD .+= wfprod .* pocc;nOBTD .+= wfprod .* nocc - #axpy!(wfprod,pocc,pOBTD);axpy!(wfprod,nocc,nOBTD) - end - end - end - end - if false - println("OBTDs: ") - for idx in idxs - i = idx.i; j=idx.j; Nij=idx.Nth - if i == j - println("_[$i] Nij $Nij") - println(" proton: ",OBTDs[Nij][1]) - println("neutron: ",OBTDs[Nij][2]);println("") - end - end - end - return nothing -end - -struct T1ifc - nth::Int64 - bi::Int64 - i::Int64 - f::Int64 - fac::Float64 -end -struct T0ifc - i::Int64 - f::Int64 - fac::Float64 -end - -function prepTBTD(tJ,idxs,p_sps,n_sps,msps_p::Array{Array{Int64,1},1},msps_n::Array{Array{Int64,1},1}, - pbits,nbits,labels,TBMEs,wfs::Array{Array{Float64,1}},tdims,Mps) - lblock=length(pbits) - lp = length(msps_p); ln = length(msps_n) - mstates = [msps_p,msps_n] - sps = [p_sps,n_sps] - loffs = [ 0, length(p_sps)] - TBTD1 = [ T1ifc[ ] , T1ifc[] ] - n = nthreads() - #@threads - for vrank =1:2 #pp:1, nn:2 - loff = loffs[vrank] - T1info = TBTD1[vrank] - vecs= [ [ [ false for i = 1:lp] for j=1:2], - [ [ false for i = 1:ln] for j=1:2]] - @inbounds for (i,ME) in enumerate(TBMEs[vrank]) - a,b,c,d,totJ = labels[vrank][i] - J2 = 2*totJ - ja,ma_s,ma_idxs = possible_mz(sps[vrank][a-loff],mstates[vrank]) - jb,mb_s,mb_idxs = possible_mz(sps[vrank][b-loff],mstates[vrank]) - jc,mc_s,mc_idxs = possible_mz(sps[vrank][c-loff],mstates[vrank]) - jd,md_s,md_idxs = possible_mz(sps[vrank][d-loff],mstates[vrank]) - @inbounds for (ic,mc) in enumerate(mc_s) - @inbounds for (id,md) in enumerate(md_s) - if c == d && mc >= md; continue;end - if abs(mc + md) > J2; continue;end - M_ani = mc + md - initialize_tvec!(vecs[vrank][1]); vecs[vrank][1][mc_idxs[ic]] = true - bit_c = bitarr_to_int(vecs[vrank][1]) - initialize_tvec!(vecs[vrank][1]); vecs[vrank][1][md_idxs[id]] = true - bit_d = bitarr_to_int(vecs[vrank][1]) - @inbounds for (ia,ma) in enumerate(ma_s) - @inbounds for (ib,mb) in enumerate(mb_s) - if a == b && ma>=mb; continue;end - if ma + mb != M_ani;continue;end - initialize_tvec!(vecs[vrank][2]);vecs[vrank][2][ma_idxs[ia]] = true - bit_a = bitarr_to_int(vecs[vrank][2]) - initialize_tvec!(vecs[vrank][2]);vecs[vrank][2][mb_idxs[ib]] = true - bit_b = bitarr_to_int(vecs[vrank][2]) - CG1 = clebschgordan(Float64,ja//2,ma//2,jb//2,mb//2,J2//2,M_ani//2) - CG2 = clebschgordan(Float64,jc//2,mc//2,jd//2,md//2,J2//2,M_ani//2) - bitlist = bit2b(bit_a,bit_b,bit_c,bit_d) - coeff = sqrt( (1.0+deltaf(a,b)) *(1.0+deltaf(c,d)) ) * CG1 * CG2 - if vrank==1 # pp - @inbounds for bi = 1:lblock - TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1] - idim = tdims[bi] - l_Nn = length(nbits[bi]) - @inbounds for (Npi,Phi) in enumerate(pbits[bi]) - TF_connectable(Phi,bitlist,TF) - if TF[1]==false; continue;end - calc_phase!(Phi,bitlist,lp,ret) - bisearch!(pbits[bi],ret[2],ridx) - Npf = ridx[1] - fac = coeff .*ret[1] - push!(T1info,T1ifc(i,bi,Npi,Npf,fac)) - end - end - else ## nn - @inbounds for bi = 1:lblock - TF=[true]; ret=[1,-1]; ridx=[-1,-1,-1] - idim = tdims[bi] - l_Nn = length(nbits[bi]) - @inbounds for (Nni,Phi) in enumerate(nbits[bi]) - TF_connectable(Phi,bitlist,TF) - if TF[1]==false; continue;end - calc_phase!(Phi,bitlist,ln,ret) - bisearch!(nbits[bi],ret[2],ridx) - Nnf = ridx[1] - fac = coeff .*ret[1] - push!(T1info,T1ifc(i,bi,Nni,Nnf,fac)) - end - end - end - end - end - end - end - end - end - - ## pn - vrank=3 - loff = loffs[2] # for neutron - delMs = Int64[] - for j = 1:lp - for i = 1:lp - push!(delMs,msps_p[i][5]-msps_p[j][5]) - end - end - unique!(delMs);sort!(delMs,rev=true) - maxDeltaM = maximum(delMs) - bfs = [ Int64[ ] for i = 1:lblock] - bis = [ Int64[ ] for i = 1:lblock] - for i = 1:lblock - for j = i:lblock - if Mps[j] - Mps[i] in delMs - push!(bfs[i],j); push!(bis[j],i) - end - end - end - - vec_ani_p = [false for i = 1:lp];vec_cre_p = [false for i = 1:lp] - vec_ani_n = [false for i = 1:ln];vec_cre_n = [false for i = 1:ln] - retM = [0,0,0] - bif_idxs= Vector{Int64}[ ] - for i = 1:lblock - for j=i:lblock - push!(bif_idxs,[i,j,idx_from_ij(i,j,lblock)]) - end - end - pjumps = [ [ [ [T0ifc(0,0,0.0)] ] ] for i=1:length(TBMEs[vrank]) ] - njumps = [ [ [ [T0ifc(0,0,0.0)] ] ] for i=1:length(TBMEs[vrank]) ] - for i = 1:length(TBMEs[vrank]) - deleteat!(pjumps[i],1) - deleteat!(njumps[i],1) - end - for (nth,ME) in enumerate(TBMEs[vrank]) - a,b,c,d,totJ,vidx = labels[vrank][nth] - J2 = 2*totJ - ja,ma_s,ma_idxs = possible_mz(p_sps[a],msps_p) - jc,mc_s,mc_idxs = possible_mz(p_sps[c],msps_p) - jb,mb_s,mb_idxs = possible_mz(n_sps[b-loff],msps_n) - jd,md_s,md_idxs = possible_mz(n_sps[d-loff],msps_n) - for (ic,mc) in enumerate(mc_s) - initialize_tvec!(vec_ani_p); vec_ani_p[mc_idxs[ic]] = true - bit_c = bitarr_to_int(vec_ani_p) - for (ia,ma) in enumerate(ma_s) - initialize_tvec!(vec_cre_p); vec_cre_p[ma_idxs[ia]] = true - bit_a = bitarr_to_int(vec_cre_p) - Mp = ma - mc - bisearch!(delMs,Mp,retM);idx = retM[1] - for (id,md) in enumerate(md_s) - if abs(mc + md) > J2; continue;end - initialize_tvec!(vec_ani_n); vec_ani_n[md_idxs[id]] = true - bit_d = bitarr_to_int(vec_ani_n) - for (ib,mb) in enumerate(mb_s) - if mb - md + Mp != 0; continue;end - initialize_tvec!(vec_cre_n); vec_cre_n[mb_idxs[ib]] = true - bit_b = bitarr_to_int(vec_cre_n) - CG1 = clebschgordan(Float64,ja//2,ma//2,jb//2,mb//2, - J2//2,(ma+mb)//2) - CG2 = clebschgordan(Float64,jc//2,mc//2,jd//2,md//2, - J2//2,(mc+md)//2) - fac = CG1*CG2 #* sqrt(tJ+1) /sqrt(J2+1.0) - pjump = [ T0ifc[] for i=1:length(bif_idxs)] - njump = [ T0ifc[] for i=1:length(bif_idxs)] - TF=[true]; ret_p=[0,0];ret_n=[0,0] - ridx=[0,0,0] - - for bi = 1:lblock - idim = tdims[bi] - l_Nn_i = length(nbits[bi]) - for bfidx = 1:length(bfs[bi]) - bf = bfs[bi][bfidx] - deltaM = Mps[bf] - Mps[bi] - tfac = fac * ifelse(bi==bf,0.5,1.0) - if (deltaM in delMs) == false;continue;end - fdim = tdims[bf] - l_Nn_f = length(nbits[bf]) - bif_idx = idx_from_ij(bi,bf,lblock) - for (Npi,pPhi) in enumerate(pbits[bi]) - TF_connectable_1(pPhi,bit_a,bit_c,TF) - if TF[1]==false; continue;end - calc_phase_1!(pPhi,bit_a,bit_c,ret_p) - bisearch!(pbits[bf],ret_p[2],ridx) - if ridx[1] == 0;continue;end - Npf = ridx[1] - push!(pjump[bif_idx], - T0ifc(Npi,Npf,tfac*ret_p[1])) - end - for (Nni,nPhi) in enumerate(nbits[bi]) - TF_connectable_1(nPhi,bit_b,bit_d,TF) - if TF[1]==false; continue;end - calc_phase_1!(nPhi,bit_b,bit_d,ret_n) - bisearch!(nbits[bf],ret_n[2],ridx) - if ridx[1] == 0;continue;end - Nnf = ridx[1] - push!(njump[bif_idx], - T0ifc(Nni,Nnf,1.0*ret_n[1])) - end - # if ma==mb==mc==md==-1 - # if a==1 && b==3 && c==1 && ; - # d==3 && totJ==1 && bi==bf - # println("V($a$b$c$d$totJ) hit ", - # labels[vrank][nth], - # " ms $ma $mb $mc $md ", - # " bi $bi bf $bf bifidx $bif_idx") - # println("pjump ", pjump[bif_idx], - # " njump ", njump[bif_idx]) - # println("") - # end - # end - end - end - push!(pjumps[nth],pjump) - push!(njumps[nth],njump) - end - end - end - end - end - return TBTD1, pjumps,njumps,bif_idxs -end - -function calcTBTD(TBTDs, - opTBTD1,pjumps,njumps, - pbits,nbits,tdims,wfs,idxs,bif_idxs, - olabels,oTBMEs,labels,to) - @inbounds @threads for k in eachindex(idxs) - tmp = idxs[k] - i = tmp.i; j=tmp.j; Nij=tmp.Nth - w_i=wfs[i]; w_j = wfs[j] - coeffs = [ zeros(Float64,length(oTBMEs)) for vrank=1:3] - #@timeit to "pp/nn" begin - ## pp/nn - vrank =1 - for tmp in opTBTD1[vrank] ## pp - nth = tmp.nth; bi = tmp.bi - nbit = nbits[bi] - pbit = pbits[bi] - lN = length(nbit) - vidx = labels[vrank][nth][end] - Npi = tmp.i;Npf = tmp.f;fac = tmp.fac - tMi = tdims[bi]+ (Npi-1)*length(nbit) - tMf = tdims[bi]+ (Npf-1)*length(nbit) - @inbounds for nidx = 1:length(nbit) - Mi = tMi+nidx; Mf = tMf+nidx - coeffs[vrank][vidx] += fac .* (w_i[Mf].*w_j[Mi]) - end - end - vrank =2 - for tmp in opTBTD1[vrank] ## nn - nth = tmp.nth; bi = tmp.bi - nbit = nbits[bi] - pbit = pbits[bi] - lN = length(nbit) - vidx = labels[vrank][nth][end] - Nni = tmp.i;Nnf = tmp.f;fac = tmp.fac - tMi = tdims[bi]+ Nni -lN - tMf = tdims[bi]+ Nnf -lN - @inbounds for pidx = 1:length(pbit) - Mi = tMi + pidx*lN; Mf = tMf + pidx*lN - coeffs[vrank][vidx] += fac .* (w_i[Mf].*w_j[Mi]) - end - end - #@timeit to "pn" begin - ### pn - vrank=3 - @inbounds for (nth,label) in enumerate(labels[vrank]) - vidx = label[end] - pjump = pjumps[nth] - njump = njumps[nth] - @inbounds for opidx in eachindex(pjump) - tmp_pj = pjump[opidx] - tmp_nj = njump[opidx] - @inbounds for bidx in eachindex(bif_idxs) - bi,bf,dummy = bif_idxs[bidx] - idim = tdims[bi]; lNi = length(nbits[bi]) - fdim = tdims[bf]; lNf = length(nbits[bf]) - pjs = tmp_pj[bidx] - njs = tmp_nj[bidx] - @inbounds for pj in pjs - Npi = pj.i;Npf = pj.f; fac_p = pj.fac - tMi = idim + (Npi-1)*lNi - tMf = fdim + (Npf-1)*lNf - @inbounds for nj in tmp_nj[bidx] - Nni = nj.i; Nnf=nj.f; fac_n = nj.fac - Mi = tMi + Nni; Mf = tMf + Nnf - coeffs[vrank][vidx] += fac_p .* fac_n .* (w_i[Mf].*w_j[Mi]) - coeffs[vrank][vidx] += fac_p .* fac_n .* (w_i[Mi].*w_j[Mf]) - end - end - end - end - end - for k = 1:3 - TBTDs[Nij] .+= coeffs[k] - end - end - return nothing -end diff --git a/src/ShellModel/lanczos_methods.jl~ b/src/ShellModel/lanczos_methods.jl~ deleted file mode 100644 index 31b855a7..00000000 --- a/src/ShellModel/lanczos_methods.jl~ +++ /dev/null @@ -1,871 +0,0 @@ -function TRL(vks,uks,Tmat,k,pbits,nbits,jocc_p,jocc_n, - SPEs,pp_2bjump,nn_2bjump,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,checkorth=true,debugmode=true) - - mdim = tdims[end] - TF=[false] - - lnJ = ls - Jvs = [zeros(Float64,mdim) for i=1:lnJ] - Jvret = [zeros(Float64,mdim)] - Jmat = zeros(Float64,lnJ,lnJ) - Jvs[1] .= vks[1] - Jtol = 1.e-7 - JTF = [false];beta_J = 1.0 - - ### - num_restart = 0 - Jvlog_before = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - Jvlog_after = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - Hvs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - HvRs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - ### - - if doubleLanczos - Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to) - vks[1] .= Jvs[1] - vks[1] .*= 1.0 / sqrt(dot(vks[1],vks[1])) - end - - elit = 1 - while TF[1]==false - for it = k:lm-1 - vk =vks[it]; vkp1 =vks[it+1] - @timeit to "operate H" begin - operate_H!(vk, vkp1, pbits,nbits,jocc_p,jocc_n,SPEs, - pp_2bjump,nn_2bjump,tdims,bis,bfs,block_tasks, - p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to) - end - if debugmode - Hvs[it] .= vkp1 - end - talpha = dot(vk,vkp1) - Tmat[it,it] = talpha - 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] - @timeit to "ReORTH" ReORTH(it,vkp1,svks) - if debugmode - HvRs[it] .= vkp1 - end - tbeta = sqrt(dot(vkp1,vkp1)) - vkp1 .*= 1.0/tbeta - #if checkorth; Check_Orthogonality(it,vks,en); end - if debugmode - print_vec("En TRL $it ",en[1]) - println(" beta $tbeta") - end - Tmat[it+1,it] = tbeta; Tmat[it,it+1] = tbeta - - if doubleLanczos - if tbeta < Jtol;TF[1]=true;elit=it;break;end - @timeit to "JJ lanczos" begin - for i = 2:length(Jvs);Jvs[i] .=0.0;end - Jmat .= 0.0;JTF[1] = false - Jvs[1] .= vkp1 - if debugmode; Jvlog_before[it] .= vkp1; end ## - Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to) - vkp1 .= Jvs[1] - if debugmode; Jvlog_after[it] .= vkp1; end ## - if TF[1];elit=it;break;end - end - #println(" beta $tbeta beta_inJ $beta_in_J") - sv = @view vks[1:it] - ReORTH(it,vkp1,sv) - vkp1 .*= 1.0 / sqrt(dot(vkp1,vkp1)) - end - end - if TF[1] == false - @timeit to "Restart" begin - num_restart += 1 - if debugmode; write_wf_hdf5(vks, Hvs, HvRs,Tmat,Jvlog_before, Jvlog_after, "JJ_TR_res$(num_restart).h5", false); end - ThickRestart(vks,uks,Tmat,lm,ls,debugmode) - end - end - k = ls+1 - end - return elit -end - -function write_wf_hdf5(vks,Hvs,HvRs,Tmat,Jv_bef, Jv_aft, fname::String, is_block) - # for debug - io = h5open(fname, "w") - mdim = length(vks[1]) - write(io,"len", length(vks)) - write(io,"mdim", mdim) - for i = 1:length(vks) - write(io, "Hvs_$i", Hvs[i][:]) - write(io,"HvRs_$i", HvRs[i][:]) - end - write(io, "Tmat",Tmat) - for i = 1:length(Jv_bef) - write(io, "Jv_bef_$i", Jv_bef[i]) - end - for i = 1:length(Jv_aft) - write(io, "Jv_aft_$i", Jv_aft[i]) - end - close(io) - return nothing -end - -function Jlanczos(Jvs,Jmat,TF,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj, - Jidxs,oPP,oNN,oPNu,oPNd,beta_J,to) - mdim = tdims[end] - lnJ = length(Jvs) - eljit=1; k = 1; inow=k - beta = 0.0 - while JTF[1] == false - for it = k:lnJ-1 - inow = it - vk = Jvs[it]; vkp1 = Jvs[it+1] - operate_J!(vk,vkp1,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd,beta_J) - axpy!(eval_jj,vk,vkp1)## vkp1 .+= eval_jj .* vk - Jmat[it,it] = dot(vk,vkp1) - teval = eigvals(@views Jmat[1:it,1:it])[1] - #println("JJ it ", @sprintf("%3i",it), " eval ", @sprintf("%12.8f",teval)) - if abs(teval-eval_jj) < Jtol - eljit=it;JTF[1] = true;break - end - - svs = @views Jvs[1:it] - ReORTH(it,vkp1,svs) - - beta = sqrt(dot(vkp1,vkp1)) - println(" betaJ $beta") - #if beta < 1.e-4;eljit=it;TF[1]=true;JTF[1]=true;break;end - vkp1 .*= 1.0/beta - Jmat[it+1,it] = beta; Jmat[it,it+1] = beta - eljit = it - end - if JTF[1]==false - ThickRestart_J(Jvs,Jvret,Jmat,eljit+1,1,eval_jj,Jtol) - k=2 - end - end - if inow > k - ThickRestart_J(Jvs,Jvret,Jmat,eljit+1,1,eval_jj,Jtol) - end - return nothing -end - -function operate_H!(wf,twf,pbits,nbits,jocc_p,jocc_n,SPEs,pp_2bjump,nn_2bjump,tdims,bis,bfs,block_tasks,p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to=nothing) - @inbounds @threads for bi in block_tasks - if bi==0; continue;end #empty job - ret = [0,0,0]; ret2 = [0,0,0] - idim = tdims[bi] - l_Np = length(pbits[bi]) - l_Nn = length(nbits[bi]) - offset = idim -l_Nn - #@timeit to "pn1" begin - for (bfidx,bf) in enumerate(bfs[bi]) - bisearch!(delMs,Mps[bf]-Mps[bi],ret) #!! bf = bfs[bi][j] - Vs = Vpn[ret[1]] - fdim = tdims[bf] - l_Nn_f = length(nbits[bf]) - p_NiNf = p_NiNfs[bi][bfidx] - n_NiNf = n_NiNfs[bi][bfidx] - off_f = fdim-l_Nn_f - for (nth,V) in enumerate(Vs) - Npifs = p_NiNf[nth]; Nnifs = n_NiNf[nth] - for Npif in Npifs - tMi = offset+ Npif.i * l_Nn - tMf = off_f + Npif.f * l_Nn_f - phase_p = Npif.phase - for Nnif in Nnifs - Mi = tMi + Nnif.i; Mf = tMf + Nnif.f - phase_n = Nnif.phase - twf[Mi] += ifelse(phase_p!=phase_n,-V,V) * wf[Mf] - end - end - end - end - #@timeit to "pn2" begin - for j = 1:length(bis[bi])-1 #### tbf = bi !!! - tbi = bis[bi][j] - bisearch!(delMs,Mps[bi]-Mps[tbi],ret) - bisearch_ord!(bfs[tbi],bi,ret2) - fdim=tdims[tbi] - l_Nn_i=length(nbits[tbi]) - p_NiNf = p_NiNfs[tbi][ret2[1]] - n_NiNf = n_NiNfs[tbi][ret2[1]] - off_f = fdim - l_Nn_i - for (nth,V) in enumerate(Vpn[ret[1]]) - Npifs = p_NiNf[nth] - Nnifs = n_NiNf[nth] - for Npif in Npifs - tMi = off_f + Npif.i *l_Nn_i #idim <-> fdim - tMf = offset + Npif.f*l_Nn - phase_p = Npif.phase - for Nnif in Nnifs - Mi = tMi + Nnif.i; Mf = tMf + Nnif.f - phase_n = Nnif.phase - twf[Mf] += ifelse(phase_p!=phase_n,-V,V) * wf[Mi] - end - end - end - end - #@timeit to "pp/nn" begin - ### pp/nn interaction - for (Npi,tinfo) in enumerate(pp_2bjump[bi]) - tMi = offset + Npi*l_Nn - for (jj,tmp) in enumerate(tinfo) - tMf = offset + l_Nn*tmp.f - fac = tmp.coef - for nidx = 1:l_Nn - Mi = tMi + nidx; Mf = tMf + nidx - twf[Mf] += fac * wf[Mi] - twf[Mi] += fac * wf[Mf] - end - end - end - for (Nni,tinfo) in enumerate(nn_2bjump[bi]) - tMi = offset + Nni - for (jj,tmp) in enumerate(tinfo) - tMf = offset + tmp.f - fac = tmp.coef - for pidx = 1:l_Np - Mi = tMi + pidx*l_Nn; Mf = tMf + pidx*l_Nn - twf[Mf] += fac * wf[Mi] - twf[Mi] += fac * wf[Mf] - end - end - end - - #@timeit to "1b" begin - ### one-body operator - for pidx = 1:l_Np - tMi = offset + pidx*l_Nn - for nidx =1:l_Nn - Mi = tMi + nidx - twf[Mi] += (dot(SPEs[1],jocc_p[bi][pidx])+ - dot(SPEs[2],jocc_n[bi][nidx])) * wf[Mi] - end - end - end - return nothing -end - -function diagonalize_T!(k::Int64,num_ev::Int64,Tmat,en::Array{Array{Float64,1}},num_history::Int64, TF::Array{Bool,1}, tol::Float64) - for ith = num_history:-1:2 - en[ith] .= en[ith-1] - end - n = minimum([num_ev,k]) - @views en[1][1:n] .= eigvals(@views Tmat[1:k,1:k] )[1:n] - if all( en[2]-en[1] .< tol) ; TF[1]=true; end - return nothing -end - -function ThickRestart(vks,uks,Tmat,lm,ls,debugmode=false,make_sure=true) - vals,vecs = eigen(@views Tmat[1:lm-1,1:lm-1]) - Tmat_before = ifelse(debugmode,copy(Tmat),[0.0]) - beta = Tmat[lm-1,lm] - Tmat .= 0.0 - @inbounds for k = 1:ls - Tmat[k,k] = vals[k] - end - @inbounds for (k,uk) in enumerate(uks) - tmp = beta .* vecs[lm-1,k] - Tmat[ls+1,k] = tmp; Tmat[k,ls+1] = tmp - uk .= 0.0 - for j=1:lm-1 - axpy!(vecs[j,k],vks[j],uk) - end - end - @inbounds for (k,uk) in enumerate(uks) - if make_sure && k > 1 - svs = @view uks[1:k-1] - ReORTH(k-1,uk,svs) - end - vks[k] .= uk .* (1.0 /sqrt(dot(uk,uk))) - end - vks[ls+1] .= vks[lm] - if make_sure - svs = @view vks[1:ls] - ReORTH(ls,vks[ls+1],svs) - end - for k = ls+2:lm - vks[k] .= 0.0 - end - # if debugmode - # io = h5open("restart.hdf5", "w") - # write(io,"lm", lm) - # write(io,"ls", ls) - # write(io, "beta", beta) - # write(io,"vals", vals) - # write(io,"vecs", vecs) - # for i = 1:ls - # write(io, "uks_$i", uks[i]) - # end - # write(io, "Tmat_before", Tmat_before) - # write(io, "Tmat", Tmat) - # close(io) - # end - return nothing -end - -""" -Thick-Restart Block Lanczos method -""" -function TRBL(q,vks,uks,Tmat,Beta_H,pbits,nbits,jocc_p,jocc_n,SPEs, - pp_2bjump,nn_2bjump,bis,bfs,block_tasks, - p_NiNfs,n_NiNfs,Mps,delMs,Vpn,tdims, - eval_jj,oPP,oNN,oPNu,oPNd,Jidxs, - num_ev,num_history,lm,ls_sub,en,tol,to; - doubleLanczos=false,verbose=false,debugmode=false) - ls = q*ls_sub - mdim = tdims[end] - if doubleLanczos && mdim < 100 - @warn "Block Lanczos with J projection is not stable for small systems.\nUse TRL (is_block=false) instead." - end - - TF=[false] - elit = 1 - inow = 1; itmin = 1; itmax = div(lm,q)-1 - rescount = 0 - Vt = zeros(Float64,q,mdim) - R = zeros(Float64,q,q) - Beta_J = zeros(Float64,q,q) - - lnJ = 20 - Jvs = [zeros(Float64,q,mdim) for i=1:lnJ] - Jvret = [zeros(Float64,q,mdim) ] - Jmat = zeros(Float64,lnJ*q,lnJ*q) - JTF = [false] - Jtol = 1.e-7 - U = zeros(Float64,mdim,lnJ*q) - Mat = zeros(Float64,mdim,q) - - Jvlog_before = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - Jvlog_after = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - Hvs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - HvRs = [ zeros(Float64,mdim) for i = 1:ifelse(debugmode,length(vks),1)] - - if doubleLanczos - @timeit to "JJ lanczos" begin - Jvs[1] .= vks[1] - bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,to,Beta_H,U,Mat) - vks[1] .= Jvs[1];Beta_J .= 0.0;JTF[1]=false - for i=1:lnJ;Jvs[i] .= 0.0;end - V = vks[1] - end - end - num_restart = 0 - while TF[1]==false - for it = itmin:itmax - inow = it - V = vks[it] - HV = vks[it+1] - @timeit to "bl_operateH" begin - bl_operate_H!(q,V,HV,pbits,nbits, - jocc_p,jocc_n,SPEs,pp_2bjump,nn_2bjump, - tdims,bis,bfs,block_tasks, - p_NiNfs,n_NiNfs,Vpn,Mps,delMs,to) - end - if debugmode - Hvs[it] .= HV[1,:] - end - BLAS.gemm!('N','T',1.0,V,HV,0.0,R) - if issymmetric(R) == false - @inbounds for i=1:q;for j=i:q - R[i,j] = R[j,i] - end;end - end - - @views Tmat[q*it-q+1:q*it,q*it-q+1:q*it] .= R - diagonalize_T!(it*q,num_ev,Tmat,en,num_history,TF,tol) - if debugmode - print_vec("En TRBL $it ",en[1]) - end - BLAS.gemm!('N','N',-1.0,R,V,1.0,HV)#mul!(HV,R,V,-1.0,1.0) - s_vks = @views vks[1:it-1] - @timeit to "ReORTH" bl_ReORTH(q,it,mdim,s_vks,HV,Vt,R) - if debugmode - HvRs[it] .= HV[1,:] - end - bl_QR!(HV',Beta_H,mdim,q)#myQR!(HV',Beta_H,mdim,q) - - if doubleLanczos - tnorm = norm(Diagonal(Beta_H),Inf) - if tnorm < Jtol - if verbose; println("Hbn norm $tnorm");end - TF[1]=true;elit=it;break - end - @timeit to "JJ lanczos" begin - for i=2:lnJ; Jvs[i] .= 0.0; end - JTF[1] = false; Jmat .= 0.0;Jvs[1] .= HV - Vt .= 0.0; R .= 0.0 - if debugmode; Jvlog_before[it] .= Jvs[1][1,:]; end - bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj,Jidxs,oPP,oNN,oPNu,oPNd,to,Beta_H,U,Mat) - HV .= Jvs[1]; Beta_J .=0.0 - if debugmode; Jvlog_after[it] .= Jvs[1][1,:]; end - end - end - if debugmode - println(" R $R betaH $Beta_H diff $(R-Beta_H)") - end - add_bl_T!(q,it,Tmat,Beta_H) - if TF[1];elit = it;break;end - end - if TF[1] == false - @timeit to "bl_Restart" begin - num_restart += 1 - if debugmode; write_wf_hdf5(vks,Hvs,HvRs,Tmat,Jvlog_before, Jvlog_after,"JJ_TRBL_res$(num_restart).h5",debugmode); end - bl_ThickRestart(q,vks,uks,Beta_H,Tmat,inow,ls_sub,mdim,Vt,debugmode) - end - itmin = ls_sub + 1 - end - end - return elit -end - -function bl_ReORTH(q,it,mdim,vks,HV,Vt,R) - #LinearAlgebra.jl: mul!(C, A, B, α, β) -> C: ABalpha + Cbeta - #BLAS.gemm!(tA, tB, alpha, A, B,beta,C) C := alpha*A*B + beta*C - # :tA,tB 'N': normal 'T': transpose - @inbounds for l = it-1:-1:1 - Vt .= vks[l] - BLAS.gemm!('N','T', 1.0,HV,Vt,0.0,R) #mul!(R,HV,Vt') - BLAS.gemm!('N','N',-1.0,R,Vt,1.0,HV) #mul!(HV,R,Vt,-1.0,1.0) - end - R .= 0.0 - return nothing -end - -# wf:V, twf:HV -function bl_operate_H!(q,wf,twf,pbits,nbits,jocc_p,jocc_n, - SPEs,pp_2bjump,nn_2bjump,tdims,bis,bfs,block_tasks, - p_NiNfs,n_NiNfs,Vpn, Mps,delMs,to=nothing) - #lblock = length(pbits) - @inbounds @threads for bi in block_tasks - if bi==0; continue;end #empty job - ret = [0,0,0]; ret2 = [0,0,0] - idim = tdims[bi] - l_Np = length(pbits[bi]) - l_Nn = length(nbits[bi]) - offset = idim -l_Nn - #@timeit to "pn1" begin - @inbounds for (bfidx,bf) in enumerate(bfs[bi]) - bisearch!(delMs,Mps[bf]-Mps[bi],ret) #!! bf = bfs[bi][j] - Vs = Vpn[ret[1]] - fdim = tdims[bf] - l_Nn_f = length(nbits[bf]) - p_NiNf = p_NiNfs[bi][bfidx] - n_NiNf = n_NiNfs[bi][bfidx] - off_f = fdim-l_Nn_f - @inbounds for (nth,V) in enumerate(Vs) - Npifs = p_NiNf[nth]; Nnifs = n_NiNf[nth] - @inbounds for Npif in Npifs - tMi = offset+ Npif.i * l_Nn - tMf = off_f + Npif.f * l_Nn_f - phase_p = Npif.phase - @inbounds for Nnif in Nnifs - Mi = tMi + Nnif.i; Mf = tMf + Nnif.f - phase_n = Nnif.phase - coeff = ifelse(phase_p!=phase_n,-V,V) - w_f = @views twf[:,Mi] - w_i = @views wf[:,Mf] - @inbounds for b=1:q - w_f[b] += coeff * w_i[b] - end - end - end - end - end - #end - #@timeit to "pn2" begin - @inbounds for j = 1:length(bis[bi])-1 #### tbf = bi !!! - tbi = bis[bi][j] - bisearch!(delMs,Mps[bi]-Mps[tbi],ret) - bisearch_ord!(bfs[tbi],bi,ret2) - fdim=tdims[tbi] - l_Nn_i=length(nbits[tbi]) - p_NiNf = p_NiNfs[tbi][ret2[1]] - n_NiNf = n_NiNfs[tbi][ret2[1]] - off_f = fdim - l_Nn_i - @inbounds for (nth,V) in enumerate(Vpn[ret[1]]) - Npifs = p_NiNf[nth] - Nnifs = n_NiNf[nth] - @inbounds for Npif in Npifs - tMi = off_f + Npif.i *l_Nn_i #idim <-> fdim - tMf = offset + Npif.f*l_Nn - phase_p = Npif.phase - @inbounds for Nnif in Nnifs - Mi = tMi + Nnif.i; Mf = tMf + Nnif.f - phase_n = Nnif.phase - coeff = ifelse(phase_p!=phase_n,-V,V) - w_f = @views twf[:,Mf] - w_i = @views wf[:,Mi] - @inbounds for b=1:q - w_f[b] += coeff * w_i[b] - end - end - end - end - end - #end - #@timeit to "pp/nn" begin - ### pp/nn interaction - @inbounds for (Npi,tinfo) in enumerate(pp_2bjump[bi]) - tMi = offset + Npi*l_Nn - @inbounds for (jj,tmp) in enumerate(tinfo) - tMf = offset + l_Nn*tmp.f - fac = tmp.coef - @inbounds for nidx = 1:l_Nn - Mi = tMi + nidx; Mf = tMf + nidx - w_f1 = @views twf[:,Mf] - w_i1 = @views wf[:,Mi] - w_f2 = @views twf[:,Mi] - w_i2 = @views wf[:,Mf] - @inbounds for b=1:q - w_f1[b] += fac * w_i1[b] - w_f2[b] += fac * w_i2[b] - end - end - end - end - @inbounds for (Nni,tinfo) in enumerate(nn_2bjump[bi]) - tMi = offset + Nni - @inbounds for (jj,tmp) in enumerate(tinfo) - tMf = offset + tmp.f - fac = tmp.coef - @inbounds for pidx = 1:l_Np - Mi = tMi + pidx*l_Nn; Mf = tMf + pidx*l_Nn - w_f1 = @views twf[:,Mf] - w_i1 = @views wf[:,Mi] - w_f2 = @views twf[:,Mi] - w_i2 = @views wf[:,Mf] - @inbounds for b=1:q - w_f1[b] += fac * w_i1[b] - w_f2[b] += fac * w_i2[b] - end - end - end - end - #end - ### one-body operator - @inbounds for pidx = 1:l_Np - tMi = offset + pidx*l_Nn - @inbounds for nidx =1:l_Nn - Mi = tMi + nidx - coeff = (dot(SPEs[1],jocc_p[bi][pidx])+ - dot(SPEs[2],jocc_n[bi][nidx])) - w_f = @views twf[:,Mi] - w_i = @views wf[:,Mi] - @inbounds for b=1:q - w_f[b] += coeff * w_i[b] - end - end - end - end - return nothing -end - -function Jcompress(q,vks,Tmat,inow,ls_sub,mdim,R,bnout,U,Mat,Beta_J;use_LAPACK=false) - lm = q*inow; ls = q*ls_sub - vals = [0.0] - vecs = [0.0;0.0] - if use_LAPACK # used for only for debug (you need "l_diag.F90" by S.Yoshida) - M = Tmat[1:lm,1:lm] - ccall((:diagonalize_double_,"l_diag.so"),Nothing, - (Ref{Int64},Ref{Float64},Ref{Float64},Ref{Float64},Ref{Int64}), - lm,M,vals,vecs,ls) - else - M = @views Tmat[1:lm,1:lm] - vals,vecs = eigen(M) # eigen(Symmetric(M,:L)) - end - - tv = @views vecs[1:ls,1:ls] - BLAS.gemm!('T','N',1.0,tv,bnout,0.0,R)# mul!(R,tv,bnout) - bnout .= R - - Tmat .= 0.0; U .= 0.0 - for i=1:inow - bv = vks[i] - for b=1:q - j = q*(i-1)+b - u = @views U[:,j] - u .= @views bv[b,:] - end - end - tU = @views U[:,1:lm] - tv = @views vecs[1:lm,1:ls] - BLAS.gemm!('T','T',1.0,tv,tU,0.0,vks[1]) - Beta_J .= @views vecs[1:ls,1:ls] - - return nothing -end - -function like_bl_ThickRestart(vks,uks,Tmat,lm,ls,debugmode) - q = 1 - ls_sub = div(ls,q) - mdim = size(vks[1]) - R = beta = Tmat[lm-1,lm] - vals,vecs = eigen(@views Tmat[1:lm-1,1:lm-1]) - Tmat .= 0.0 - - for k = 1:ls; Tmat[k,k] = vals[k]; end - tv = @views vecs[lm-1,1:ls] - r = R .* tv - - Tmat[ls+1,1:ls] .= r - Tmat[1:ls,ls+1] .= r - - for k = 1:ls - uk = uks[k] #.= 0.0 - uk .= 0.0 - for j=1:lm-1 - vk = vks[j] - fac = vecs[j,k] - uk .+= fac .* vk - end - end - for i = 1:ls - vks[i] .= uks[i] - println("norm vks[$i] $(norm(vks[i]))") - end - vks[ls+1] .= vks[lm] - for i= ls_sub+2:length(vks) - vks[i] .= 0.0 - end - for i=2:ls+1 - v = vks[i] - for j = i-1:-1:1 - fac = dot(vks[j],v) - v .-= fac .* vks[j] - end - w = sqrt(dot(v,v)) - v .*= 1.0/w - println(" => norm vks[$i] $(norm(vks[i]))") - end - return nothing -end - -function bl_ThickRestart(q,vks,uks,R,Tmat,inow,ls_sub,mdim,Vt,debug_mode) - lm = q*inow - ls = q*ls_sub - r = zeros(Float64,q,ls) - vals,vecs = eigen(@views Tmat[1:lm,1:lm]) - Tmat_before = copy(Tmat) - - println("βm1 ",Tmat[lm-1,lm], " ", Tmat[lm-2,lm-1]) - - Tmat .= 0.0 - for k = 1:ls; Tmat[k,k] = vals[k]; end - tv = @views vecs[lm-q+1:lm,1:ls] - mul!(r,R,tv) #BLAS.gemm!('N','N',1.0,R,tv,0.0,r) - - println("R = beta? $R") - println("r in blockRestart $r") - println("tv $tv") - - Tmat[ls+1:ls+q,1:ls] .= r - Tmat[1:ls,ls+1:ls+q] .= r' - for bi = 1:ls_sub - uk = uks[bi] #.= 0.0 - uk .= 0.0 - for b=1:q - k = q*(bi-1) + b - for j=1:inow - vk = vks[j] - for bj = 1:q - v = @views vk[bj,:] - idx = q*(j-1) +bj - fac = vecs[idx,k] - @views u_vv = @views uk[b,:] - axpy!(fac, v, u_vv) - end - end - end - end - for i = 1:ls_sub - vks[i] .= uks[i] - end - vks[ls_sub+1] .= vks[inow+1] - for i= ls_sub+2:length(vks) - vks[i] .= 0.0 - end - for i=1:ls_sub+1 - V = vks[i] - if i>1 - for j = i-1:-1:1 - Vt .= vks[j] - mul!(R,V,Vt') - mul!(V,R,Vt,-1.0,1.0) - end - bl_QR!(V',R,mdim,q) - end - end - # if debug_mode - # io = h5open("restart_blocck.hdf5", "w") - # write(io,"lm", lm) - # write(io,"ls", ls) - # write(io,"vals", vals) - # write(io,"vecs", vecs) - # for i = 1:ls - # write(io, "uks_$i", uks[i][1,:]) - # end - # write(io, "Tmat_before", Tmat_before) - # write(io, "Tmat", Tmat) - # close(io) - # end - return nothing -end - -function bl_JJ_Lanczos(q,Jvs,Jmat,Vt,R,Beta_J,JTF,Jtol,Jvret,pbits,nbits,tdims,eval_jj, - Jidxs,oPP,oNN,oPNu,oPNd,to,bnout,U,Mat;verbose=false) - mdim = tdims[end] - lnJ = length(Jvs) - itmin = 1;itmax = lnJ-1; inow = 0 - rescount = 0 - while JTF[1] == false - for it = itmin:itmax - inow = it - V = Jvs[it]; JV = Jvs[it+1] - bl_operate_J!(q,V,JV,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd) - axpy!(eval_jj,V,JV)##JV .+= eval_jj .* V - BLAS.gemm!('N','T',1.0,V,JV,0.0,R) - @inbounds for j=1:q;for i=1:j - R[j,i] = R[i,j] - end;end - @views Jmat[q*it-q+1:q*it,q*it-q+1:q*it] .= R - tJmat = @views Jmat[1:it*q,1:it*q] - - jeval = eigvals(tJmat)[1:q] - if verbose; print_vec("$it jeval=> ",jeval;long=true);end - if jeval[1] - eval_jj < -Jtol - if verbose;println("neg JJ @$it jeval");end - JTF[1] = true;break - end - if all( abs.(jeval.-eval_jj) .< Jtol) - if verbose;println("J converged @$it");end - JTF[1] = true;break - end - BLAS.gemm!('N','N',-1.0,R,V,1.0,JV) #mul!(JV,R,V,-1.0,1.0) - s_vs = @views Jvs[1:it-1] - bl_ReORTH(q,it,mdim,s_vs,JV,Vt,R) - bl_QR!(JV',Beta_J,mdim,q) - - tnorm = norm(Diagonal(Beta_J),Inf) - add_bl_T!(q,it,Jmat,Beta_J) - if tnorm < Jtol - println("Jbn norm $tnorm") - JTF[1]=true;break - end - end - if rescount == 20;println("JJ not converged");return nothing;end - if JTF[1]==false - tJmat = @views Jmat[1:inow*q,1:inow*q] - bl_ThickRestart(q,Jvs,Jvret,Beta_J,Jmat,inow,1,mdim,Vt,false) - itmin = 2;rescount += 1 - end - end - if inow > itmin - Jcompress(q,Jvs,Jmat,inow,1,mdim,R,bnout,U,Mat,Beta_J) - end - return nothing -end - -function bl_operate_J!(q,Rvec,Jv,pbits,nbits,tdims,Jidxs,oPP,oNN,oPNu,oPNd,beta_J=1.0) - lblock=length(pbits) - @inbounds @threads for bi in Jidxs - if bi==0;continue;end - idim = tdims[bi] - lNn = length(nbits[bi]) - lNp = length(pbits[bi]) - opPP = oPP[bi] - opNN = oNN[bi] - offset = idim-lNn - @inbounds for tmp in opPP - Npi =tmp.Mi; Npf=tmp.Mf; fac=tmp.fac * beta_J - tMi = offset + Npi*lNn - tMf = offset + Npf*lNn - @inbounds for nidx = 1:lNn - Mi = tMi+nidx; Mf = tMf+nidx - w_f1 = @views Jv[:,Mf] - w_i1 = @views Rvec[:,Mi] - w_f2 = @views Jv[:,Mi] - w_i2 = @views Rvec[:,Mf] - @inbounds for b=1:q - w_f1[b] += fac * w_i1[b] - w_f2[b] += fac * w_i2[b] - end - end - end - @inbounds for tmp in opNN #nn - Nni =tmp.Mi; Nnf=tmp.Mf; fac=tmp.fac * beta_J - tMi = offset + Nni - tMf = offset + Nnf - @inbounds for pidx = 1:lNp - Mi = tMi+pidx*lNn; Mf = tMf+pidx*lNn - w_f1 = @views Jv[:,Mf] - w_i1 = @views Rvec[:,Mi] - w_f2 = @views Jv[:,Mi] - w_i2 = @views Rvec[:,Mf] - @inbounds for b=1:q - w_f1[b] += fac * w_i1[b] - w_f2[b] += fac * w_i2[b] - end - end - end - end - @inbounds @threads for bi = 2:lblock - operator = oPNd[bi] - bf = bi-1 - l_Nn_i = length(nbits[bi]) - l_Nn_f = length(nbits[bf]) - off_i = tdims[bi] - l_Nn_i - off_f = tdims[bf] - l_Nn_f - @inbounds for top in operator - pj = top.pjump - nj = top.njump - fac =top.fac * beta_J - @inbounds for tmp_p in pj - phase_p=tmp_p.phase - tMi = off_i + tmp_p.i * l_Nn_i - tMf = off_f + tmp_p.f * l_Nn_f - @inbounds for tmp_n in nj - phase_n=tmp_n.phase - Mi = tMi + tmp_n.i - Mf = tMf + tmp_n.f - coeff = ifelse(phase_p!=phase_n,-fac,fac) - w_f = @views Jv[:,Mf]; w_i = @views Rvec[:,Mi] - @inbounds for b=1:q - w_f[b] += coeff * w_i[b] - end - end - end - end - end - @inbounds @threads for bi = 1:lblock-1 - operator = oPNu[bi] - bf = bi+1 - l_Nn_i = length(nbits[bi]) - l_Nn_f = length(nbits[bf]) - off_i = tdims[bi] - l_Nn_i - off_f = tdims[bf] - l_Nn_f - @inbounds for top in operator - pj = top.pjump - nj = top.njump - fac = top.fac * beta_J - @inbounds for tmp_p in pj - phase_p=tmp_p.phase - tMi = off_i + tmp_p.i * l_Nn_i - tMf = off_f + tmp_p.f * l_Nn_f - @inbounds for tmp_n in nj - Nni = tmp_n.i; Nnf = tmp_n.f; phase_n=tmp_n.phase - Mi = tMi + tmp_n.i - Mf = tMf + tmp_n.f - coeff = ifelse(phase_p!=phase_n,-fac,fac) - w_f = @views Jv[:,Mf]; w_i = @views Rvec[:,Mi] - @inbounds for b=1:q - w_f[b] += coeff * w_i[b] - end - end - end - end - end - return nothing -end diff --git a/test/HFMBPT_IMSRG_test.jl~ b/test/HFMBPT_IMSRG_test.jl~ deleted file mode 100644 index 49381ea8..00000000 --- a/test/HFMBPT_IMSRG_test.jl~ +++ /dev/null @@ -1,76 +0,0 @@ -@testset "HFMBPT & VS-IMSRG calculations" begin - hw = 20; emax=2 - nuc = "He4"; core = "He4"; vspace="p-shell" - nucs = ["He4"] - sntf = "tbme_em500n3lo_barehw20emax2.snt.bin" - ## HF-MBPT from snt/snt.bin - @testset "HFMBPT results under bare EM500,hw20,e2,nmesh50" begin - Eref = [1.493, -5.805, 0.395] - HFobj1 = hf_main(nucs,sntf,hw,emax;return_obj=true,verbose=true) - Es1 = [HFobj1.E0, HFobj1.EMP2, HFobj1.EMP3] - println("Eref $Eref") - println("Es1 $Es1") - @testset "HF" begin - @test (Es1[1] - Eref[1])^2 < 1.e-4 - end - @testset "EMP2" begin - @test (Es1[2] - Eref[2])^2 < 1.e-4 - end - @testset "EMP3" begin - @test (Es1[3] - Eref[3])^2 < 1.e-4 - end - @testset "snt & snt.bin must give identical results" begin - tsntf = replace(sntf,".snt.bin" => ".snt") - HFobj2 = hf_main(nucs,tsntf,hw,emax;return_obj=true,fn_params="parameters/optional_parameters.jl") - Es2 = [HFobj2.E0, HFobj2.EMP2, HFobj2.EMP3] - @test ((HFobj1.E0-HFobj2.E0)^2 + (HFobj1.EMP2-HFobj2.EMP2)^2 + (HFobj1.EMP3-HFobj2.EMP3)^2) < 1.e-6 - end - end - Eref = -4.05225276 - @testset "IMSRG results under bare EM500,hw20,e2,nmesh50" begin - IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,return_obj=true,fn_params="parameters/optional_parameters.jl") - Es = IMSRGobj.H.zerobody[1] - @test abs(Eref-Es[1]) < 1.e-6 - end - @testset "VSIMSRG results under bare EM500,hw20,e2,nmesh50" begin - IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,corenuc=core,ref="nuc",valencespace=vspace,return_obj=true,fn_params="parameters/optional_parameters.jl") - Es = IMSRGobj.H.zerobody[1] - @test abs(Eref - Es[1]) < 1.e-6 - end - - @testset "testing HFMBPT/IMSRG for Rp2 with BetaCM=1.0" begin - nucs = ["He4"] - sntf = "tbme_emn500n4lo_2n3n_srg10.0hw20emax2.snt.bin" - hw = 20 - emax = 2 - IMSRGobj = hf_main(nucs,sntf,hw,emax;doIMSRG=true,Operators=["Rp2"],return_obj=true,fn_params="parameters/optional_parameters.jl") - Rp2_ref = 1.559825^2 - @test abs(IMSRGobj.ExpectationValues["Rp2"] - Rp2_ref)^2 < 1.e-6 - end - - @testset "testing DMD" begin - nuc = "He4"; emax = 2; hw = 20 - s_pred = [30.0, 50.0] - ds = 0.25; smin = 15.0 - # generating Omega - sntf = "tbme_emn500n4lo_2n3n_srg10.0hw20emax2.snt.bin" - pid = getpid() - hf_main([nuc],sntf,hw,emax;doIMSRG=true,Hsample=1,fn_params="parameters/optional_parameters_forDMD.jl") - fn_exact = [ "flowOmega/omega_vec_$(pid)$(nuc)_s$(strip(@sprintf("%6.2f",s))).h5" for s in s_pred] - fns = [ "flowOmega/omega_vec_$(pid)$(nuc)_s$(strip(@sprintf("%6.2f",s))).h5" for s = 15.0:ds:20.0] - # generating DMD vectors - trank = 10 - dmd_main(emax, nuc, fns, trank, smin, ds; s_pred=s_pred, fn_exact=fn_exact) - # restart from DMD - for s in s_pred - s_str = strip(@sprintf("%6.2f",s)) - println("s = $s") - println("DMD:") - fns_conv = ["flowOmega/omega_dmdvec_e$(emax)_$(nuc)_s$(s_str).h5"] - hf_main([nuc],sntf,hw,emax;doIMSRG=true, - restart_from_files=[fns_conv,[]], - fn_params="parameters/optional_parameters_forDMD.jl") - end - end - -end \ No newline at end of file diff --git a/test/ShellModel_test.jl b/test/ShellModel_test.jl index 0ca3435b..02659020 100644 --- a/test/ShellModel_test.jl +++ b/test/ShellModel_test.jl @@ -38,7 +38,7 @@ end @testset "Testing EC on shell model" begin - target_nuc = "O18" + target_nuc = "O20" num_ev = 3 targetJ = 0 sntpath = "interaction_file/random_snts/" From cf81371098664080e1821ff3232f82df5fb15e6c Mon Sep 17 00:00:00 2001 From: SotaYoshida Date: Tue, 5 Mar 2024 13:47:07 +0900 Subject: [PATCH 3/4] v0.4.2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index bb3dd5bb..4fd2e4ac 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "NuclearToolkit" uuid = "89bb3bae-bcec-43ae-87b7-9dd181dc6334" authors = ["SotaYoshida and contributors"] -version = "0.4.1" +version = "0.4.2" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" From 953cffc77c6305f1d544bb137ca8a329b88630f6 Mon Sep 17 00:00:00 2001 From: SotaYoshida Date: Tue, 5 Mar 2024 04:50:37 +0000 Subject: [PATCH 4/4] Update log_sample_script.txt --- example/log_sample_script.txt | 2 + ...hell_coreO16refO16_O16_hw20e4_Delta0.0.snt | 175 ++++++++++++++++++ ..._coreO16refO16_O16_hw20e4_Delta0.0_Rp2.snt | 175 ++++++++++++++++++ 3 files changed, 352 insertions(+) create mode 100644 vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0.snt create mode 100644 vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0_Rp2.snt diff --git a/example/log_sample_script.txt b/example/log_sample_script.txt index 2eed9bf1..f72c683a 100644 --- a/example/log_sample_script.txt +++ b/example/log_sample_script.txt @@ -1,3 +1,4 @@ +option in optional_parameters.jl will be used. --- chiEFTparameters used --- n_mesh = 50 pmax_fm = 5.0 @@ -31,6 +32,7 @@ size of dWS (jmax 9 lmax 40 e2max 8 Nnmax 20): # of sp states 15 # of channels 2bstate 55 #TBME = 34320 E(2H): bare = -2.224578 srg = -2.224578 Diff.1.484e-10 +option in optional_parameters.jl will be used. size of dWS (jmax 9 lmax 40 e2max 8 Nnmax 20): dtri 4.46 MB dcgm0 1.11 MB d6j_int 1.11 MB d6j_lj 0.28 MB d9j_lsj 1.76 MB dictHOB 0.44 MB diff --git a/vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0.snt b/vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0.snt new file mode 100644 index 00000000..f149e463 --- /dev/null +++ b/vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0.snt @@ -0,0 +1,175 @@ +!input interaction => tbme_em500n3lo_srg2.0hw20emax4.snt.bin +!Op:Hamiltonian, zerobody: -156.59729163 + 3 3 8 8 + 1 1 0 1 -1 + 2 0 2 3 -1 + 3 0 2 5 -1 + 4 1 0 1 1 + 5 0 2 3 1 + 6 0 2 5 1 + 6 0 20 + 1 1 -0.785963 + 4 4 -4.880650 + 2 2 7.287169 + 5 5 3.369814 + 3 3 -2.024052 + 6 6 -6.294431 + 158 0 20.00000 + 2 3 2 3 4 -1.652141 + 2 3 3 3 4 1.589057 + 3 3 3 3 4 -0.040063 + 1 6 1 6 3 -2.582206 + 1 6 3 4 3 -2.746191 + 1 6 2 5 3 0.042460 + 1 6 2 6 3 -0.720335 + 1 6 3 5 3 0.918921 + 1 6 3 6 3 -1.582865 + 3 4 3 4 3 -2.622435 + 2 5 3 4 3 -0.004356 + 2 6 3 4 3 -0.928007 + 3 4 3 5 3 0.767474 + 3 4 3 6 3 -1.599253 + 2 5 2 5 3 -3.638236 + 2 5 2 6 3 -1.751884 + 2 5 3 5 3 1.790877 + 2 5 3 6 3 0.562478 + 2 6 2 6 3 -0.721422 + 2 6 3 5 3 0.566882 + 2 6 3 6 3 -1.720089 + 3 5 3 5 3 -0.777046 + 3 5 3 6 3 1.731746 + 3 6 3 6 3 -1.089275 + 1 1 1 1 0 -1.760383 + 1 1 2 2 0 -1.263229 + 1 1 3 3 0 -1.508961 + 2 2 2 2 0 -0.384073 + 2 2 3 3 0 -3.955770 + 3 3 3 3 0 -1.697512 + 1 5 1 5 2 -1.467643 + 1 5 1 6 2 0.655391 + 1 5 2 4 2 -0.518974 + 1 5 3 4 2 -2.604268 + 1 5 2 5 2 -0.197157 + 1 5 2 6 2 -0.999313 + 1 5 3 5 2 -1.814110 + 1 5 3 6 2 -0.859748 + 1 6 1 6 2 -1.225061 + 1 6 2 4 2 2.585762 + 1 6 3 4 2 -0.601848 + 1 6 2 5 2 -0.849113 + 1 6 2 6 2 1.444724 + 1 6 3 5 2 1.112417 + 1 6 3 6 2 -0.730121 + 2 4 2 4 2 -1.464434 + 2 4 3 4 2 -0.650044 + 2 4 2 5 2 0.190557 + 2 4 2 6 2 -1.819322 + 2 4 3 5 2 -1.014451 + 2 4 3 6 2 0.875514 + 3 4 3 4 2 -1.257007 + 2 5 3 4 2 -0.872846 + 2 6 3 4 2 -1.122293 + 3 4 3 5 2 -1.472628 + 3 4 3 6 2 -0.737218 + 2 5 2 5 2 -0.409408 + 2 5 2 6 2 0.796247 + 2 5 3 5 2 -0.793732 + 2 5 3 6 2 -0.962897 + 2 6 2 6 2 -2.940835 + 2 6 3 5 2 -2.535385 + 2 6 3 6 2 0.215100 + 3 5 3 5 2 -2.974266 + 3 5 3 6 2 -0.230675 + 3 6 3 6 2 -1.665831 + 1 4 1 4 1 -3.605260 + 1 4 1 5 1 -0.077959 + 1 4 2 4 1 0.098823 + 1 4 2 5 1 -0.270265 + 1 4 2 6 1 -2.297611 + 1 4 3 5 1 2.337913 + 1 4 3 6 1 -1.007728 + 1 5 1 5 1 -2.876169 + 1 5 2 4 1 2.747981 + 1 5 2 5 1 1.146474 + 1 5 2 6 1 0.958555 + 1 5 3 5 1 -1.163835 + 1 5 3 6 1 -0.350453 + 2 4 2 4 1 -2.877172 + 2 4 2 5 1 -1.153929 + 2 4 2 6 1 -1.127880 + 2 4 3 5 1 0.969235 + 2 4 3 6 1 0.364332 + 2 5 2 5 1 -1.046355 + 2 5 2 6 1 0.151921 + 2 5 3 5 1 -0.121601 + 2 5 3 6 1 2.421661 + 2 6 2 6 1 -3.225087 + 2 6 3 5 1 3.005688 + 2 6 3 6 1 -2.895372 + 3 5 3 5 1 -3.353336 + 3 5 3 6 1 2.902879 + 3 6 3 6 1 -0.984511 + 4 4 4 4 0 -2.286403 + 4 4 5 5 0 -1.355113 + 4 4 6 6 0 -1.605889 + 5 5 5 5 0 -0.886480 + 5 5 6 6 0 -4.154933 + 6 6 6 6 0 -2.304438 + 4 5 4 5 2 -0.885813 + 4 5 4 6 2 -1.940868 + 4 5 5 5 2 -0.224157 + 4 5 5 6 2 0.807250 + 4 5 6 6 2 -1.228542 + 4 6 4 6 2 -1.803225 + 4 6 5 5 2 -1.223378 + 4 6 5 6 2 0.308123 + 4 6 6 6 2 -1.015357 + 5 5 5 5 2 -0.359773 + 5 5 5 6 2 1.129547 + 5 5 6 6 2 -0.934619 + 5 6 5 6 2 -0.395232 + 5 6 6 6 2 0.277493 + 6 6 6 6 2 -1.652928 + 5 6 5 6 4 -2.171276 + 5 6 6 6 4 1.676396 + 6 6 6 6 4 -0.521973 + 4 6 4 6 3 0.127249 + 4 6 5 6 3 0.162911 + 5 6 5 6 3 -0.206016 + 1 2 1 2 2 -0.449692 + 1 2 1 3 2 -1.829618 + 1 2 2 2 2 -0.190902 + 1 2 2 3 2 0.755254 + 1 2 3 3 2 -1.157465 + 1 3 1 3 2 -1.281861 + 1 3 2 2 2 -1.152612 + 1 3 2 3 2 0.285414 + 1 3 3 3 2 -0.949715 + 2 2 2 2 2 0.031828 + 2 2 2 3 2 1.065174 + 2 2 3 3 2 -0.888374 + 2 3 2 3 2 0.054567 + 2 3 3 3 2 0.290712 + 3 3 3 3 2 -1.100273 + 2 6 2 6 4 -4.200387 + 2 6 3 5 4 -1.991830 + 2 6 3 6 4 1.179002 + 3 5 3 5 4 -4.254746 + 3 5 3 6 4 -1.202928 + 3 6 3 6 4 -0.531971 + 1 3 1 3 3 0.548646 + 1 3 2 3 3 0.143468 + 2 3 2 3 3 0.183872 + 1 4 1 4 0 -2.339914 + 1 4 2 5 0 -1.399456 + 1 4 3 6 0 -1.642406 + 2 5 2 5 0 -1.050066 + 2 5 3 6 0 -4.202401 + 3 6 3 6 0 -2.483902 + 1 2 1 2 1 0.213977 + 1 2 2 3 1 -0.206459 + 2 3 2 3 1 0.082890 + 3 6 3 6 5 -5.140433 + 4 5 4 5 1 -0.156485 + 4 5 5 6 1 -0.202198 + 5 6 5 6 1 -0.347085 diff --git a/vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0_Rp2.snt b/vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0_Rp2.snt new file mode 100644 index 00000000..df8914ae --- /dev/null +++ b/vsimsrg_sd-shell_coreO16refO16_O16_hw20e4_Delta0.0_Rp2.snt @@ -0,0 +1,175 @@ +!input interaction => tbme_em500n3lo_srg2.0hw20emax4.snt.bin +!Op:Rp2, zerobody: 4.40732392 + 3 3 8 8 + 1 1 0 1 -1 + 2 0 2 3 -1 + 3 0 2 5 -1 + 4 1 0 1 1 + 5 0 2 3 1 + 6 0 2 5 1 + 6 0 20 + 1 1 1.074234 + 4 4 0.034391 + 2 2 1.084849 + 5 5 0.128713 + 3 3 0.930611 + 6 6 -0.017233 + 158 0 20.00000 + 2 3 2 3 4 -0.012151 + 2 3 3 3 4 0.011734 + 3 3 3 3 4 -0.013068 + 1 6 1 6 3 -0.022129 + 1 6 3 4 3 -0.007805 + 1 6 2 5 3 -0.013907 + 1 6 2 6 3 -0.004977 + 1 6 3 5 3 0.002220 + 1 6 3 6 3 -0.003843 + 3 4 3 4 3 -0.007162 + 2 5 3 4 3 0.018294 + 2 6 3 4 3 0.002798 + 3 4 3 5 3 -0.013780 + 3 4 3 6 3 0.002692 + 2 5 2 5 3 -0.010598 + 2 5 2 6 3 -0.018863 + 2 5 3 5 3 -0.002932 + 2 5 3 6 3 0.000929 + 2 6 2 6 3 -0.011713 + 2 6 3 5 3 -0.002660 + 2 6 3 6 3 -0.012623 + 3 5 3 5 3 0.005874 + 3 5 3 6 3 0.005710 + 3 6 3 6 3 -0.004765 + 1 1 1 1 0 -0.021291 + 1 1 2 2 0 0.006633 + 1 1 3 3 0 0.006258 + 2 2 2 2 0 0.023407 + 2 2 3 3 0 -0.022942 + 3 3 3 3 0 0.021405 + 1 5 1 5 2 -0.009117 + 1 5 1 6 2 0.006712 + 1 5 2 4 2 -0.004740 + 1 5 3 4 2 -0.011588 + 1 5 2 5 2 0.003576 + 1 5 2 6 2 -0.001932 + 1 5 3 5 2 -0.004848 + 1 5 3 6 2 -0.002810 + 1 6 1 6 2 -0.012561 + 1 6 2 4 2 0.023620 + 1 6 3 4 2 -0.003574 + 1 6 2 5 2 -0.005889 + 1 6 2 6 2 0.007638 + 1 6 3 5 2 0.007153 + 1 6 3 6 2 -0.001670 + 2 4 2 4 2 -0.004249 + 2 4 3 4 2 -0.009646 + 2 4 2 5 2 -0.003421 + 2 4 2 6 2 -0.003379 + 2 4 3 5 2 0.002095 + 2 4 3 6 2 -0.004630 + 3 4 3 4 2 -0.000643 + 2 5 3 4 2 0.002356 + 2 6 3 4 2 -0.007076 + 3 4 3 5 2 0.005166 + 3 4 3 6 2 0.000968 + 2 5 2 5 2 -0.000240 + 2 5 2 6 2 0.004727 + 2 5 3 5 2 -0.003786 + 2 5 3 6 2 -0.003002 + 2 6 2 6 2 -0.023023 + 2 6 3 5 2 -0.014308 + 2 6 3 6 2 0.004422 + 3 5 3 5 2 -0.009437 + 3 5 3 6 2 0.003708 + 3 6 3 6 2 -0.013830 + 1 4 1 4 1 -0.021971 + 1 4 1 5 1 -0.013245 + 1 4 2 4 1 0.006064 + 1 4 2 5 1 -0.007666 + 1 4 2 6 1 -0.007861 + 1 4 3 5 1 -0.007499 + 1 4 3 6 1 0.004544 + 1 5 1 5 1 -0.007624 + 1 5 2 4 1 0.001600 + 1 5 2 5 1 -0.001485 + 1 5 2 6 1 -0.002671 + 1 5 3 5 1 0.008824 + 1 5 3 6 1 -0.000025 + 2 4 2 4 1 -0.003157 + 2 4 2 5 1 0.006159 + 2 4 2 6 1 -0.004579 + 2 4 3 5 1 -0.003329 + 2 4 3 6 1 -0.003181 + 2 5 2 5 1 0.001831 + 2 5 2 6 1 -0.012446 + 2 5 3 5 1 -0.009776 + 2 5 3 6 1 0.009136 + 2 6 2 6 1 -0.023254 + 2 6 3 5 1 0.000466 + 2 6 3 6 1 -0.012572 + 3 5 3 5 1 0.029905 + 3 5 3 6 1 0.001043 + 3 6 3 6 1 0.011394 + 4 4 4 4 0 -0.000772 + 4 4 5 5 0 0.002011 + 4 4 6 6 0 0.004767 + 5 5 5 5 0 0.004237 + 5 5 6 6 0 0.002668 + 6 6 6 6 0 0.006513 + 4 5 4 5 2 -0.000431 + 4 5 4 6 2 -0.000157 + 4 5 5 5 2 0.000865 + 4 5 5 6 2 0.000155 + 4 5 6 6 2 0.002223 + 4 6 4 6 2 0.000202 + 4 6 5 5 2 0.000556 + 4 6 5 6 2 -0.001496 + 4 6 6 6 2 0.002202 + 5 5 5 5 2 -0.000246 + 5 5 5 6 2 -0.000229 + 5 5 6 6 2 0.000469 + 5 6 5 6 2 0.000090 + 5 6 6 6 2 0.000988 + 6 6 6 6 2 -0.000461 + 5 6 5 6 4 -0.002305 + 5 6 6 6 4 -0.000532 + 6 6 6 6 4 -0.001709 + 4 6 4 6 3 -0.002161 + 4 6 5 6 3 0.000029 + 5 6 5 6 3 -0.001836 + 1 2 1 2 2 -0.000245 + 1 2 1 3 2 -0.016418 + 1 2 2 2 2 0.008551 + 1 2 2 3 2 0.007554 + 1 2 3 3 2 0.001122 + 1 3 1 3 2 -0.014895 + 1 3 2 2 2 -0.004942 + 1 3 2 3 2 -0.010295 + 1 3 3 3 2 -0.002576 + 2 2 2 2 2 0.001635 + 2 2 2 3 2 0.010933 + 2 2 3 3 2 -0.004905 + 2 3 2 3 2 -0.000857 + 2 3 3 3 2 -0.000778 + 3 3 3 3 2 -0.021133 + 2 6 2 6 4 -0.034126 + 2 6 3 5 4 -0.013050 + 2 6 3 6 4 0.011511 + 3 5 3 5 4 -0.011986 + 3 5 3 6 4 0.002336 + 3 6 3 6 4 -0.009781 + 1 3 1 3 3 -0.008996 + 1 3 2 3 3 0.007152 + 2 3 2 3 3 -0.006881 + 1 4 1 4 0 -0.015198 + 1 4 2 5 0 0.004509 + 1 4 3 6 0 0.005077 + 2 5 2 5 0 0.013756 + 2 5 3 6 0 -0.013010 + 3 6 3 6 0 0.010568 + 1 2 1 2 1 -0.004454 + 1 2 2 3 1 0.000099 + 2 3 2 3 1 0.007595 + 3 6 3 6 5 -0.025277 + 4 5 4 5 1 -0.001474 + 4 5 5 6 1 -0.000763 + 5 6 5 6 1 0.001204