diff --git a/src/tucker_outer.jl b/src/tucker_outer.jl index 78ee87a..e1cd316 100644 --- a/src/tucker_outer.jl +++ b/src/tucker_outer.jl @@ -428,8 +428,8 @@ function tucker_cepa_solve(ref_vector::BSTstate{T,N,R}, cepa_vector::BSTstate, c Ecepa = (e0[1] + ecorr[1])/(1+SxC[1]) #@printf(" %s %18.12f\n",cepa_shift, (e0 + ecorr)/(1+SxC)) @printf("Iter: %4d %18.12f %18.12f \n",it,Ec ,Ecepa-e0) - if abs(Ec - (Ecepa-e0)) < 1e-6 - @printf(" Converged %s %18.12f\n",cepa_shift, (e0 + ecorr)/(1+SxC)) + if abs(Ec - (Ecepa-e0)) < 1e-6 + @printf(" Converged %s %18.12f\n",cepa_shift, (e0[1] + ecorr[1])/(1+SxC[1])) break end Ec = Ecepa - e0 @@ -1213,13 +1213,15 @@ function _add_results!(sig_cts, data, compress_twice, thresh, N) end -function do_fois_ci(ref::BSTstate, cluster_ops, clustered_ham; +function do_fois_ci(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; H0 = "Hcmf", max_iter = 50, nbody = 4, thresh_foi = 1e-6, + thresh_pt = 1e-5, tol = 1e-5, - verbose = true) + pt = false, + verbose = true) where {T,N,R} @printf("\n-------------------------------------------------------\n") @printf(" Do CI in FOIS\n") @printf(" H0 = %-s\n", H0) @@ -1240,33 +1242,74 @@ function do_fois_ci(ref::BSTstate, cluster_ops, clustered_ham; println() println(" Compute FOIS. Reference space dim = ", length(ref_vec)) @time pt1_vec = build_compressed_1st_order_state(ref_vec, cluster_ops, clustered_ham, nbody=nbody, thresh=thresh_foi) - - @printf(" Nick: %12.8f\n", sqrt(orth_dot(pt1_vec,pt1_vec))) + for i in 1:R + @printf(" Nick: %12.8f\n", sqrt.(orth_dot(pt1_vec,pt1_vec))[i]) + end project_out!(pt1_vec, ref) # # Compress FOIS - norm1 = sqrt(orth_dot(pt1_vec, pt1_vec)) + norm1 = sqrt.(orth_dot(pt1_vec, pt1_vec)) dim1 = length(pt1_vec) pt1_vec = compress(pt1_vec, thresh=thresh_foi) - norm2 = sqrt(orth_dot(pt1_vec, pt1_vec)) + norm2 = sqrt.(orth_dot(pt1_vec, pt1_vec)) dim2 = length(pt1_vec) @printf(" %-50s%10i → %-10i (thresh = %8.1e)\n", "FOIS Compressed from: ", dim1, dim2, thresh_foi) - @printf(" %-50s%10.2e → %-10.2e (thresh = %8.1e)\n", "Norm of |1>: ",norm1, norm2, thresh_foi) - @printf(" %-50s%10.6f\n", "Overlap between <1|0>: ", nonorth_dot(pt1_vec, ref_vec, verbose=0)) + for i in 1:R + @printf(" %-50s%10.2e → %-10.2e (thresh = %8.1e)\n", "Norm of |1>: ",norm1[i], norm2[i], thresh_foi) + end + for i in 1:R + @printf(" %-50s%10.6f\n", "Overlap between <1|0>: ", nonorth_dot(pt1_vec, ref_vec, verbose=0)[i]) + end nonorth_add!(ref_vec, pt1_vec) # # Solve for first order wavefunction println(" Compute CI energy in the space = ", length(ref_vec)) - pt1_vec, e_pt2= hylleraas_compressed_mp2(pt1_vec, ref_vec, cluster_ops, clustered_ham; tol=tol, max_iter=max_iter, H0=H0) - eci, ref_vec = ci_solve(ref_vec, cluster_ops, clustered_ham, conv_thres=tol) - @printf(" E(Ref) = %12.8f\n", e0[1]) - @printf(" E(CI) tot = %12.8f\n", eci[1]) - return eci[1], ref_vec + if pt==true + pt1_vec, e_pt2= hylleraas_compressed_mp2(pt1_vec, ref_vec, cluster_ops, clustered_ham; tol=tol, max_iter=max_iter, H0=H0) + nonorth_add!(ref_vec, pt1_vec) + ref_vec = compress(ref_vec, thresh=thresh_pt) + end + eci, ref_vec = ci_solve(ref_vec, cluster_ops, clustered_ham;) + for i in 1:R + @printf(" E(Ref) for %ith state = %12.8f\n",i, e0[i]) + @printf(" E(CI) tot for %ith state = %12.8f\n",i, eci[i]) + end + return eci, ref_vec end - - + +""" + do_fois_cepa(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; + max_iter=20, + cepa_shift="cepa", + cepa_mit=30, + nbody=4, + thresh_foi=1e-6, + tol=1e-5, + compress_type="matvec", + prescreen=false, + verbose=true) where {T,N,R} + +Perform Coupled Electron Pair Approximation (CEPA) calculations. + +# Arguments + - `ref::BSTstate{T,N,R}`: The reference `BSTstate` object representing the wavefunction. + - `cluster_ops`: Operators related to the cluster. + - `clustered_ham`: Clustered Hamiltonian. + - `max_iter::Int`: Maximum number of iterations for the CEPA solver. Default is 20. + - `cepa_shift::String`: Type of CEPA calculation. Default is "cepa". + - `cepa_mit::Int`: Maximum number of iterations for the CEPA method. Default is 30. + - `nbody::Int`: Number of bodies for the first-order interaction space. Default is 4. + - `thresh_foi::Float64`: Threshold for the first-order interaction space. Default is 1e-6. + - `tol::Float64`: Convergence threshold. Default is 1e-5. + - `compress_type::String`: Type of compression for the first-order interaction space. Default is "matvec". + - `prescreen::Bool`: Whether to prescreen interactions. Default is false. + - `verbose::Bool`: Whether to print verbose output. Default is true. + +# Returns + - `Vector{Float64}`: A vector containing the CEPA correlation energies. +""" function do_fois_cepa(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; @@ -1310,14 +1353,14 @@ function do_fois_cepa(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; pt1_vec, e_pt2 = hylleraas_compressed_mp2(pt1_vec, ref_vec, cluster_ops, clustered_ham; tol=tol, do_pt=true) end - display(pt1_vec) + # display(pt1_vec) # # Compress FOIS - norm1 = orth_dot(pt1_vec, pt1_vec) + # norm1 = orth_dot(pt1_vec, pt1_vec) dim1 = length(pt1_vec) pt1_vec = compress(pt1_vec, thresh=thresh_foi) - norm2 = orth_dot(pt1_vec, pt1_vec) + # norm2 = orth_dot(pt1_vec, pt1_vec) dim2 = length(pt1_vec) @printf(" %-50s%10i → %-10i (thresh = %8.1e)\n", "FOIS Compressed from: ", dim1, dim2, thresh_foi) #@printf(" %-50s%10.2e → %-10.2e (thresh = %8.1e)\n", "Norm of |1>: ",norm1, norm2, thresh_foi) @@ -1326,6 +1369,119 @@ function do_fois_cepa(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; [@printf("%10.6f", ovlp[r]) for r in 1:R] println() + # + # Solve CEPA + println() + cepa_vec = deepcopy(pt1_vec) + # display(cepa_vec) + # display(ref_vec) + e_cepa_vec=[] + println(" Do CEPA: Dim = ", length(cepa_vec)) + for i in 1:R + ref_vec_i=FermiCG.BSTstate(ref_vec,i) + display(ref_vec_i) + cepa_vec_i=FermiCG.BSTstate(cepa_vec,i) + zero!(cepa_vec_i) + println(" Do CEPA: Dim = ", length(cepa_vec_i)) + @time e_cepa, x_cepa = tucker_cepa_solve(ref_vec_i, cepa_vec_i, cluster_ops, clustered_ham, cepa_shift, cepa_mit, tol=tol, max_iter=max_iter, verbose=verbose) + + @printf(" E(cepa) corr = %12.8f\n", e_cepa[1]) + @printf(" X(cepa) norm = %12.8f\n", sqrt(orth_dot(x_cepa, x_cepa)[1])) + # nonorth_add!(x_cepa, ref_vec_i) + # orthonormalize!(x_cepa) + push!(e_cepa_vec, e_cepa[1]) + end + return e_cepa_vec +end + +""" + do_fois_cepa(ref::BSTstate{T,N,1}, cluster_ops, clustered_ham; + max_iter=20, + cepa_shift="cepa", + cepa_mit=30, + nbody=4, + thresh_foi=1e-6, + tol=1e-5, + compress_type="matvec", + prescreen=false, + verbose=true) where {T,N} + +Perform Coupled Electron Pair Approximation (CEPA) calculations. + +# Arguments + - `ref::BSTstate{T,N,1}`: The reference `BSTstate` object representing the wavefunction with only one root. + - `cluster_ops`: Operators related to the cluster. + - `clustered_ham`: Clustered Hamiltonian. + - `max_iter::Int`: Maximum number of iterations for the CEPA solver. Default is 20. + - `cepa_shift::String`: Type of CEPA calculation. Default is "cepa". + - `cepa_mit::Int`: Maximum number of iterations for the CEPA method. Default is 30. + - `nbody::Int`: Number of bodies for the first-order interaction space. Default is 4. + - `thresh_foi::Float64`: Threshold for the first-order interaction space. Default is 1e-6. + - `tol::Float64`: Convergence threshold. Default is 1e-5. + - `compress_type::String`: Type of compression for the first-order interaction space. Default is "matvec". + - `prescreen::Bool`: Whether to prescreen interactions. Default is false. + - `verbose::Bool`: Whether to print verbose output. Default is true. + +# Returns + - `Tuple{Float64, BSTstate{T,N,1}}`: The CEPA correlation energy and the updated wavefunction. +""" + +function do_fois_cepa(ref::BSTstate{T,N,1}, cluster_ops, clustered_ham; + max_iter=20, + cepa_shift="cepa", + cepa_mit=30, + nbody=4, + thresh_foi=1e-6, + tol=1e-5, + compress_type="matvec", + prescreen=false, + verbose=true) where {T,N} + @printf("\n-------------------------------------------------------\n") + @printf(" Do CEPA\n") + @printf(" thresh_foi = %-8.1e\n", thresh_foi) + @printf(" nbody = %-i\n", nbody) + @printf("\n") + @printf(" Length of Reference = %-i\n", length(ref)) + @printf(" Calculation type = %s\n", cepa_shift) + @printf(" Compression type = %s\n", compress_type) + @printf("\n-------------------------------------------------------\n") + + # + # Solve variationally in reference space + println() + ref_vec = deepcopy(ref) + @printf(" Solve zeroth-order problem. Dimension = %10i\n", length(ref_vec)) + @time e0, ref_vec = ci_solve(ref_vec, cluster_ops, clustered_ham, conv_thresh=tol) + + # + # Get First order wavefunction + println() + println(" Compute FOIS. Reference space dim = ", length(ref_vec)) + @time pt1_vec = build_compressed_1st_order_state(ref_vec, cluster_ops, clustered_ham, nbody=nbody, thresh=thresh_foi, prescreen=prescreen) + + project_out!(pt1_vec, ref) + + if compress_type == "pt_vec" + println() + println(" Compute PT vector. Reference space dim = ", length(ref_vec)) + pt1_vec, e_pt2 = hylleraas_compressed_mp2(pt1_vec, ref_vec, cluster_ops, clustered_ham; tol=tol, do_pt=true) + end + + + # + # Compress FOIS + # norm1 = orth_dot(pt1_vec, pt1_vec) + dim1 = length(pt1_vec) + pt1_vec = compress(pt1_vec, thresh=thresh_foi) + # norm2 = orth_dot(pt1_vec, pt1_vec) + dim2 = length(pt1_vec) + @printf(" %-50s%10i → %-10i (thresh = %8.1e)\n", "FOIS Compressed from: ", dim1, dim2, thresh_foi) + #@printf(" %-50s%10.2e → %-10.2e (thresh = %8.1e)\n", "Norm of |1>: ",norm1, norm2, thresh_foi) + @printf(" %-50s", "Overlap between <1|0>: ") + ovlp = nonorth_dot(pt1_vec, ref_vec, verbose=0) + [@printf("%10.6f", ovlp[1])] + println() + # # Solve CEPA println() diff --git a/src/type_BSTstate.jl b/src/type_BSTstate.jl index a6030a2..e3050c2 100644 --- a/src/type_BSTstate.jl +++ b/src/type_BSTstate.jl @@ -148,7 +148,34 @@ function BSTstate(v::BSTstate{TT,N,RR}; T=TT, R=RR) where {TT,N,RR} return w #=}}}=# end +""" + +Constructor - create copy, of a particular root +# Arguments +- `v`: input `BSTstate` object +- `T`: data type of new state +- `R`: specific root in new state +# Returns +- `BSTstate` +""" +function BSTstate(v::BSTstate{T,N,R}, root) where {T,N,R} + #={{{=# + data = OrderedDict{FockConfig{N},OrderedDict{TuckerConfig{N},Tucker{T,N,1}} }() + + w = BSTstate{T,N,1}(v.clusters, data, v.p_spaces, v.q_spaces) + for (fock, tconfigs) in v.data + add_fockconfig!(w, fock) + for (tconfig, tuck) in tconfigs + core=(tuck.core[root],) + factors = ntuple(i->tuck.factors[i],N) + w[fock][tconfig] = Tucker{T, N, 1}(core, factors) + end + end + return w +#=}}}=# +end +# return Tucker{T,NN,R}(cores, ntuple(i->t.factors[i],NN)) """ BSTstate(clusters::Vector{MOCluster},