cleaned cepa in tpsci and spt pt2 version and qdpt pt2 functions are…

… ready to commit and pull request

src/tpsci_outer.jl

@@ -1,7 +1,6 @@
 using TimerOutputs
 using BlockDavidson
 using BenchmarkTools
-using Printf
 
 """
     build_full_H(ci_vector::TPSCIstate, cluster_ops, clustered_ham::ClusteredOperator)
@@ -1052,82 +1051,3 @@ end
 #    return eval.(a)
 #end
 
-function do_fois_ci(ref::TPSCIstate{T,N,R}, cluster_ops, clustered_ham;
-                    H0          = "Hcmf",
-                    max_iter    = 50,
-                    nbody       = 4,
-                    thresh_foi  = 1e-6,
-                    tol         = 1e-5,
-                    thresh_clip = 1e-6,
-                    threaded    =false,
-                    prescreen   = false,
-                    compress    = false,
-                    pt          =false,
-                    verbose     = true) where {T,N,R}
-    @printf("\n-------------------------------------------------------\n")
-    @printf(" Do CI in FOIS\n")
-    @printf("   H0                      = %-s\n", H0)
-    @printf("   thresh_foi              = %-8.1e\n", thresh_foi)
-    @printf("   nbody                   = %-i\n", nbody)
-    @printf("\n")
-    @printf("   Length of Reference     = %-i\n", length(ref))
-    @printf("\n-------------------------------------------------------\n")
-
-# 
-    # Solve variationally in reference space
-    ref_vec = deepcopy(ref)
-    @printf(" Solve zeroth-order problem. Dimension = %10i\n", length(ref_vec))
-    @time e0, ref_vec = tps_ci_direct(ref_vec, cluster_ops, clustered_ham, conv_thresh=tol)
-
-
-    #
-    # Get First order wavefunction
-    println()
-    println(" Compute FOIS. Reference space dim = ", length(ref_vec))
-    # pt1_vec= deepcopy(ref_vec)
-    # pt1_vec=matvec(pt1_vec)
-    if threaded == true
-        pt1_vec = open_matvec_thread(ref_vec, cluster_ops, clustered_ham, nbody=nbody, thresh=thresh_foi, prescreen=prescreen)
-    else
-        pt1_vec = open_matvec_serial(ref_vec, cluster_ops, clustered_ham, nbody=nbody, thresh=thresh_foi, prescreen=prescreen)
-    end
-    for i in 1:R
-        @printf("Arnab: %12.8f\n", sqrt.(orth_dot(pt1_vec,pt1_vec))[i])
-    end
-    project_out!(pt1_vec, ref)
-    # Compress FOIS
-    if compress==true
-        norm1 = sqrt.(orth_dot(pt1_vec, pt1_vec))
-        dim1 = length(pt1_vec)
-        clip!(pt1_vec, thresh=thresh_clip) #does clip! function do the compression? or have to write a compress function.
-        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)
-        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
-    end
-    for i in 1:R
-        @printf(" %-50s%10.6f\n", "Overlap between <1|0>: ", overlap(pt1_vec, ref_vec)[i])
-    end
-
-    add!(ref_vec, pt1_vec)
-    # Solve for first order wavefunction 
-    println(" Compute CI energy in the space = ", length(ref_vec))
-
-    eci, ref_vec = tps_ci_direct(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
-    if pt==true
-        e_pt2,pt1_vec= compute_pt1_wavefunction(ref_vec, cluster_ops, clustered_ham;  H0=H0,verbose=verbose)  
-        for i in 1:R
-            @printf(" E(PT2)  for %ith state   = %12.8f\n",i, e_pt2[i])
-        end
-    end
-    return eci, ref_vec 
-    # println("debugging")
-    # error()
-end
-

src/tpsci_pt2_energy.jl

@@ -282,3 +282,79 @@ function _sum_pt2(sig_v, e2, Hd, E0, R)
 end
 
 
+"""
+    compute_qdpt_energy2(ci_vector::TPSCIstate{T,N,R}, cluster_ops, clustered_ham::ClusteredOperator; 
+        nbody=4, 
+        H0="Hcmf",
+        E0=nothing, #pass in <0|H0|0>, or compute it
+        thresh_foi=1e-8, 
+        prescreen=true,
+        threaded = true,
+        verbose=1) where {T,N,R}
+        this function computes the second-order energy correction using the Quasidegenerate perturbation theory
+    args::
+    TPSCIstate{T,N,R} : TPSCIstate object
+    cluster_ops::Dict{Int64,Dict{String,Any}} : cluster operators
+    clustered_ham::ClusteredOperator : clustered Hamiltonian
+    nbody::Int64 : number of clusters
+    H0::String : zeroth-order Hamiltonian
+    E0::Union{Nothing,Vector{T}} : <0|H0|0>
+    thresh_foi::Float64 : threshold for first-order interaction space
+    prescreen::Bool : prescreening
+    threaded::Bool : use threading
+    verbose::Int64 : verbosity level
+    Return:
+    corrected_energies::Vector{T} : Pt2 corrected energies
+"""
+
+function compute_qdpt_energy(ci_vector_in::TPSCIstate{T,N,R}, cluster_ops, clustered_ham::ClusteredOperator;
+                             nbody=4, 
+                             H0="Hcmf",
+                             E0=nothing, # pass in <0|H0|0>, or compute it
+                             thresh_foi=1e-8, 
+                             prescreen=true,
+                             threaded = true,
+                             verbose=1) where {T,N,R}
+    println(" |..................................do  QDPT......................................")
+    println(" thresh_foi    :", thresh_foi) 
+    println(" prescreen     :", prescreen) 
+    println(" H0            :", H0) 
+    println(" nbody         :", nbody) 
+
+    ci_vector = deepcopy(ci_vector_in)
+    orthonormalize!(ci_vector)
+
+    if E0 == nothing
+        println(" Compute <0|H|0>:")
+        E0 = compute_expectation_value_parallel(ci_vector, cluster_ops, clustered_ham)
+    end
+
+    if threaded == true
+        @time sig = open_matvec_thread(ci_vector, cluster_ops, clustered_ham, nbody=nbody, thresh=thresh_foi, prescreen=prescreen)
+    else
+        sig = open_matvec_serial(ci_vector, cluster_ops, clustered_ham, nbody=nbody, thresh=thresh_foi, prescreen=prescreen)
+    end
+    project_out!(sig, ci_vector)
+    # Here, `a` and `b` index into the model space configurations in `ci_vector`.
+    # `x` indexes into the external space configurations in `sig`.
+    println(" Compute diagonal elements of fois ")
+    @time Hd = compute_diagonal(sig, cluster_ops, clustered_ham)
+    H_ax = build_full_H_parallel(ci_vector, sig,  cluster_ops, clustered_ham)
+    # H_xb=build_full_H_parallel(sig, ci_vector,  cluster_ops, clustered_ham)
+    H_xb=H_ax'
+    H2 = zeros(size(H_ax)[1], size(H_xb)[2])
+    # display(H2)
+    n = length(Hd)
+    I = Diagonal(ones(n))
+    H2=(H_ax * pinv((I*E0[1])*I - Diagonal(Hd))* H_xb)
+    # Add zeroth-order Hamiltonian (H^0) contributions
+    H0=build_full_H_parallel(ci_vector, ci_vector,  cluster_ops, clustered_ham)
+    H_total = H2 + H0
+    # Diagonalize the resulting Hamiltonian to get the corrected energies
+    corrected_energies = eigen(H_total).values
+    @printf(" %5s %12s %12s\n", "Root", "E(0)", "E(2)")
+    for r in 1:R
+        @printf(" %5s %12.8f %12.8f\n", r, E0[r], corrected_energies[r])
+    end
+    return corrected_energies
+end

src/type_TPSCIstate.jl

@@ -581,27 +581,7 @@ function orthonormalize!(s::TPSCIstate{T,N,R}) where {T,N,R}
     set_vector!(s,v0)
 end
 #=}}}=#
-"""
-    orth_dot(ts1::FermiCG.TPSCIstate, ts2::FermiCG.TPSCIstate)
-
-Dot product between `ts2` and `ts1`
 
-"""
-function orth_dot(ts1::TPSCIstate{T,N,R}, ts2::TPSCIstate{T,N,R}) where {T,N,R}
-    #={{{=#
-    overlap = zeros(T,R)
-    for (fock,configs) in ts2.data
-        haskey(ts1, fock) || continue
-        for (config,coeffs) in configs
-            haskey(ts1[fock], config) || continue
-            for r in 1:R
-                overlap[r] += sum(ts1[fock][config][r] .* ts2[fock][config][r])
-            end
-        end
-    end
-    return overlap
-    #=}}}=#
-end
 
 """
     clip!(s::TPSCIstate; thresh=1e-5)
@@ -670,27 +650,6 @@ function extract_roots(v::TPSCIstate{T,N,R}, roots) where {T,N,R}
     return out
 end
 
-"""
-    function extract_chosen_root(v::TPSCIstate{T,N,R}, root::Int)
-
-Extracts the chosen root to give a new `TPSCIstate` with only that root.
-"""
-function extract_chosen_root(v::TPSCIstate{T,N,R}, root::Int) where {T,N,R}
-    if root < 1 || root > R
-        throw(ArgumentError("Root index $root is out of bounds. It should be between 1 and $R."))
-    end
-    out = TPSCIstate(v.clusters, T=T, R=1)
-    for (fock, configs) in v.data
-        add_fockconfig!(out, fock)
-        # display(configs)
-        for (config, coeffs) in configs
-            # display(coeffs)
-            out[fock][config] =MVector{1, T}(coeffs[root])
-        end
-    end
-
-    return out
-end
 
 
 """