diff --git a/docs/src/installation_instructions.md b/docs/src/installation_instructions.md index 245b6f35..cb31eb07 100644 --- a/docs/src/installation_instructions.md +++ b/docs/src/installation_instructions.md @@ -5,21 +5,15 @@ A Julia package for course-grained electronic structure calculations ### Installation 1. Download -```julia +``` git clone https://github.com/nmayhall-vt/FermiCG.git -cd FermiCG/ +cd FermiCG ``` - - -2. Create python virtual environment which will hold the PYSCF executable - -```julia -cd src/python -virtualenv -p python3 venv -source venv/bin/activate -pip install -r requirements.txt -cd ../../ -julia --project=./ -julia> using Pkg; Pkg.build("PyCall") +2. Testing +``` +julia --project=./ -tauto +using FermiCG +using Pkg +Pkg.test() ``` diff --git a/examples/tetracene_dimer/rhf.ipynb b/examples/tetracene_dimer/rhf.ipynb new file mode 100644 index 00000000..ef0b8d42 --- /dev/null +++ b/examples/tetracene_dimer/rhf.ipynb @@ -0,0 +1,327 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Initial guess from minao.\n", + "\n", + "\n", + "******** ********\n", + "method = RHF\n", + "initial guess = minao\n", + "damping factor = 0\n", + "level_shift factor = 0\n", + "DIIS = \n", + "diis_start_cycle = 1\n", + "diis_space = 8\n", + "SCF conv_tol = 1e-09\n", + "SCF conv_tol_grad = None\n", + "SCF max_cycles = 200\n", + "direct_scf = True\n", + "direct_scf_tol = 1e-13\n", + "chkfile to save SCF result = /var/folders/jm/nd85t1m57dv98qh_jtqr2clh0000gn/T/tmpdnsov0fw\n", + "max_memory 4000 MB (current use 0 MB)\n", + "Set gradient conv threshold to 3.16228e-05\n", + "Initial guess from minao.\n", + "init E= -1373.8028556562\n", + " HOMO = -0.155360683388173 LUMO = -0.0462263741041608\n", + "cycle= 1 E= -1360.24317894862 delta_E= 13.6 |g|= 0.719 |ddm|= 11.3\n", + " HOMO = -0.143564650973908 LUMO = 0.138616575609906\n", + "cycle= 2 E= -1360.45923594589 delta_E= -0.216 |g|= 0.214 |ddm|= 1.43\n", + " HOMO = -0.16857409289756 LUMO = 0.12156954691401\n", + "cycle= 3 E= -1360.47370291225 delta_E= -0.0145 |g|= 0.126 |ddm|= 0.515\n", + " HOMO = -0.161685188896432 LUMO = 0.130732383404331\n", + "cycle= 4 E= -1360.4797695373 delta_E= -0.00607 |g|= 0.0139 |ddm|= 0.195\n", + " HOMO = -0.161052562766623 LUMO = 0.132842939686475\n", + "cycle= 5 E= -1360.47997484479 delta_E= -0.000205 |g|= 0.007 |ddm|= 0.0363\n", + " HOMO = -0.161719238589572 LUMO = 0.13333514441787\n", + "cycle= 6 E= -1360.48004918981 delta_E= -7.43e-05 |g|= 0.0028 |ddm|= 0.0275\n", + " HOMO = -0.162116206364872 LUMO = 0.13340292370013\n", + "cycle= 7 E= -1360.48005930686 delta_E= -1.01e-05 |g|= 0.000758 |ddm|= 0.0121\n", + " HOMO = -0.162172259276158 LUMO = 0.133406742943205\n", + "cycle= 8 E= -1360.48005985105 delta_E= -5.44e-07 |g|= 0.000404 |ddm|= 0.00231\n", + " HOMO = -0.162160553391006 LUMO = 0.133425925701453\n", + "cycle= 9 E= -1360.480060043 delta_E= -1.92e-07 |g|= 0.000209 |ddm|= 0.00137\n", + " HOMO = -0.162173349764571 LUMO = 0.133407030042586\n", + "cycle= 10 E= -1360.48006011093 delta_E= -6.79e-08 |g|= 7.63e-05 |ddm|= 0.0011\n", + " HOMO = -0.16217164410856 LUMO = 0.133406153044998\n", + "cycle= 11 E= -1360.48006011748 delta_E= -6.56e-09 |g|= 2.19e-05 |ddm|= 0.000359\n", + " HOMO = -0.162171319204681 LUMO = 0.133405720724885\n", + "cycle= 12 E= -1360.48006011796 delta_E= -4.82e-10 |g|= 6.26e-06 |ddm|= 9.42e-05\n", + " HOMO = -0.162171197951934 LUMO = 0.133405667622305\n", + "Extra cycle E= -1360.480060118 delta_E= -3.27e-11 |g|= 4.12e-06 |ddm|= 1.48e-05\n", + "converged SCF energy = -1360.480060118\n" + ] + } + ], + "source": [ + "from functools import reduce\n", + "import numpy as np\n", + "import scipy\n", + "\n", + "import pyscf\n", + "from pyscf import fci\n", + "from pyscf import gto, scf, ao2mo, lo, tdscf, cc\n", + "\n", + "\n", + "def tda_denisty_matrix(td, state_id):\n", + " '''\n", + " Taking the TDA amplitudes as the CIS coefficients, calculate the density\n", + " matrix (in AO basis) of the excited states\n", + " '''\n", + " cis_t1 = td.xy[state_id][0]\n", + " dm_oo =-np.einsum('ia,ka->ik', cis_t1.conj(), cis_t1)\n", + " dm_vv = np.einsum('ia,ic->ac', cis_t1, cis_t1.conj())\n", + "\n", + " # The ground state density matrix in mo_basis\n", + " mf = td._scf\n", + " dm = np.diag(mf.mo_occ)\n", + "\n", + " # Add CIS contribution\n", + " nocc = cis_t1.shape[0]\n", + " # Note that dm_oo and dm_vv correspond to spin-up contribution. \"*2\" to\n", + " # include the spin-down contribution\n", + " dm[:nocc,:nocc] += dm_oo * 2\n", + " dm[nocc:,nocc:] += dm_vv * 2\n", + "\n", + " # Transform density matrix to AO basis\n", + " mo = mf.mo_coeff\n", + " dm = np.einsum('pi,ij,qj->pq', mo, dm, mo.conj())\n", + " return dm\n", + "\n", + "mol = gto.Mole()\n", + "mol.atom = '''\n", + "C 3.86480 -1.02720 1.27180\n", + "C 3.60170 -0.34020 0.07040\n", + "C 2.93000 -1.06360 2.30730\n", + "C 4.54410 -0.28950 -0.97720\n", + "C 2.32170 0.33090 -0.09430\n", + "C 3.20030 -1.73210 3.54060\n", + "C 1.65210 -0.39240 2.14450\n", + "C 4.27330 0.38740 -2.16700\n", + "C 2.05960 1.01870 -1.29550\n", + "C 1.37970 0.28170 0.95340\n", + "C 2.27260 -1.73270 4.55160\n", + "C 0.71620 -0.41900 3.22490\n", + "C 2.99590 1.05910 -2.32920\n", + "C 5.20900 0.41800 -3.24760\n", + "C 1.01870 -1.06880 4.39470\n", + "C 2.72630 1.73410 -3.55840\n", + "C 4.90670 1.07300 -4.41500\n", + "C 3.65410 1.73950 -4.56900\n", + "C -3.86650 1.01870 -1.29550\n", + "C -2.93020 1.05910 -2.32920\n", + "C -3.60440 0.33090 -0.09430\n", + "C -3.19980 1.73410 -3.55840\n", + "C -1.65280 0.38740 -2.16700\n", + "C -2.32430 -0.34020 0.07040\n", + "C -4.54630 0.28170 0.95340\n", + "C -2.27200 1.73950 -4.56900\n", + "C -0.71710 0.41800 -3.24760\n", + "C -1.38200 -0.28950 -0.97720\n", + "C -2.06130 -1.02720 1.27180\n", + "C -4.27400 -0.39240 2.14450\n", + "C -1.01930 1.07300 -4.41500\n", + "C -2.99610 -1.06360 2.30730\n", + "C -5.20980 -0.41900 3.22490\n", + "C -2.72580 -1.73210 3.54060\n", + "C -4.90730 -1.06880 4.39470\n", + "C -3.65350 -1.73270 4.55160\n", + "H 4.82300 -1.53290 1.39770\n", + "H 5.49910 -0.80290 -0.85660\n", + "H 4.15900 -2.23700 3.66390\n", + "H 1.10180 1.52560 -1.42170\n", + "H 0.42460 0.79440 0.83100\n", + "H 2.50000 -2.24040 5.48840\n", + "H -0.23700 0.09640 3.10140\n", + "H 6.16210 -0.09790 -3.12730\n", + "H 0.29870 -1.07700 5.21470\n", + "H 1.76850 2.24120 -3.67870\n", + "H 5.62580 1.08320 -5.23570\n", + "H 3.42730 2.25190 -5.50340\n", + "H -4.82430 1.52560 -1.42170\n", + "H -4.15760 2.24120 -3.67870\n", + "H -5.50150 0.79440 0.83100\n", + "H -2.49880 2.25190 -5.50340\n", + "H 0.23610 -0.09790 -3.12730\n", + "H -0.42700 -0.80290 -0.85660\n", + "H -1.10300 -1.53290 1.39770\n", + "H -0.30030 1.08320 -5.23570\n", + "H -6.16310 0.09640 3.10140\n", + "H -1.76710 -2.23700 3.66390\n", + "H -5.62740 -1.07700 5.21470\n", + "H -3.42610 -2.24040 5.48840\n", + "'''\n", + "\n", + "\n", + "mol.basis = 'sto-3g'\n", + "mol.spin = 0\n", + "mol.build()\n", + "\n", + "#mf = scf.ROHF(mol).x2c()\n", + "mf = scf.RHF(mol)\n", + "mf.verbose = 4\n", + "mf.get_init_guess(mol, key='minao')\n", + "mf.conv_tol = 1e-9\n", + "#mf.level_shift = .1\n", + "#mf.diis_start_cycle = 4\n", + "#mf.diis_space = 10\n", + "mf.run(max_cycle=200)\n", + "\n", + "\n", + "n_triplets = 2\n", + "n_singlets = 2\n", + "\n", + "avg_rdm1 = mf.make_rdm1()\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "\n", + "******** for ********\n", + "nstates = 2 singlet\n", + "deg_eia_thresh = 1.000e-03\n", + "wfnsym = None\n", + "conv_tol = 1e-09\n", + "eigh lindep = 1e-12\n", + "eigh level_shift = 0\n", + "eigh max_space = 50\n", + "eigh max_cycle = 100\n", + "chkfile = /var/folders/jm/nd85t1m57dv98qh_jtqr2clh0000gn/T/tmpdnsov0fw\n", + "max_memory 4000 MB (current use 0 MB)\n", + "\n", + "\n", + "Excited State energies (eV)\n", + "[4.37757462 4.50864348]\n", + "\n", + "** Singlet excitation energies and oscillator strengths **\n", + "Excited State 1: 4.37757 eV 283.23 nm f=0.7064\n", + " 117 -> 123 -0.10590\n", + " 118 -> 124 -0.10581\n", + " 119 -> 122 -0.45930\n", + " 120 -> 121 -0.49153\n", + "Excited State 2: 4.50864 eV 274.99 nm f=0.0000\n", + " 119 -> 121 0.47837\n", + " 120 -> 122 0.46727\n", + "\n", + "** Transition electric dipole moments (AU) **\n", + "state X Y Z Dip. S. Osc.\n", + " 1 2.3095 -1.1112 0.1355 6.5870 0.7064\n", + " 2 -0.0015 0.0006 0.0007 0.0000 0.0000\n", + "\n", + "** Transition velocity dipole moments (imaginary part, AU) **\n", + "state X Y Z Dip. S. Osc.\n", + " 1 -0.1016 0.0712 -0.0291 0.0162 0.0673\n", + " 2 0.0001 0.0000 0.0001 0.0000 0.0000\n", + "\n", + "** Transition magnetic dipole moments (imaginary part, AU) **\n", + "state X Y Z\n", + " 1 -0.0003 -0.0000 0.0022\n", + " 2 -0.0411 -0.1064 0.4355\n" + ] + } + ], + "source": [ + "# compute singlets\n", + "mytd = tdscf.TDA(mf)\n", + "mytd.singlet = True \n", + "mytd = mytd.run(nstates=n_singlets)\n", + "mytd.analyze()\n", + "for i in range(mytd.nroots):\n", + " avg_rdm1 += tda_denisty_matrix(mytd, i)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# compute triplets \n", + "mytd = tdscf.TDA(mf)\n", + "mytd.singlet = False \n", + "mytd = mytd.run(nstates=n_triplets)\n", + "mytd.analyze()\n", + "for i in range(mytd.nroots):\n", + " avg_rdm1 += tda_denisty_matrix(mytd, i)\n", + "\n", + "# normalize\n", + "avg_rdm1 = avg_rdm1 / (n_singlets + n_triplets + 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "S = mf.get_ovlp()\n", + "F = mf.get_fock()\n", + "np.save(\"fock_mat18\", F)\n", + "np.save(\"overlap_mat18\", S)\n", + "np.save(\"density_mat18\", mf.make_rdm1())\n", + "np.save(\"rhf_mo_coeffs18\", mf.mo_coeff)\n", + "np.save(\"cis_sa_density_mat18\", avg_rdm1)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/tetracene_dimer/spt_direct.jl b/examples/tetracene_dimer/spt_direct.jl index 444e932f..8b4d1894 100644 --- a/examples/tetracene_dimer/spt_direct.jl +++ b/examples/tetracene_dimer/spt_direct.jl @@ -39,24 +39,9 @@ FermiCG.add_double_excitons!(ci_vector,FermiCG.FockConfig(init_fspace),nroots) fspace_0 = FermiCG.FockConfig(init_fspace) # FermiCG.add_1electron_transfers!(ci_vector, fspace_0, 1) FermiCG.add_spin_flip_states!(ci_vector, fspace_0,1) -# # Spin-flip states -# ## ba -# tmp_fspace = FermiCG.replace(fspace_0, (1,2), ([4,2],[2,4])) -# FermiCG.add_fockconfig!(ci_vector, tmp_fspace) -# ci_vector[tmp_fspace][FermiCG.TuckerConfig((1:1,1:1))]=FermiCG.Tucker(tuple([zeros(Float64, 1, 1) for _ in 1:nroots]...)) -# ## ab -# tmp_fspace = FermiCG.replace(fspace_0, (1,2), ([2,4],[4,2])) -# FermiCG.add_fockconfig!(ci_vector, tmp_fspace) -# ci_vector[tmp_fspace][FermiCG.TuckerConfig((1:1,1:1))]=FermiCG.Tucker(tuple([zeros(Float64, 1, 1) for _ in 1:nroots]...)) display(ci_vector.data) FermiCG.eye!(ci_vector) display(ci_vector) -# σ = FermiCG.build_compressed_1st_order_state(ci_vector, cluster_ops, clustered_ham, -# nbody=4, -# thresh=1e-3) -# σ = FermiCG.compress(σ, thresh=1e-5) -# v2 = BSTstate(σ,R=10) -# FermiCG.eye!(v2) e_ci, v2 = FermiCG.ci_solve(ci_vector, cluster_ops, clustered_ham); e_var, v_var = FermiCG.block_sparse_tucker(v2, cluster_ops, clustered_ham, max_iter = 200, @@ -70,6 +55,8 @@ e_var, v_var = FermiCG.block_sparse_tucker(v2, cluster_ops, clustered_ham, resolve_ss = false, tol_tucker = 1e-4, solver = "davidson") +@time ept2 = FermiCG.compute_pt2_energy(v_var, cluster_ops, clustered_ham, thresh_foi=1e-6,prescreen = true,compress_twice = true) +@time ept2 = FermiCG.compute_pt2_energy2(v_var, cluster_ops, clustered_ham, thresh_foi=1e-6,prescreen = true,compress_twice = true) e_var, v_var = FermiCG.block_sparse_tucker(v_var, cluster_ops, clustered_ham, max_iter = 200, nbody = 4, diff --git a/src/tpsci_pt2_energy.jl b/src/tpsci_pt2_energy.jl index 2d3db0d7..9ae527cd 100644 --- a/src/tpsci_pt2_energy.jl +++ b/src/tpsci_pt2_energy.jl @@ -128,7 +128,11 @@ function compute_pt2_energy(ci_vector_in::TPSCIstate{T,N,R}, cluster_ops, cluste flush(stdout) + #tmp = Int(round(length(jobs_vec)/100)) tmp = Int(round(length(jobs_vec)/100)) + if tmp == 0 + tmp += 1 + end verbose < 1 || println(" |----------------------------------------------------------------------------------------------------|") verbose < 1 || println(" |0% 100%|") verbose < 1 || print(" |") @@ -278,3 +282,80 @@ 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 + Approximation: Taking ∂E_S=0 in denominator of the perturbation theory expression + 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 \ No newline at end of file diff --git a/src/tucker_outer.jl b/src/tucker_outer.jl index 78ee87a2..014c7983 100644 --- a/src/tucker_outer.jl +++ b/src/tucker_outer.jl @@ -257,11 +257,11 @@ After solving, the Energy can be obtained as: """ function tucker_cepa_solve(ref_vector::BSTstate{T,N,R}, cepa_vector::BSTstate, cluster_ops, clustered_ham, cepa_shift="cepa", - cepa_mit = 50; - tol=1e-5, - cache=true, - max_iter=30, - verbose=false) where {T,N,R} + cepa_mit = 50; + tol =1e-5, + cache =true, + max_iter =30, + verbose =false) where {T,N,R} #={{{=# sig = deepcopy(ref_vector) @@ -270,7 +270,7 @@ function tucker_cepa_solve(ref_vector::BSTstate{T,N,R}, cepa_vector::BSTstate, c e0 = nonorth_dot(ref_vector, sig) length(e0) == 1 || error("Only one state at a time please", e0) e0 = e0[1] - @printf(" Reference Energy: %12.8f\n",e0[1]) + @printf(" Reference Energy: %12.8f\n",e0) n_clusters = length(cepa_vector.clusters) @@ -332,29 +332,33 @@ function tucker_cepa_solve(ref_vector::BSTstate{T,N,R}, cepa_vector::BSTstate, c Ec = 0 Ecepa = 0 - if cepa_shift == "cepa" - cepa_mit = 1 - end + # if cepa_shift == "cepa" + + # end for it in 1:cepa_mit - + println("CEPA cycle: ", it) bv = -get_vector(b) #n_clusters = 8 if cepa_shift == "cepa" - shift = 0.0 - elseif cepa_shift == "acpf" - - shift = Ec * 2.0 / n_clusters - elseif cepa_shift == "aqcc" - shift = (1.0 - (n_clusters-3.0)*(n_clusters - 2.0)/(n_clusters * ( n_clusters-1.0) )) * Ec - elseif cepa_shift == "cisd" - shift = Ec - else - println() - println("NYI: cepa_shift is not available:",cepa_shift) - println() - exit() - end - eshift = e0+shift + cepa_mit = 1 + shift = 0.0 + elseif cepa_shift == "acpf" + + shift = Ec * 2.0 / n_clusters + elseif cepa_shift == "aqcc" + shift = (1.0 - (n_clusters-3.0)*(n_clusters - 2.0)/(n_clusters * ( n_clusters-1.0) )) * Ec + elseif cepa_shift == "cisd" + shift = Ec + else + println() + println("NYI: cepa_shift is not available:",cepa_shift) + println() + exit() + end + + eshift = e0+shift + # display(shift) + # display(eshift) bv .= bv .+ get_vector(Sx)* (eshift) function mymatvec(v) @@ -428,10 +432,10 @@ 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)) - break - end + if abs(Ec - (Ecepa-e0)) < 1e-10 + @printf(" Converged %s %18.12f\n",cepa_shift, (e0[1] + ecorr[1])/(1+SxC[1])) + break + end Ec = Ecepa - e0 end @@ -1213,13 +1217,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 +1246,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 +1357,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 +1373,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/tucker_pt.jl b/src/tucker_pt.jl index a35b4685..95fec3bc 100644 --- a/src/tucker_pt.jl +++ b/src/tucker_pt.jl @@ -344,7 +344,11 @@ function form_1body_operator_diagonal!(sig::BSTstate{T,N,R}, Fdiag::BSTstate{T,N # # Rotate the original vector into this basis - transform_basis!(sig[fock][tconfig], rotations) + if length(rotations)==N + transform_basis!(sig[fock][tconfig], rotations) + else + continue + end end end @@ -723,7 +727,11 @@ function compute_pt2_energy(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; push!(e2_thread, zeros(T,R)) end - tmp = ceil(length(jobs_vec)/100) + #tmp = ceil(length(jobs_vec)/100) + tmp = Int(round(length(jobs_vec)/100)) + if tmp == 0 + tmp += 1 + end verbose < 2 || println(" |----------------------------------------------------------------------------------------------------|") verbose < 2 || println(" |0% 100%|") verbose < 2 || print(" |") @@ -927,5 +935,283 @@ function _pt2_job(sig_fock, job, ket::BSTstate{T,N,R}, cluster_ops, clustered_ha return ecorr end +""" + compute_pt2_energy2(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; + H0 = "Hcmf", + nbody = 4, + thresh_foi = 1e-6, + max_number = nothing, + opt_ref = true, + ci_tol = 1e-6, + verbose = true) where {T,N,R} + +Directly compute the PT2 energy, in parallel, with each thread handling a single `FockConfig` at a time. +""" +function compute_pt2_energy2(ref::BSTstate{T,N,R}, cluster_ops, clustered_ham; + H0 = "Hcmf", + nbody = 4, + thresh_foi = 1e-6, + max_number = nothing, + opt_ref = true, + ci_tol = 1e-6, + verbose = 1, + prescreen = false, + compress_twice = false) where {T,N,R} + println() + println(" |...................................BST-PT2............................................") + verbose < 1 || println(" H0 : ", H0 ) + verbose < 1 || println(" nbody : ", nbody ) + verbose < 1 || println(" thresh_foi : ", thresh_foi ) + verbose < 1 || println(" max_number : ", max_number ) + verbose < 1 || println(" opt_ref : ", opt_ref ) + verbose < 1 || println(" ci_tol : ", ci_tol ) + verbose < 1 || println(" verbose : ", verbose ) + verbose < 1 || @printf("\n") + verbose < 1 || @printf(" %-50s", "Length of Reference: ") + verbose < 1 || @printf("%10i\n", length(ref)) + + lk = ReentrantLock() + + # + # Solve variationally in reference space + ref_vec = deepcopy(ref) + clusters = ref_vec.clusters + + E0 = zeros(T,R) + + if opt_ref + @printf(" %-50s\n", "Solve zeroth-order problem: ") + time = @elapsed E0, ref_vec = ci_solve(ref_vec, cluster_ops, clustered_ham, conv_thresh=ci_tol) + @printf(" %-50s%10.6f seconds\n", "Diagonalization time: ",time) + else + @printf(" %-50s", "Compute zeroth-order energy: ") + flush(stdout) + @time E0 = compute_expectation_value(ref_vec, cluster_ops, clustered_ham) + end + + # + # get E0 = <0|H0|0> + clustered_ham_0 = extract_1body_operator(clustered_ham, op_string = H0) + @printf(" %-50s", "Compute <0|H0|0>: ") + @time F0 = compute_expectation_value(ref_vec, cluster_ops, clustered_ham_0) + + if verbose > 0 + @printf(" %5s %12s %12s\n", "Root", "<0|H|0>", "<0|F|0>") + for r in 1:R + @printf(" %5s %12.8f %12.8f\n",r, E0[r], F0[r]) + end + end + + # + # define batches (FockConfigs present in resolvant) + jobs = Dict{FockConfig{N},Vector{Tuple}}() + for (fock_ket, configs_ket) in ref_vec.data + for (ftrans, terms) in clustered_ham + fock_x = ftrans + fock_ket + + # + # check to make sure this fock config doesn't have negative or too many electrons in any cluster + all(f[1] >= 0 for f in fock_x) || continue + all(f[2] >= 0 for f in fock_x) || continue + all(f[1] <= length(clusters[fi]) for (fi,f) in enumerate(fock_x)) || continue + all(f[2] <= length(clusters[fi]) for (fi,f) in enumerate(fock_x)) || continue + + job_input = (terms, fock_ket, configs_ket) + if haskey(jobs, fock_x) + push!(jobs[fock_x], job_input) + else + jobs[fock_x] = [job_input] + end + end + end + + + jobs_vec = [] + for (fock_x, job) in jobs + push!(jobs_vec, (fock_x, job)) + end + + println(" Number of jobs: ", length(jobs_vec)) + println(" Number of threads: ", Threads.nthreads()) + BLAS.set_num_threads(1) + flush(stdout) + + #ham_0s = Vector{ClusteredOperator}() + #for t in Threads.nthreads() + # push!(ham_0s, extract_1body_operator(clustered_ham, op_string = H0) ) + #end + + + e2_thread = Vector{Vector{Float64}}() + for tid in 1:Threads.nthreads() + push!(e2_thread, zeros(T,R)) + end + + #tmp = ceil(length(jobs_vec)/100) + tmp = Int(round(length(jobs_vec)/100)) + if tmp == 0 + tmp += 1 + end + verbose < 2 || println(" |----------------------------------------------------------------------------------------------------|") + verbose < 2 || println(" |0% 100%|") + verbose < 2 || print(" |") + #@profilehtml @Threads.threads for job in jobs_vec + nprinted = 0 + alloc = @allocated t = @elapsed begin + + @Threads.threads for (jobi,job) in collect(enumerate(jobs_vec)) + #for (jobi,job) in collect(enumerate(jobs_vec)) + fock_sig = job[1] + tid = Threads.threadid() + e2_thread[tid] .+= _pt2_job2(fock_sig, job[2], ref_vec, cluster_ops, clustered_ham, clustered_ham_0, + nbody, verbose, thresh_foi, max_number, E0, F0, prescreen, compress_twice) + if verbose > 1 + if jobi%tmp == 0 + begin + lock(lk) + try + print("-") + nprinted += 1 + flush(stdout) + finally + unlock(lk) + end + end + end + end + end + end + flush(stdout) + verbose < 2 || for i in nprinted+1:100 + print("-") + end + verbose < 2 || println("|") + flush(stdout) + + @printf(" %-48s%10.1f s Allocated: %10.1e GB\n", "Time spent computing E2: ",t,alloc*1e-9) + ecorr = sum(e2_thread) + + E2 = zeros(R) + for r in 1:R + E2[r] = E0[r] + ecorr[r] + @printf(" State %3i: %-35s%14.8f\n", r, "E(PT2) corr: ", ecorr[r]) + end + + @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], E2[r]) + end + println(" ......................................................................................|") + + return E2 +end + + +function _pt2_job2(sig_fock, job, ket::BSTstate{T,N,R}, cluster_ops, clustered_ham, clustered_ham_0, + nbody, verbose, thresh, max_number, E0, F0, prescreen, compress_twice) where {T,N,R} + + tconfigs_to_process = Dict{TuckerConfig{N}, Vector{Vector{Any}}}() + + ecorr = Vector{T}([0.0 for i in 1:R]) + + for jobi in job + + terms, ket_fock, ket_tconfigs = jobi + + for term in terms + + length(term.clusters) <= nbody || continue + for (ket_tconfig, ket_tuck) in ket_tconfigs + available = [] + for ci in term.clusters + tmp = [] + if haskey(ket.p_spaces[ci.idx], sig_fock[ci.idx]) + push!(tmp, ket.p_spaces[ci.idx][sig_fock[ci.idx]]) + end + if haskey(ket.q_spaces[ci.idx], sig_fock[ci.idx]) + push!(tmp, ket.q_spaces[ci.idx][sig_fock[ci.idx]]) + end + push!(available, tmp) + end + for prod in Iterators.product(available...) + sig_tconfig = [ket_tconfig.config...] + for cidx in 1:length(term.clusters) + ci = term.clusters[cidx] + sig_tconfig[ci.idx] = prod[cidx] + end + sig_tconfig = TuckerConfig(sig_tconfig) + # println(sig_tconfig) + # println([term, ket_fock, ket_tconfig]) + if haskey(tconfigs_to_process, sig_tconfig) + push!(tconfigs_to_process[sig_tconfig], [term, ket_fock, ket_tconfig]) + else + tconfigs_to_process[sig_tconfig] = [[term, ket_fock, ket_tconfig]] + end + end + end + end + end + curr_tuck = Vector{Any}() + for (sig_tconfig, terms_to_process) in tconfigs_to_process + #curr_tuck=Tucker{T,N,R}() #initialization problem is there as in every iteration there is a new dimension of Tucker core and factor + #if I don't want to store tucker in a list, though it is not affecting the time and memory allocations + for (term, ket_fock, ket_tconfig) in terms_to_process + ket_tuck = ket[ket_fock][ket_tconfig] + check_term(term, sig_fock, sig_tconfig, ket_fock, ket_tconfig) || continue + + + if prescreen + bound = calc_bound(term, cluster_ops, + sig_fock, sig_tconfig, + ket_fock, ket_tconfig, ket_tuck, + prescreen=thresh) + bound == true || continue + end + + sig_tuck = form_sigma_block_expand(term, cluster_ops, + sig_fock, sig_tconfig, + ket_fock, ket_tconfig, ket_tuck, + max_number=max_number, + prescreen=thresh) + #compress new addition + sig_tuck = compress(sig_tuck, thresh=thresh) + + # println(sig_tuck) + if length(curr_tuck) == 0 + curr_tuck = sig_tuck + else + curr_tuck=nonorth_add([curr_tuck, sig_tuck]) + + end + ##curr_tuck=nonorth_add([curr_tuck, sig_tuck]) + end + + if length(curr_tuck) == 0 + continue + end + if norm(curr_tuck) < thresh + continue + end + + + #compress new addition + curr_tuck = compress(curr_tuck, thresh=thresh) + if compress_twice + curr_tuck = compress(curr_tuck, thresh=thresh) + end + + sig = BSTstate(ket.clusters, ket.p_spaces, ket.q_spaces, T=T, R=R) + add_fockconfig!(sig, sig_fock) + sig[sig_fock][sig_tconfig] = curr_tuck + + # Compute PT2 energy for this job + _, _, ecorr_i = compute_pt1_wavefunction(sig, ket, cluster_ops, clustered_ham, clustered_ham_0, E0, F0, verbose=0) + # println(ecorr_i) + ecorr += ecorr_i + end + + + return ecorr +end diff --git a/src/type_BSTstate.jl b/src/type_BSTstate.jl index e6d56829..dd3d50eb 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}, @@ -654,7 +681,7 @@ end Pretty print """ -function print_fock_occupations(s::BSTstate; thresh=1e-3) +function print_fock_occupations(s::BSTstate; thresh=1e-3,root=1) #={{{=# println() @@ -947,4 +974,46 @@ function project_into_new_basis(v1::BSTstate{T,N,R}, v2::BSTstate{T,N,R}) where end end return out -end \ No newline at end of file +end + +""" + ct_table(s::BSTstate; ne_cluster=10, nroots=1) + +Prints total weight of charge transfer in each root in table formate +# Arguments +- `s::BSTstate` +- `ne_cluster`: Int, number of total electrons in each cluster +- `nroots`: Total number of roots +""" +function ct_table(s::BSTstate{T,N,R}; ne_cluster=10, nroots=1) where {T,N,R} + @printf(" -----------------------\n") + @printf(" --- CHARGE TRANSFER ---\n") + @printf(" -----------------------\n") + @printf(" %-15s%-10s\n", "Root", "Total CT") + @printf(" %-15s%-10s\n", "-------", "---------") + for root in 1:nroots + ct = 0 + for (fock,configs) in s.data + # println(fock) + prob = 0 + is_ct = false + + for cluster in 1:length(s.clusters) + if sum(fock[cluster]) != ne_cluster + is_ct = true + end + end + if is_ct + prob = 0 + for (tconfig, tucker) in configs + # println(tucker.core) + # println(tucker.factors) + prob += sum(tucker.core[root] .^ 2) + end + end + ct += prob + end + @printf(" %-15i%-10.5f", root, ct) + println() + end +end diff --git a/test/_testdata_cmf_cr2_19.jld2 b/test/_testdata_cmf_cr2_19.jld2 new file mode 100644 index 00000000..9fffda5a Binary files /dev/null and b/test/_testdata_cmf_cr2_19.jld2 differ diff --git a/test/runtests.jl b/test/runtests.jl index ab59eab7..0a39c348 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,5 +16,6 @@ Random.seed!(1234567) include("test_bs.jl") include("test_schmidt.jl") include("test_openshell.jl") + include("test_qdpt.jl") end diff --git a/test/test_qdpt.jl b/test/test_qdpt.jl new file mode 100644 index 00000000..29da1115 --- /dev/null +++ b/test/test_qdpt.jl @@ -0,0 +1,41 @@ +using QCBase +using RDM +using FermiCG +using Printf +using Test +using LinearAlgebra +using Profile +using JLD2 + +# load data +@load "_testdata_cmf_cr2_19.jld2" + +M=20 +init_fspace = FockConfig(init_fspace) + +cluster_bases = FermiCG.compute_cluster_eigenbasis_spin(ints, clusters, d1, [3,3,3,3,3], init_fspace, max_roots=M, verbose=1); + +clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters); +cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints); + +FermiCG.add_cmf_operators!(cluster_ops, cluster_bases, ints, d1.a, d1.b); + +nroots=4 + + +ci_vector = FermiCG.TPSCIstate(clusters, init_fspace, R=nroots) + +ci_vector = FermiCG.add_spin_focksectors(ci_vector) + +eci, v = FermiCG.tps_ci_direct(ci_vector, cluster_ops, clustered_ham); + + +e2a = FermiCG.compute_pt2_energy(v, cluster_ops, clustered_ham) +e_qdpt=FermiCG.compute_qdpt_energy(v, cluster_ops, clustered_ham, nbody=4, H0="Hcmf", thresh_foi=1e-8, verbose=1, threaded=true) + +println(" PT2 - ROCMF-PT2") +display(e2a) +println(" PT2 - ROCMF-QDPT") +display(e_qdpt) +@test isapprox(e_qdpt, [-108.13989647019216, -108.13975457069576, -108.1394893503325, -108.1390882840833], atol=1e-12) +