From 57e592870cb3fa3e4cc650433adc7f28398d9983 Mon Sep 17 00:00:00 2001 From: eginer Date: Wed, 5 Oct 2022 17:01:27 +0200 Subject: [PATCH] added tc_bi_ortho --- src/tc_bi_ortho/EZFIO.cfg | 11 + src/tc_bi_ortho/NEED | 6 + src/tc_bi_ortho/e_corr_bi_ortho.irp.f | 104 +++++ src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f | 92 +++++ src/tc_bi_ortho/normal_ordered.irp.f | 319 ++++++++++++++++ src/tc_bi_ortho/print_tc_wf.irp.f | 104 +++++ src/tc_bi_ortho/psi_det_tc_sorted.irp.f | 157 ++++++++ src/tc_bi_ortho/psi_left_qmc.irp.f | 44 +++ src/tc_bi_ortho/psi_r_l_prov.irp.f | 217 +++++++++++ .../save_bitcpsileft_for_qmcchem.irp.f | 61 +++ src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f | 15 + src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f | 33 ++ src/tc_bi_ortho/select_dets_bi_ortho.irp.f | 61 +++ src/tc_bi_ortho/slater_tc.irp.f | 359 ++++++++++++++++++ src/tc_bi_ortho/slater_tc_3e.irp.f | 293 ++++++++++++++ src/tc_bi_ortho/symmetrized_3_e_int.irp.f | 111 ++++++ src/tc_bi_ortho/tc_bi_ortho.irp.f | 61 +++ src/tc_bi_ortho/tc_bi_ortho_prop.irp.f | 24 ++ src/tc_bi_ortho/tc_cisd_sc2.irp.f | 24 ++ src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f | 110 ++++++ src/tc_bi_ortho/tc_h_eigvectors.irp.f | 179 +++++++++ src/tc_bi_ortho/tc_hmat.irp.f | 41 ++ src/tc_bi_ortho/tc_prop.irp.f | 268 +++++++++++++ src/tc_bi_ortho/test_normal_order.irp.f | 73 ++++ src/tc_bi_ortho/test_tc_bi_ortho.irp.f | 131 +++++++ src/tc_bi_ortho/test_tc_fock.irp.f | 169 +++++++++ 26 files changed, 3067 insertions(+) create mode 100644 src/tc_bi_ortho/EZFIO.cfg create mode 100644 src/tc_bi_ortho/NEED create mode 100644 src/tc_bi_ortho/e_corr_bi_ortho.irp.f create mode 100644 src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f create mode 100644 src/tc_bi_ortho/normal_ordered.irp.f create mode 100644 src/tc_bi_ortho/print_tc_wf.irp.f create mode 100644 src/tc_bi_ortho/psi_det_tc_sorted.irp.f create mode 100644 src/tc_bi_ortho/psi_left_qmc.irp.f create mode 100644 src/tc_bi_ortho/psi_r_l_prov.irp.f create mode 100644 src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f create mode 100644 src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f create mode 100644 src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f create mode 100644 src/tc_bi_ortho/select_dets_bi_ortho.irp.f create mode 100644 src/tc_bi_ortho/slater_tc.irp.f create mode 100644 src/tc_bi_ortho/slater_tc_3e.irp.f create mode 100644 src/tc_bi_ortho/symmetrized_3_e_int.irp.f create mode 100644 src/tc_bi_ortho/tc_bi_ortho.irp.f create mode 100644 src/tc_bi_ortho/tc_bi_ortho_prop.irp.f create mode 100644 src/tc_bi_ortho/tc_cisd_sc2.irp.f create mode 100644 src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f create mode 100644 src/tc_bi_ortho/tc_h_eigvectors.irp.f create mode 100644 src/tc_bi_ortho/tc_hmat.irp.f create mode 100644 src/tc_bi_ortho/tc_prop.irp.f create mode 100644 src/tc_bi_ortho/test_normal_order.irp.f create mode 100644 src/tc_bi_ortho/test_tc_bi_ortho.irp.f create mode 100644 src/tc_bi_ortho/test_tc_fock.irp.f diff --git a/src/tc_bi_ortho/EZFIO.cfg b/src/tc_bi_ortho/EZFIO.cfg new file mode 100644 index 00000000..a34d2134 --- /dev/null +++ b/src/tc_bi_ortho/EZFIO.cfg @@ -0,0 +1,11 @@ +[psi_l_coef_bi_ortho] +interface: ezfio +doc: Coefficients for the left wave function +type: double precision +size: (determinants.n_det,determinants.n_states) + +[psi_r_coef_bi_ortho] +interface: ezfio +doc: Coefficients for the right wave function +type: double precision +size: (determinants.n_det,determinants.n_states) diff --git a/src/tc_bi_ortho/NEED b/src/tc_bi_ortho/NEED new file mode 100644 index 00000000..9a0c20ef --- /dev/null +++ b/src/tc_bi_ortho/NEED @@ -0,0 +1,6 @@ +bi_ort_ints +bi_ortho_mos +tc_keywords +non_hermit_dav +dav_general_mat +tc_scf diff --git a/src/tc_bi_ortho/e_corr_bi_ortho.irp.f b/src/tc_bi_ortho/e_corr_bi_ortho.irp.f new file mode 100644 index 00000000..ec66a8b5 --- /dev/null +++ b/src/tc_bi_ortho/e_corr_bi_ortho.irp.f @@ -0,0 +1,104 @@ + use bitmasks ! you need to include the bitmasks_module.f90 features + BEGIN_PROVIDER [ double precision, e_tilde_00] + implicit none + double precision :: hmono,htwoe,hthree,htot + call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,htot) + e_tilde_00 = htot + END_PROVIDER + + BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth] +&BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth_single] +&BEGIN_PROVIDER [ double precision, e_pt2_tc_bi_orth_double] + implicit none + integer :: i,degree + double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e + e_pt2_tc_bi_orth = 0.d0 + e_pt2_tc_bi_orth_single = 0.d0 + e_pt2_tc_bi_orth_double = 0.d0 + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree == 1 .or. degree == 2)then + call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0) + delta_e = e_tilde_00 - e_i0 + coef_pt1 = htilde_ij / delta_e + call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) + e_pt2_tc_bi_orth += coef_pt1 * htilde_ij + if(degree == 1)then + e_pt2_tc_bi_orth_single += coef_pt1 * htilde_ij + else +! print*,'coef_pt1, e_pt2',coef_pt1,coef_pt1 * htilde_ij + e_pt2_tc_bi_orth_double += coef_pt1 * htilde_ij + endif + endif + enddo + END_PROVIDER + + BEGIN_PROVIDER [ double precision, e_tilde_bi_orth_00] + implicit none + double precision :: hmono,htwoe,hthree,htilde_ij + call htilde_mu_mat_bi_ortho(HF_bitmask,HF_bitmask,N_int,hmono,htwoe,hthree,e_tilde_bi_orth_00) + e_tilde_bi_orth_00 += nuclear_repulsion + END_PROVIDER + + BEGIN_PROVIDER [ double precision, e_corr_bi_orth ] +&BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj ] +&BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth ] +&BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth ] + implicit none + integer :: i,degree + double precision :: hmono,htwoe,hthree,htilde_ij + + e_corr_bi_orth = 0.d0 + e_corr_single_bi_orth = 0.d0 + e_corr_double_bi_orth = 0.d0 + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) + if(degree == 1)then + e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) + else if(degree == 2)then + e_corr_double_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) +! print*,'coef_wf , e_cor',reigvec_tc_bi_orth(i,1)/reigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) + endif + enddo + e_corr_bi_orth_proj = e_corr_single_bi_orth + e_corr_double_bi_orth + e_corr_bi_orth = eigval_right_tc_bi_orth(1) - e_tilde_bi_orth_00 + END_PROVIDER + + BEGIN_PROVIDER [ double precision, e_tc_left_right ] + implicit none + integer :: i,j + double precision :: hmono,htwoe,hthree,htilde_ij,accu + e_tc_left_right = 0.d0 + accu = 0.d0 + do i = 1, N_det + accu += reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(i,1) + do j = 1, N_det + call htilde_mu_mat_bi_ortho(psi_det(1,1,j),psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) + e_tc_left_right += htilde_ij * reigvec_tc_bi_orth(i,1) * leigvec_tc_bi_orth(j,1) + enddo + enddo + e_tc_left_right *= 1.d0/accu + e_tc_left_right += nuclear_repulsion + + END_PROVIDER + + +BEGIN_PROVIDER [ double precision, coef_pt1_bi_ortho, (N_det)] + implicit none + integer :: i,degree + double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree==0)then + coef_pt1_bi_ortho(i) = 1.d0 + else + call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0) + delta_e = e_tilde_00 - e_i0 + coef_pt1 = htilde_ij / delta_e + coef_pt1_bi_ortho(i)= coef_pt1 + endif + enddo +END_PROVIDER diff --git a/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f new file mode 100644 index 00000000..b7129d36 --- /dev/null +++ b/src/tc_bi_ortho/h_tc_bi_ortho_psi.irp.f @@ -0,0 +1,92 @@ +subroutine htc_bi_ortho_calc_tdav(v, u, N_st, sze) + + use bitmasks + + BEGIN_DOC + ! Application of H_TC on a vector + ! + ! v(i,istate) = \sum_j u(j,istate) H_TC(i,j), with: + ! H_TC(i,j) = < Di | H_TC | Dj > + ! + END_DOC + + implicit none + + integer, intent(in) :: N_st, sze + double precision, intent(in) :: u(sze,N_st) + double precision, intent(inout) :: v(sze,N_st) + + integer :: i, j, istate + double precision :: htot + + PROVIDE N_int + PROVIDE psi_det + + + ! TODO : transform it with the bi-linear representation in terms of alpha-beta. + + i = 1 + j = 1 + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot) + + v = 0.d0 + !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & + !$OMP SHARED(N_st, sze, N_int, psi_det, u, v) & + !$OMP PRIVATE(istate, i, j, htot) + do istate = 1, N_st + do i = 1, sze + do j = 1, sze + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot) + v(i,istate) = v(i,istate) + htot * u(j,istate) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end + +subroutine htcdag_bi_ortho_calc_tdav(v, u, N_st, sze) + + use bitmasks + + BEGIN_DOC + ! Application of (H_TC)^dagger on a vector + ! + ! v(i,istate) = \sum_j u(j,istate) H_TC(j,i), with: + ! H_TC(i,j) = < Di | H_TC | Dj > + ! + END_DOC + + implicit none + + integer, intent(in) :: N_st, sze + double precision, intent(in) :: u(sze,N_st) + double precision, intent(inout) :: v(sze,N_st) + + integer :: i, j, istate + double precision :: htot + + PROVIDE N_int + PROVIDE psi_det + + i = 1 + j = 1 + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,j), N_int, htot) + + v = 0.d0 + + !$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(dynamic,8) & + !$OMP SHARED(N_st, sze, N_int, psi_det, u, v) & + !$OMP PRIVATE(istate, i, j, htot) + do istate = 1, N_st + do i = 1, sze + do j = 1, sze + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot) + v(i,istate) = v(i,istate) + htot * u(j,istate) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end + diff --git a/src/tc_bi_ortho/normal_ordered.irp.f b/src/tc_bi_ortho/normal_ordered.irp.f new file mode 100644 index 00000000..81f5fb2c --- /dev/null +++ b/src/tc_bi_ortho/normal_ordered.irp.f @@ -0,0 +1,319 @@ +BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! Normal ordering of the three body interaction on the HF density + END_DOC + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + + integer :: i,h1,p1,h2,p2 + integer :: hh1,hh2,pp1,pp2 + integer :: Ne(2) + integer, allocatable :: occ(:,:) + integer(bit_kind), allocatable :: key_i_core(:,:) + double precision :: hthree_aba,hthree_aaa,hthree_aab + double precision :: wall0,wall1 + + PROVIDE N_int + + allocate( occ(N_int*bit_kind_size,2) ) + allocate( key_i_core(N_int,2) ) + + if(core_tc_op) then + do i = 1, N_int + key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1)) + key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core,occ,Ne,N_int) + else + call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + endif + + normal_two_body_bi_orth = 0.d0 + print*,'Providing normal_two_body_bi_orth ...' + call wall_time(wall0) + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aba, hthree_aab, hthree_aaa) & + !$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth) + !$OMP DO SCHEDULE (static) + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + do hh2 = 1, n_act_orb + h2 = list_act(hh2) + do pp2 = 1, n_act_orb + p2 = list_act(pp2) + ! opposite spin double excitations + call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aba) + ! same spin double excitations with opposite spin contributions + if(h1h2 + ! same spin double excitations with same spin contributions + if(Ne(2).ge.3)then + call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2 + else + hthree_aaa = 0.d0 + endif + else + call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab) + if(Ne(2).ge.3)then + call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa) + else + hthree_aaa = 0.d0 + endif + endif + normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aba + hthree_aab + hthree_aaa) + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + + call wall_time(wall1) + print*,'Wall time for normal_two_body_bi_orth ',wall1-wall0 + + deallocate( occ ) + deallocate( key_i_core ) + +END_PROVIDER + + + +subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + integer, intent(in) :: Nint, h1, h2, p1, p2 + integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) + double precision, intent(out) :: hthree + integer :: ii, i + double precision :: int_direct, int_exc_12, int_exc_13, integral + + !!!! double alpha/beta + hthree = 0.d0 + do ii = 1, Ne(2) ! purely closed shell part + i = occ(ii,2) + call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral) + int_direct = -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral) + int_exc_13 = -1.d0 * integral + call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral) + int_exc_12 = -1.d0 * integral + hthree += 2.d0 * int_direct - 1.d0 * ( int_exc_13 + int_exc_12) + enddo + do ii = Ne(2) + 1, Ne(1) ! purely open-shell part + i = occ(ii,1) + call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral) + int_direct = -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral) + int_exc_13 = -1.d0 * integral + call give_integrals_3_body_bi_ort(p2, i,p1,i,h2,h1,integral) + int_exc_12 = -1.d0 * integral + hthree += 1.d0 * int_direct - 0.5d0* ( int_exc_13 + int_exc_12) + enddo + +end subroutine give_aba_contraction + + + +BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_ab, (mo_num, mo_num, mo_num, mo_num)] + + BEGIN_DOC + ! Normal ordered two-body sector of the three-body terms for opposite spin double excitations + END_DOC + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + integer :: h1, p1, h2, p2, i + integer :: hh1, hh2, pp1, pp2 + integer :: Ne(2) + integer, allocatable :: occ(:,:) + integer(bit_kind), allocatable :: key_i_core(:,:) + double precision :: hthree + + PROVIDE N_int + + allocate( key_i_core(N_int,2) ) + allocate( occ(N_int*bit_kind_size,2) ) + + if(core_tc_op)then + do i = 1, N_int + key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1)) + key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core,occ,Ne,N_int) + else + call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int) + endif + normal_two_body_bi_orth_ab = 0.d0 + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + do pp1 = 1, n_act_orb + p1 = list_act(pp1) + do hh2 = 1, n_act_orb + h2 = list_act(hh2) + do pp2 = 1, n_act_orb + p2 = list_act(pp2) + call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree) + normal_two_body_bi_orth_ab(p2,h2,p1,h1) = hthree + enddo + enddo + enddo + enddo + + deallocate( key_i_core ) + deallocate( occ ) + +END_PROVIDER + + + +BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth_aa_bb, (n_act_orb, n_act_orb, n_act_orb, n_act_orb)] + + BEGIN_DOC + ! Normal ordered two-body sector of the three-body terms for same spin double excitations + END_DOC + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + integer :: i,ii,j,h1,p1,h2,p2 + integer :: hh1,hh2,pp1,pp2 + integer :: Ne(2) + integer, allocatable :: occ(:,:) + integer(bit_kind), allocatable :: key_i_core(:,:) + double precision :: hthree_aab, hthree_aaa + + PROVIDE N_int + + allocate( key_i_core(N_int,2) ) + allocate( occ(N_int*bit_kind_size,2) ) + + if(core_tc_op)then + do i = 1, N_int + key_i_core(i,1) = xor(ref_bitmask(i,1),core_bitmask(i,1)) + key_i_core(i,2) = xor(ref_bitmask(i,2),core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core, occ, Ne, N_int) + else + call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int) + endif + + normal_two_body_bi_orth_aa_bb = 0.d0 + do hh1 = 1, n_act_orb + h1 = list_act(hh1) + do pp1 = 1 , n_act_orb + p1 = list_act(pp1) + do hh2 = 1, n_act_orb + h2 = list_act(hh2) + do pp2 = 1 , n_act_orb + p2 = list_act(pp2) + if(h1h2 + if(Ne(2).ge.3)then + call give_aaa_contraction(N_int, h2, h1, p1, p2, Ne, occ, hthree_aaa) ! exchange h1<->h2 + else + hthree_aaa = 0.d0 + endif + else + call give_aab_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aab) + if(Ne(2).ge.3)then + call give_aaa_contraction(N_int, h1, h2, p1, p2, Ne, occ, hthree_aaa) + else + hthree_aaa = 0.d0 + endif + endif + normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) = hthree_aab + hthree_aaa + enddo + enddo + enddo + enddo + + deallocate( key_i_core ) + deallocate( occ ) + +END_PROVIDER + + + +subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + integer, intent(in) :: Nint, h1, h2, p1, p2 + integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) + double precision, intent(out) :: hthree + integer :: ii,i + double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23 + double precision :: integral,int_exc_l,int_exc_ll + + hthree = 0.d0 + do ii = 1, Ne(2) ! purely closed shell part + i = occ(ii,2) + call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral) + int_direct = -1.d0 * integral + call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral) + int_exc_l = -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral) + int_exc_ll= -1.d0 * integral + call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral) + int_exc_12= -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral) + int_exc_13= -1.d0 * integral + call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral) + int_exc_23= -1.d0 * integral + + hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 ) + enddo + do ii = Ne(2)+1,Ne(1) ! purely open-shell part + i = occ(ii,1) + call give_integrals_3_body_bi_ort(i ,p2,p1,i,h2,h1,integral) + int_direct = -1.d0 * integral + call give_integrals_3_body_bi_ort(p2,p1,i ,i,h2,h1,integral) + int_exc_l = -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,i ,p2,i,h2,h1,integral) + int_exc_ll= -1.d0 * integral + call give_integrals_3_body_bi_ort(p2,i ,p1,i,h2,h1,integral) + int_exc_12= -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,p2, i,i,h2,h1,integral) + int_exc_13= -1.d0 * integral + call give_integrals_3_body_bi_ort(i ,p1,p2,i,h2,h1,integral) + int_exc_23= -1.d0 * integral + + hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll -( int_exc_12+ int_exc_13+ int_exc_23 )) + enddo + +end subroutine give_aaa_contraction + + + +subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree) + implicit none + use bitmasks ! you need to include the bitmasks_module.f90 features + integer, intent(in) :: Nint, h1, h2, p1, p2 + integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2) + double precision, intent(out) :: hthree + integer :: ii, i + double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23 + double precision :: integral, int_exc_l, int_exc_ll + + hthree = 0.d0 + do ii = 1, Ne(2) ! purely closed shell part + i = occ(ii,2) + call give_integrals_3_body_bi_ort(p2,p1,i,h2,h1,i,integral) + int_direct = -1.d0 * integral + call give_integrals_3_body_bi_ort(p1,p2,i,h2,h1,i,integral) + int_exc_23= -1.d0 * integral + hthree += 1.d0 * int_direct - int_exc_23 + enddo + +end subroutine give_aab_contraction diff --git a/src/tc_bi_ortho/print_tc_wf.irp.f b/src/tc_bi_ortho/print_tc_wf.irp.f new file mode 100644 index 00000000..58a733a7 --- /dev/null +++ b/src/tc_bi_ortho/print_tc_wf.irp.f @@ -0,0 +1,104 @@ +program print_tc_bi_ortho + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid +! if(three_body_h_tc)then +! call provide_all_three_ints_bi_ortho +! endif +! call routine + call write_l_r_wf +end + +subroutine write_l_r_wf + implicit none + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.tc_wf' + i_unit_output = getUnitAndOpen(output,'w') + integer :: i + print*,'Writing the left-right wf' + do i = 1, N_det + write(i_unit_output,*)i,psi_l_coef_sorted_bi_ortho_left(i),psi_r_coef_sorted_bi_ortho_right(i) + enddo + + +end + +subroutine routine + implicit none + integer :: i,degree + integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2 + double precision :: hmono,htwoe,hthree,htilde_ij,coef_pt1,e_i0,delta_e,e_pt2 + double precision :: contrib_pt,e_corr,coef,contrib,phase + double precision :: accu_positive,accu_positive_pt, accu_positive_core,accu_positive_core_pt + e_pt2 = 0.d0 + accu_positive = 0.D0 + accu_positive_pt = 0.D0 + accu_positive_core = 0.d0 + accu_positive_core_pt = 0.d0 + + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree == 1 .or. degree == 2)then + call htilde_mu_mat_bi_ortho(psi_det(1,1,i),HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call htilde_mu_mat_bi_ortho(psi_det(1,1,i),psi_det(1,1,i),N_int,hmono,htwoe,hthree,e_i0) + delta_e = e_tilde_00 - e_i0 + coef_pt1 = htilde_ij / delta_e + + call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) + contrib_pt = coef_pt1 * htilde_ij + e_pt2 += contrib_pt + + coef = psi_r_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1) + contrib = coef * htilde_ij + e_corr += contrib + call get_excitation(HF_bitmask,psi_det(1,1,i),exc,degree,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + print*,'*********' + if(degree==1)then + print*,'s1',s1 + print*,'h1,p1 = ',h1,p1 + else if(degree ==2)then + print*,'s1',s1 + print*,'h1,p1 = ',h1,p1 + print*,'s2',s2 + print*,'h2,p2 = ',h2,p2 + endif + print*,'coef_pt1 = ',coef_pt1 + print*,'coef = ',coef + print*,'contrib_pt ',contrib_pt + print*,'contrib = ',contrib + if(contrib.gt.0.d0)then + accu_positive += contrib + if(h1==1.or.h2==1)then + accu_positive_core += contrib + endif + if(dabs(contrib).gt.1.d-5)then + print*,'Found a positive contribution to correlation energy !!' + endif + endif + if(contrib_pt.gt.0.d0)then + accu_positive_pt += contrib_pt + if(h2==1.or.h1==1)then + accu_positive_core_pt += contrib_pt + endif + endif + endif + enddo + print*,'' + print*,'' + print*,'Total correlation energy = ',e_corr + print*,'Total correlation energy PT = ',e_pt2 + print*,'Positive contribution to ecorr = ',accu_positive + print*,'Positive contribution to ecorr PT = ',accu_positive_pt + print*,'Pure core contribution = ',accu_positive_core + print*,'Pure core contribution PT = ',accu_positive_core_pt +end diff --git a/src/tc_bi_ortho/psi_det_tc_sorted.irp.f b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f new file mode 100644 index 00000000..212c8588 --- /dev/null +++ b/src/tc_bi_ortho/psi_det_tc_sorted.irp.f @@ -0,0 +1,157 @@ +use bitmasks + +BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_tc, (psi_det_size) ] + implicit none + BEGIN_DOC + ! Contribution of determinants to the state-averaged density. + END_DOC + integer :: i,j,k + double precision :: f + + psi_average_norm_contrib_tc(:) = 0.d0 + do k=1,N_states + do i=1,N_det + psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i) + & + dabs(psi_l_coef_bi_ortho(i,k)*psi_r_coef_bi_ortho(i,k))*state_average_weight(k) + enddo + enddo + f = 1.d0/sum(psi_average_norm_contrib_tc(1:N_det)) + do i=1,N_det + psi_average_norm_contrib_tc(i) = psi_average_norm_contrib_tc(i)*f + enddo +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc, (psi_det_size,N_states) ] +&BEGIN_PROVIDER [ double precision, psi_average_norm_contrib_sorted_tc, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, psi_det_sorted_tc_order, (psi_det_size) ] + implicit none + BEGIN_DOC + ! Wave function sorted by determinants contribution to the norm (state-averaged) + ! + ! psi_det_sorted_tc_order(i) -> k : index in psi_det + END_DOC + integer :: i,j,k + integer, allocatable :: iorder(:) + allocate ( iorder(N_det) ) + do i=1,N_det + psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_tc(i) + iorder(i) = i + enddo + call dsort(psi_average_norm_contrib_sorted_tc,iorder,N_det) + do i=1,N_det + do j=1,N_int + psi_det_sorted_tc(j,1,i) = psi_det(j,1,iorder(i)) + psi_det_sorted_tc(j,2,i) = psi_det(j,2,iorder(i)) + enddo + psi_average_norm_contrib_sorted_tc(i) = -psi_average_norm_contrib_sorted_tc(i) + psi_det_sorted_tc_order(iorder(i)) = i + enddo + double precision :: accu + do k=1,N_states + accu = 0.d0 + do i=1,N_det + psi_coef_sorted_tc(i,k) = dsqrt(dabs(psi_l_coef_bi_ortho(iorder(i),k)*psi_r_coef_bi_ortho(iorder(i),k))) + accu += psi_coef_sorted_tc(i,k)**2 + enddo + accu = 1.d0/dsqrt(accu) + do i=1,N_det + psi_coef_sorted_tc(i,k) *= accu + enddo + enddo + + psi_det_sorted_tc(:,:,N_det+1:psi_det_size) = 0_bit_kind + psi_coef_sorted_tc(N_det+1:psi_det_size,:) = 0.d0 + psi_average_norm_contrib_sorted_tc(N_det+1:psi_det_size) = 0.d0 + psi_det_sorted_tc_order(N_det+1:psi_det_size) = 0 + + deallocate(iorder) + +END_PROVIDER + + BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho, (psi_det_size, N_states)] +&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho, (psi_det_size, N_states)] + BEGIN_DOC + ! psi_r_coef_sorted_bi_ortho : right coefficients corresponding to psi_det_sorted_tc + ! psi_l_coef_sorted_bi_ortho : left coefficients corresponding to psi_det_sorted_tc + END_DOC + implicit none + integer :: i, j, k + psi_r_coef_sorted_bi_ortho = 0.d0 + psi_l_coef_sorted_bi_ortho = 0.d0 + do i = 1, N_det + psi_r_coef_sorted_bi_ortho(i,1) = psi_r_coef_bi_ortho(psi_det_sorted_tc_order(i),1) + psi_l_coef_sorted_bi_ortho(i,1) = psi_l_coef_bi_ortho(psi_det_sorted_tc_order(i),1) + enddo + +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_coef_sorted_tc_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! Determinants on which we apply $\langle i|H|psi \rangle$ for perturbation. + ! They are sorted by determinants interpreted as integers. Useful + ! to accelerate the search of a random determinant in the wave + ! function. + END_DOC + + call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, size(psi_coef,1), & + psi_det_sorted_tc_bit, psi_coef_sorted_tc_bit, N_states) + +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_right, (N_int,2,N_det) ] +&BEGIN_PROVIDER [double precision, psi_r_coef_sorted_bi_ortho_right, (N_det)] + implicit none + BEGIN_DOC + ! psi_det_sorted_tc_right : Slater determinants sorted by decreasing value of |right- coefficients| + ! + ! psi_r_coef_sorted_bi_ortho_right : right wave function according to psi_det_sorted_tc_right + END_DOC + integer, allocatable :: iorder(:) + double precision, allocatable :: coef(:) + integer :: i,j + allocate ( iorder(N_det) , coef(N_det)) + do i=1,N_det + coef(i) = -dabs(psi_r_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1)) + iorder(i) = i + enddo + call dsort(coef,iorder,N_det) + do i=1,N_det + do j=1,N_int + psi_det_sorted_tc_right(j,1,i) = psi_det(j,1,iorder(i)) + psi_det_sorted_tc_right(j,2,i) = psi_det(j,2,iorder(i)) + enddo + psi_r_coef_sorted_bi_ortho_right(i) = psi_r_coef_bi_ortho(iorder(i),1)/psi_r_coef_bi_ortho(iorder(1),1) + enddo +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_tc_left, (N_int,2,N_det) ] +&BEGIN_PROVIDER [double precision, psi_l_coef_sorted_bi_ortho_left, (N_det)] + implicit none + BEGIN_DOC + ! psi_det_sorted_tc_left : Slater determinants sorted by decreasing value of |LEFTt- coefficients| + ! + ! psi_r_coef_sorted_bi_ortho_left : LEFT wave function according to psi_det_sorted_tc_left + END_DOC + integer, allocatable :: iorder(:) + double precision, allocatable :: coef(:) + integer :: i,j + allocate ( iorder(N_det) , coef(N_det)) + do i=1,N_det + coef(i) = -dabs(psi_l_coef_bi_ortho(i,1)/psi_r_coef_bi_ortho(1,1)) + iorder(i) = i + enddo + call dsort(coef,iorder,N_det) + do i=1,N_det + do j=1,N_int + psi_det_sorted_tc_left(j,1,i) = psi_det(j,1,iorder(i)) + psi_det_sorted_tc_left(j,2,i) = psi_det(j,2,iorder(i)) + enddo + psi_l_coef_sorted_bi_ortho_left(i) = psi_l_coef_bi_ortho(iorder(i),1)/psi_l_coef_bi_ortho(iorder(1),1) + enddo +END_PROVIDER diff --git a/src/tc_bi_ortho/psi_left_qmc.irp.f b/src/tc_bi_ortho/psi_left_qmc.irp.f new file mode 100644 index 00000000..25048f82 --- /dev/null +++ b/src/tc_bi_ortho/psi_left_qmc.irp.f @@ -0,0 +1,44 @@ + +! --- + +BEGIN_PROVIDER [ double precision, psi_bitcleft_bilinear_matrix_values, (N_det,N_states) ] + + BEGIN_DOC + ! Sparse coefficient matrix if the wave function is expressed in a bilinear form : + ! $D_\alpha^\dagger.C.D_\beta$ + ! + ! Rows are $\alpha$ determinants and columns are $\beta$. + ! + ! Order refers to psi_det + END_DOC + + use bitmasks + + implicit none + integer :: k, l + + if(N_det .eq. 1) then + + do l = 1, N_states + psi_bitcleft_bilinear_matrix_values(1,l) = 1.d0 + enddo + + else + + do l = 1, N_states + do k = 1, N_det + psi_bitcleft_bilinear_matrix_values(k,l) = psi_l_coef_bi_ortho(k,l) + enddo + enddo + + PROVIDE psi_bilinear_matrix_order + do l = 1, N_states + call dset_order(psi_bitcleft_bilinear_matrix_values(1,l), psi_bilinear_matrix_order, N_det) + enddo + + endif + +END_PROVIDER + +! --- + diff --git a/src/tc_bi_ortho/psi_r_l_prov.irp.f b/src/tc_bi_ortho/psi_r_l_prov.irp.f new file mode 100644 index 00000000..551da858 --- /dev/null +++ b/src/tc_bi_ortho/psi_r_l_prov.irp.f @@ -0,0 +1,217 @@ +use bitmasks + +BEGIN_PROVIDER [ double precision, psi_l_coef_bi_ortho, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file + ! is empty. + END_DOC + + integer :: i,k, N_int2 + logical :: exists + character*(64) :: label + + PROVIDE read_wf N_det mo_label ezfio_filename nproc + psi_l_coef_bi_ortho = 0.d0 + do i=1,min(N_states,N_det) + psi_l_coef_bi_ortho(i,i) = 1.d0 + enddo + + if (mpi_master) then + if (read_wf) then + call ezfio_has_tc_bi_ortho_psi_l_coef_bi_ortho(exists) +! if (exists) then +! call ezfio_has_tc_bi_ortho_mo_label(exists) +! if (exists) then +! call ezfio_get_tc_bi_ortho_mo_label(label) +! exists = (label == mo_label) +! endif +! endif + + if (exists) then + + double precision, allocatable :: psi_l_coef_bi_ortho_read(:,:) + allocate (psi_l_coef_bi_ortho_read(N_det,N_states)) + print *, 'Read psi_l_coef_bi_ortho', N_det, N_states + call ezfio_get_tc_bi_ortho_psi_l_coef_bi_ortho(psi_l_coef_bi_ortho_read) + do k=1,N_states + do i=1,N_det + psi_l_coef_bi_ortho(i,k) = psi_l_coef_bi_ortho_read(i,k) + enddo + enddo + deallocate(psi_l_coef_bi_ortho_read) + + endif + endif + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( psi_l_coef_bi_ortho, size(psi_l_coef_bi_ortho), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read psi_l_coef_bi_ortho with MPI' + endif + IRP_ENDIF +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, psi_r_coef_bi_ortho, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! The wave function coefficients. Initialized with Hartree-Fock if the |EZFIO| file + ! is empty. + END_DOC + + integer :: i,k, N_int2 + logical :: exists + character*(64) :: label + + PROVIDE read_wf N_det mo_label ezfio_filename nproc + psi_r_coef_bi_ortho = 0.d0 + do i=1,min(N_states,N_det) + psi_r_coef_bi_ortho(i,i) = 1.d0 + enddo + + if (mpi_master) then + if (read_wf) then + call ezfio_has_tc_bi_ortho_psi_r_coef_bi_ortho(exists) +! if (exists) then +! call ezfio_has_tc_bi_ortho_mo_label(exists) +! if (exists) then +! call ezfio_get_tc_bi_ortho_mo_label(label) +! exists = (label == mo_label) +! endif +! endif + + if (exists) then + + double precision, allocatable :: psi_r_coef_bi_ortho_read(:,:) + allocate (psi_r_coef_bi_ortho_read(N_det,N_states)) + print *, 'Read psi_r_coef_bi_ortho', N_det, N_states + call ezfio_get_tc_bi_ortho_psi_r_coef_bi_ortho(psi_r_coef_bi_ortho_read) + do k=1,N_states + do i=1,N_det + psi_r_coef_bi_ortho(i,k) = psi_r_coef_bi_ortho_read(i,k) + enddo + enddo + deallocate(psi_r_coef_bi_ortho_read) + + endif + endif + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( psi_r_coef_bi_ortho, size(psi_r_coef_bi_ortho), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read psi_r_coef_bi_ortho with MPI' + endif + IRP_ENDIF +END_PROVIDER + + +subroutine save_tc_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psilcoef,psircoef) + implicit none + BEGIN_DOC + ! Save the wave function into the |EZFIO| file + END_DOC + use bitmasks + include 'constants.include.F' + integer, intent(in) :: ndet,nstates,dim_psicoef + integer(bit_kind), intent(in) :: psidet(N_int,2,ndet) + double precision, intent(in) :: psilcoef(dim_psicoef,nstates) + double precision, intent(in) :: psircoef(dim_psicoef,nstates) + integer*8, allocatable :: psi_det_save(:,:,:) + double precision, allocatable :: psil_coef_save(:,:) + double precision, allocatable :: psir_coef_save(:,:) + + double precision :: accu_norm + integer :: i,j,k, ndet_qp_edit + + if (mpi_master) then + ndet_qp_edit = min(ndet,N_det_qp_edit) + + call ezfio_set_determinants_N_int(N_int) + call ezfio_set_determinants_bit_kind(bit_kind) + call ezfio_set_determinants_N_det(ndet) + call ezfio_set_determinants_N_det_qp_edit(ndet_qp_edit) + call ezfio_set_determinants_n_states(nstates) + call ezfio_set_determinants_mo_label(mo_label) + + allocate (psi_det_save(N_int,2,ndet)) + do i=1,ndet + do j=1,2 + do k=1,N_int + psi_det_save(k,j,i) = transfer(psidet(k,j,i),1_8) + enddo + enddo + enddo + call ezfio_set_determinants_psi_det(psi_det_save) + call ezfio_set_determinants_psi_det_qp_edit(psi_det_save) + deallocate (psi_det_save) + + allocate (psil_coef_save(ndet,nstates),psir_coef_save(ndet,nstates)) + do k=1,nstates + do i=1,ndet + psil_coef_save(i,k) = psilcoef(i,k) + psir_coef_save(i,k) = psircoef(i,k) + enddo + enddo + + call ezfio_set_tc_bi_ortho_psi_l_coef_bi_ortho(psil_coef_save) + call ezfio_set_tc_bi_ortho_psi_r_coef_bi_ortho(psir_coef_save) + deallocate (psil_coef_save,psir_coef_save) + +! allocate (psi_coef_save(ndet_qp_edit,nstates)) +! do k=1,nstates +! do i=1,ndet_qp_edit +! psi_coef_save(i,k) = psicoef(i,k) +! enddo +! enddo +! +! call ezfio_set_determinants_psi_coef_qp_edit(psi_coef_save) +! deallocate (psi_coef_save) + + call write_int(6,ndet,'Saved determinantsi and psi_r/psi_l coef') + endif +end + +subroutine save_tc_bi_ortho_wavefunction + implicit none + call save_tc_wavefunction_general(N_det,N_states,psi_det,size(psi_l_coef_bi_ortho, 1),psi_l_coef_bi_ortho,psi_r_coef_bi_ortho) + call routine_save_right_bi_ortho +end + +subroutine routine_save_right_bi_ortho + implicit none + double precision, allocatable :: coef_tmp(:,:) + integer :: i + N_states = 1 + allocate(coef_tmp(N_det, N_states)) + do i = 1, N_det + coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1) + enddo + call save_wavefunction_general_unormalized(N_det,N_states,psi_det,size(coef_tmp,1),coef_tmp(1,1)) +end + +subroutine routine_save_left_right_bi_ortho + implicit none + double precision, allocatable :: coef_tmp(:,:) + integer :: i,n_states_tmp + n_states_tmp = 2 + allocate(coef_tmp(N_det, n_states_tmp)) + do i = 1, N_det + coef_tmp(i,1) = psi_r_coef_bi_ortho(i,1) + coef_tmp(i,2) = psi_l_coef_bi_ortho(i,1) + enddo + call save_wavefunction_general_unormalized(N_det,n_states_tmp,psi_det,size(coef_tmp,1),coef_tmp(1,1)) +end + diff --git a/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f new file mode 100644 index 00000000..60201f5f --- /dev/null +++ b/src/tc_bi_ortho/save_bitcpsileft_for_qmcchem.irp.f @@ -0,0 +1,61 @@ +program save_bitcpsileft_for_qmcchem + + integer :: iunit + logical :: exists + double precision :: e_ref + + print *, ' ' + print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + print *, ' call save_for_qmcchem before ' + print *, ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~' + print *, ' ' + + call write_lr_spindeterminants() + + e_ref = 0.d0 + iunit = 13 + open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write') + call ezfio_has_fci_energy_pt2(exists) + + if(.not.exists) then + call ezfio_has_fci_energy(exists) + + if(.not.exists) then + call ezfio_has_tc_scf_bitc_energy(exists) + if(exists) then + call ezfio_get_tc_scf_bitc_energy(e_ref) + endif + endif + + endif + write(iunit,*) e_ref + close(iunit) + +end + +! -- + +subroutine write_lr_spindeterminants() + + use bitmasks + + implicit none + + integer :: k, l + double precision, allocatable :: buffer(:,:) + + PROVIDE psi_bitcleft_bilinear_matrix_values + + allocate(buffer(N_det,N_states)) + do l = 1, N_states + do k = 1, N_det + buffer(k,l) = psi_bitcleft_bilinear_matrix_values(k,l) + enddo + enddo + call ezfio_set_spindeterminants_psi_left_coef_matrix_values(buffer) + deallocate(buffer) + +end subroutine write_lr_spindeterminants + +! --- + diff --git a/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f b/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f new file mode 100644 index 00000000..5eb3c069 --- /dev/null +++ b/src/tc_bi_ortho/save_lr_bi_ortho_states.irp.f @@ -0,0 +1,15 @@ +program tc_bi_ortho + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call routine_save_left_right_bi_ortho +! call test +end diff --git a/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f new file mode 100644 index 00000000..9b7b9f5a --- /dev/null +++ b/src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f @@ -0,0 +1,33 @@ + program tc_natorb_bi_ortho + implicit none + BEGIN_DOC + ! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call save_tc_natorb + end + + subroutine save_tc_natorb + implicit none + print*,'Saving the natorbs ' + provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao + call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao) + call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao) + call save_ref_determinant_nstates_1 + call ezfio_set_determinants_read_wf(.False.) + end + + subroutine save_ref_determinant_nstates_1 + implicit none + use bitmasks + double precision :: buffer(1,N_states) + buffer = 0.d0 + buffer(1,1) = 1.d0 + call save_wavefunction_general(1,1,ref_bitmask,1,buffer) + end diff --git a/src/tc_bi_ortho/select_dets_bi_ortho.irp.f b/src/tc_bi_ortho/select_dets_bi_ortho.irp.f new file mode 100644 index 00000000..e6bf3d6e --- /dev/null +++ b/src/tc_bi_ortho/select_dets_bi_ortho.irp.f @@ -0,0 +1,61 @@ +program tc_bi_ortho + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + !!!!!!!!!!!!!!! WARNING NO 3-BODY + !!!!!!!!!!!!!!! WARNING NO 3-BODY + three_body_h_tc = .False. + touch three_body_h_tc + !!!!!!!!!!!!!!! WARNING NO 3-BODY + !!!!!!!!!!!!!!! WARNING NO 3-BODY + + call routine_test +! call test +end + +subroutine routine_test + implicit none + use bitmasks ! you need to include the bitmasks_module.f90 features + integer :: i,n_good,degree + integer(bit_kind), allocatable :: dets(:,:,:) + integer, allocatable :: iorder(:) + double precision, allocatable :: coef(:),coef_new(:,:) + double precision :: thr + allocate(coef(N_det), iorder(N_det)) + do i = 1, N_det + iorder(i) = i + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree==1)then + coef(i) = -0.5d0 + else + coef(i) = -dabs(coef_pt1_bi_ortho(i)) + endif + enddo + call dsort(coef,iorder,N_det) + !thr = save_threshold + thr = 1d-15 + n_good = 0 + do i = 1, N_det + if(dabs(coef(i)).gt.thr)then + n_good += 1 + endif + enddo + print*,'n_good = ',n_good + allocate(dets(N_int,2,n_good),coef_new(n_good,n_states)) + do i = 1, n_good + dets(:,:,i) = psi_det(:,:,iorder(i)) + coef_new(i,:) = psi_coef(iorder(i),:) + enddo + call save_wavefunction_general(n_good,n_states,dets,n_good,coef_new) + + +end diff --git a/src/tc_bi_ortho/slater_tc.irp.f b/src/tc_bi_ortho/slater_tc.irp.f new file mode 100644 index 00000000..45115a40 --- /dev/null +++ b/src/tc_bi_ortho/slater_tc.irp.f @@ -0,0 +1,359 @@ +!!!!!! +subroutine htilde_mu_mat_bi_ortho_tot(key_j, key_i, Nint, htot) + + BEGIN_DOC + ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + double precision, intent(out) :: htot + double precision :: hmono,htwoe,hthree + integer :: degree + call get_excitation_degree(key_j, key_i, degree, Nint) + if(degree.gt.2)then + htot = 0.d0 + else + call htilde_mu_mat_bi_ortho(key_j,key_i, Nint, hmono,htwoe,hthree,htot) + endif + +end subroutine htilde_mu_mat_tot + +! -- + +subroutine htilde_mu_mat_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot) + implicit none + use bitmasks + BEGIN_DOC + ! where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis + !! + ! Returns the detail of the matrix element in terms of single, two and three electron contribution. + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2),key_j(Nint,2) + double precision, intent(out) :: hmono,htwoe,hthree,htot + integer :: degree + + hmono = 0.d0 + htwoe= 0.d0 + htot = 0.d0 + hthree = 0.D0 + call get_excitation_degree(key_i, key_j, degree, Nint) + if(degree.gt.2)return + if(degree == 0)then + call diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot) + else if (degree == 1)then + call single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + else if(degree == 2)then + call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + endif + if(three_body_h_tc) then + if(degree == 2) then + if(.not.double_normal_ord) then + call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + endif + else if(degree == 1)then + call single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + else if(degree == 0)then + call diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) + endif + endif + htot = hmono + htwoe + hthree + if(degree==0)then + htot += nuclear_repulsion + endif + +end + +subroutine diag_htilde_mu_mat_bi_ortho(Nint, key_i, hmono, htwoe, htot) + + BEGIN_DOC + ! diagonal element of htilde ONLY FOR ONE- AND TWO-BODY TERMS + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + double precision, intent(out) :: hmono,htwoe,htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + double precision :: get_mo_two_e_integral_tc_int + integer(bit_kind) :: key_i_core(Nint,2) + +! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e +! +! PROVIDE mo_integrals_erf_map core_energy nuclear_repulsion core_bitmask +! PROVIDE core_fock_operator +! +! PROVIDE j1b_gauss + +! if(core_tc_op)then +! print*,'core_tc_op not already taken into account for bi ortho' +! print*,'stopping ...' +! stop +! do i = 1, Nint +! key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) +! key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) +! enddo +! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) +! hmono = core_energy - nuclear_repulsion +! else + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + hmono = 0.d0 +! endif + htwoe= 0.d0 + htot = 0.d0 + + do ispin = 1, 2 + do i = 1, Ne(ispin) ! + ii = occ(i,ispin) + hmono += mo_bi_ortho_tc_one_e(ii,ii) + +! if(j1b_gauss .eq. 1) then +! print*,'j1b not implemented for bi ortho TC' +! print*,'stopping ....' +! stop +! !hmono += mo_j1b_gauss_hermI (ii,ii) & +! ! + mo_j1b_gauss_hermII (ii,ii) & +! ! + mo_j1b_gauss_nonherm(ii,ii) +! endif + +! if(core_tc_op)then +! print*,'core_tc_op not already taken into account for bi ortho' +! print*,'stopping ...' +! stop +! hmono += core_fock_operator(ii,ii) ! add the usual Coulomb - Exchange from the core +! endif + enddo + enddo + + + ! alpha/beta two-body + ispin = 1 + jspin = 2 + do i = 1, Ne(ispin) ! electron 1 (so it can be associated to mu(r1)) + ii = occ(i,ispin) + do j = 1, Ne(jspin) ! electron 2 + jj = occ(j,jspin) + htwoe += mo_bi_ortho_tc_two_e(jj,ii,jj,ii) + enddo + enddo + + ! alpha/alpha two-body + do i = 1, Ne(ispin) + ii = occ(i,ispin) + do j = i+1, Ne(ispin) + jj = occ(j,ispin) + htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii) + enddo + enddo + + ! beta/beta two-body + do i = 1, Ne(jspin) + ii = occ(i,jspin) + do j = i+1, Ne(jspin) + jj = occ(j,jspin) + htwoe += mo_bi_ortho_tc_two_e(ii,jj,ii,jj) - mo_bi_ortho_tc_two_e(ii,jj,jj,ii) + enddo + enddo + htot = hmono + htwoe + +end + + + +subroutine double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + + BEGIN_DOC + ! for double excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: hmono, htwoe, htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + integer :: other_spin(2) + integer(bit_kind) :: key_i_core(Nint,2) + double precision :: get_mo_two_e_integral_tc_int,phase + +! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e + + other_spin(1) = 2 + other_spin(2) = 1 + + call get_excitation_degree(key_i, key_j, degree, Nint) + + hmono = 0.d0 + htwoe= 0.d0 + htot = 0.d0 + + if(degree.ne.2)then + return + endif + +! if(core_tc_op)then +! print*,'core_tc_op not already taken into account for bi ortho' +! print*,'stopping ...' +! stop +! do i = 1, Nint +! key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) +! key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) +! enddo +! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) +! else + call bitstring_to_list_ab(key_i, occ, Ne, Nint) +! endif + call get_double_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) + + if(s1.ne.s2)then + ! opposite spin two-body +! key_j, key_i + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + if(double_normal_ord.and.+Ne(1).gt.2)then + htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ??? + endif + else + ! same spin two-body + ! direct terms + htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) + ! exchange terms + htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) + if(double_normal_ord.and.+Ne(1).gt.2)then + htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ??? + htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ??? + endif + endif + htwoe *= phase + htot = htwoe + +end + + +subroutine single_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot) + + BEGIN_DOC + ! for single excitation ONLY FOR ONE- AND TWO-BODY TERMS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2) + double precision, intent(out) :: hmono, htwoe, htot + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2), i, j, ii, jj, ispin, jspin, k, kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: get_mo_two_e_integral_tc_int, phase + double precision :: direct_int, exchange_int_12, exchange_int_23, exchange_int_13 + integer :: other_spin(2) + integer(bit_kind) :: key_j_core(Nint,2), key_i_core(Nint,2) + +! PROVIDE mo_two_e_integrals_tc_int_in_map mo_bi_ortho_tc_two_e +! +! PROVIDE core_bitmask core_fock_operator mo_integrals_erf_map + +! PROVIDE j1b_gauss + + other_spin(1) = 2 + other_spin(2) = 1 + + hmono = 0.d0 + htwoe= 0.d0 + htot = 0.d0 + call get_excitation_degree(key_i, key_j, degree, Nint) + if(degree.ne.1)then + return + endif +! if(core_tc_op)then +! print*,'core_tc_op not already taken into account for bi ortho' +! print*,'stopping ...' +! stop +! do i = 1, Nint +! key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) +! key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) +! key_j_core(i,1) = xor(key_j(i,1),core_bitmask(i,1)) +! key_j_core(i,2) = xor(key_j(i,2),core_bitmask(i,2)) +! enddo +! call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) +! else + call bitstring_to_list_ab(key_i, occ, Ne, Nint) +! endif + + call get_single_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc,1,h1,p1,h2,p2,s1,s2) + + hmono = mo_bi_ortho_tc_one_e(p1,h1) * phase + +! if(j1b_gauss .eq. 1) then +! print*,'j1b not implemented for bi ortho TC' +! print*,'stopping ....' +! stop +! !hmono += ( mo_j1b_gauss_hermI (h1,p1) & +! ! + mo_j1b_gauss_hermII (h1,p1) & +! ! + mo_j1b_gauss_nonherm(h1,p1) ) * phase +! endif + +! if(core_tc_op)then +! print*,'core_tc_op not already taken into account for bi ortho' +! print*,'stopping ...' +! stop +! hmono += phase * core_fock_operator(h1,p1) +! endif + + ! alpha/beta two-body + ispin = other_spin(s1) + if(s1==1)then + ! single alpha + do i = 1, Ne(ispin) ! electron 2 + ii = occ(i,ispin) + htwoe += mo_bi_ortho_tc_two_e(ii,p1,ii,h1) + enddo + else + ! single beta + do i = 1, Ne(ispin) ! electron 1 + ii = occ(i,ispin) + htwoe += mo_bi_ortho_tc_two_e(p1,ii,h1,ii) + enddo + endif +! ! same spin two-body + do i = 1, Ne(s1) + ii = occ(i,s1) + ! (h1p1|ii ii) - (h1 ii| p1 ii) + htwoe += mo_bi_ortho_tc_two_e(ii,p1,ii,h1) - mo_bi_ortho_tc_two_e(p1,ii,ii,h1) + enddo + + htwoe *= phase + htot = hmono + htwoe + +end + + diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f new file mode 100644 index 00000000..7c2c9c86 --- /dev/null +++ b/src/tc_bi_ortho/slater_tc_3e.irp.f @@ -0,0 +1,293 @@ +subroutine provide_all_three_ints_bi_ortho + implicit none + BEGIN_DOC +! routine that provides all necessary three-electron integrals + END_DOC + if(three_body_h_tc)then + PROVIDE three_e_3_idx_direct_bi_ort three_e_3_idx_cycle_1_bi_ort three_e_3_idx_cycle_2_bi_ort + PROVIDE three_e_3_idx_exch23_bi_ort three_e_3_idx_exch13_bi_ort three_e_3_idx_exch12_bi_ort + PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_cycle_2_bi_ort + PROVIDE three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort three_e_4_idx_exch12_bi_ort + endif +if(.not.double_normal_ord)then + PROVIDE three_e_5_idx_direct_bi_ort three_e_5_idx_cycle_1_bi_ort three_e_5_idx_cycle_2_bi_ort + PROVIDE three_e_5_idx_exch23_bi_ort three_e_5_idx_exch13_bi_ort three_e_5_idx_exch12_bi_ort +else + PROVIDE normal_two_body_bi_orth +endif +end + +subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree) + + BEGIN_DOC + ! diagonal element of htilde ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + double precision, intent(out) :: hthree + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm + integer(bit_kind) :: key_i_core(Nint,2) + double precision :: direct_int, exchange_int + double precision :: sym_3_e_int_from_6_idx_tensor + double precision :: three_e_diag_parrallel_spin + + if(core_tc_op)then + do i = 1, Nint + key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) + key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core,occ,Ne,Nint) + else + call bitstring_to_list_ab(key_i,occ,Ne,Nint) + endif + hthree = 0.d0 + + if(Ne(1)+Ne(2).ge.3)then +!! ! alpha/alpha/beta three-body + do i = 1, Ne(1) + ii = occ(i,1) + do j = i+1, Ne(1) + jj = occ(j,1) + do m = 1, Ne(2) + mm = occ(m,2) +! direct_int = three_body_ints_bi_ort(mm,jj,ii,mm,jj,ii) USES THE 6-IDX TENSOR +! exchange_int = three_body_ints_bi_ort(mm,jj,ii,mm,ii,jj) USES THE 6-IDX TENSOR + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) ! USES 3-IDX TENSOR + hthree += direct_int - exchange_int + enddo + enddo + enddo + + ! beta/beta/alpha three-body + do i = 1, Ne(2) + ii = occ(i,2) + do j = i+1, Ne(2) + jj = occ(j,2) + do m = 1, Ne(1) + mm = occ(m,1) + direct_int = three_e_3_idx_direct_bi_ort(mm,jj,ii) + exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,ii) + hthree += direct_int - exchange_int + enddo + enddo + enddo + + ! alpha/alpha/alpha three-body + do i = 1, Ne(1) + ii = occ(i,1) ! 1 + do j = i+1, Ne(1) + jj = occ(j,1) ! 2 + do m = j+1, Ne(1) + mm = occ(m,1) ! 3 +! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR + hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS + enddo + enddo + enddo + + ! beta/beta/beta three-body + do i = 1, Ne(2) + ii = occ(i,2) ! 1 + do j = i+1, Ne(2) + jj = occ(j,2) ! 2 + do m = j+1, Ne(2) + mm = occ(m,2) ! 3 +! ref = sym_3_e_int_from_6_idx_tensor(mm,jj,ii,mm,jj,ii) USES THE 6 IDX TENSOR + hthree += three_e_diag_parrallel_spin(mm,jj,ii) ! USES ONLY 3-IDX TENSORS + enddo + enddo + enddo + endif + +end + + +subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + + BEGIN_DOC + ! for single excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + double precision, intent(out) :: hthree + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2),i,j,ii,jj,ispin,jspin,k,kk + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: direct_int,phase,exchange_int,three_e_single_parrallel_spin + double precision :: sym_3_e_int_from_6_idx_tensor + integer :: other_spin(2) + integer(bit_kind) :: key_j_core(Nint,2),key_i_core(Nint,2) + + other_spin(1) = 2 + other_spin(2) = 1 + + + hthree = 0.d0 + call get_excitation_degree(key_i,key_j,degree,Nint) + if(degree.ne.1)then + return + endif + if(core_tc_op)then + do i = 1, Nint + key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) + key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) + key_j_core(i,1) = xor(key_j(i,1),core_bitmask(i,1)) + key_j_core(i,2) = xor(key_j(i,2),core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) + else + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + endif + + call get_single_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2) + + ! alpha/alpha/beta three-body +! print*,'IN SLAT RULES' + if(Ne(1)+Ne(2).ge.3)then + ! hole of spin s1 :: contribution from purely other spin + ispin = other_spin(s1) ! ispin is the other spin than s1 + do i = 1, Ne(ispin) ! i is the orbitals of the other spin than s1 + ii = occ(i,ispin) + do j = i+1, Ne(ispin) ! j has the same spin than s1 + jj = occ(j,ispin) + ! is == ispin in ::: s1 is is s1 is is s1 is is s1 is is + ! < h1 j i | p1 j i > - < h1 j i | p1 i j > + ! +! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor +! exchange_int = three_body_ints_bi_ort(jj,ii,p1,ii,jj,h1) ! USES THE 6-IDX tensor + direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1) + exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1) + hthree += direct_int - exchange_int + enddo + enddo + + ! hole of spin s1 :: contribution from mixed other spin / same spin + do i = 1, Ne(ispin) ! other spin + ii = occ(i,ispin) ! other spin + do j = 1, Ne(s1) ! same spin + jj = occ(j,s1) ! same spin +! direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor +! exchange_int = three_body_ints_bi_ort(jj,ii,p1,h1,ii,jj) ! exchange the spin s1 :: 6-IDX tensor + direct_int = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1) + exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1) + ! < h1 j i | p1 j i > - < h1 j i | j p1 i > + hthree += direct_int - exchange_int +! print*,'h1,p1,ii,jj = ',h1,p1,ii,jj +! print*,direct_int, exchange_int + enddo + enddo +! + ! hole of spin s1 :: PURE SAME SPIN CONTRIBUTIONS !!! + do i = 1, Ne(s1) + ii = occ(i,s1) + do j = i+1, Ne(s1) + jj = occ(j,s1) +! ref = sym_3_e_int_from_6_idx_tensor(jj,ii,p1,jj,ii,h1) + hthree += three_e_single_parrallel_spin(jj,ii,p1,h1) ! USES THE 4-IDX TENSOR + enddo + enddo + endif + hthree *= phase + +end + +subroutine double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree) + + BEGIN_DOC + ! for double excitation ONLY FOR THREE-BODY TERMS WITH BI ORTHONORMAL ORBITALS + !! + !! WARNING !! + ! + ! Non hermitian !! + END_DOC + + use bitmasks + + implicit none + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_j(Nint,2),key_i(Nint,2) + double precision, intent(out) :: hthree + integer :: occ(Nint*bit_kind_size,2) + integer :: Ne(2),i,j,ii,jj,ispin,jspin,m,mm + integer :: degree,exc(0:2,2,2) + integer :: h1, p1, h2, p2, s1, s2 + double precision :: phase + integer :: other_spin(2) + integer(bit_kind) :: key_i_core(Nint,2) + double precision :: direct_int,exchange_int,sym_3_e_int_from_6_idx_tensor + double precision :: three_e_double_parrallel_spin + + other_spin(1) = 2 + other_spin(2) = 1 + + call get_excitation_degree(key_i, key_j, degree, Nint) + + hthree = 0.d0 + + if(degree.ne.2)then + return + endif + + if(core_tc_op)then + do i = 1, Nint + key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1)) + key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2)) + enddo + call bitstring_to_list_ab(key_i_core, occ, Ne, Nint) + else + call bitstring_to_list_ab(key_i, occ, Ne, Nint) + endif + call get_double_excitation(key_i, key_j, exc, phase, Nint) + call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) + + + if(Ne(1)+Ne(2).ge.3)then + if(s1==s2)then ! same spin excitation + ispin = other_spin(s1) +! print*,'htilde ij' + do m = 1, Ne(ispin) ! direct(other_spin) - exchange(s1) + mm = occ(m,ispin) +!! direct_int = three_body_ints_bi_ort(mm,p2,p1,mm,h2,h1) +!! exchange_int = three_body_ints_bi_ort(mm,p2,p1,mm,h1,h2) + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch12_bi_ort(mm,p2,h2,p1,h1) +! print*,direct_int,exchange_int + hthree += direct_int - exchange_int + enddo + do m = 1, Ne(s1) ! pure contribution from s1 + mm = occ(m,s1) + hthree += three_e_double_parrallel_spin(mm,p2,h2,p1,h1) + enddo + else ! different spin excitation + do m = 1, Ne(s1) + mm = occ(m,s1) ! + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1) + hthree += direct_int - exchange_int + enddo + do m = 1, Ne(s2) + mm = occ(m,s2) ! + direct_int = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) + exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1) + hthree += direct_int - exchange_int + enddo + endif + endif + hthree *= phase + end diff --git a/src/tc_bi_ortho/symmetrized_3_e_int.irp.f b/src/tc_bi_ortho/symmetrized_3_e_int.irp.f new file mode 100644 index 00000000..e4f7ca93 --- /dev/null +++ b/src/tc_bi_ortho/symmetrized_3_e_int.irp.f @@ -0,0 +1,111 @@ +subroutine give_all_perm_for_three_e(n,l,k,m,j,i,idx_list,phase) + implicit none + BEGIN_DOC + ! returns all the list of permutting indices for the antimmetrization of + ! + ! (k^dagger l^dagger n^dagger m j i) when all indices have the same spins + ! + ! idx_list(:,i) == list of the 6 indices corresponding the permutation "i" + ! + ! phase(i) == phase of the permutation "i" + ! + ! there are in total 6 permutations with different indices + END_DOC + integer, intent(in) :: n,l,k,m,j,i + integer, intent(out) :: idx_list(6,6) + double precision :: phase(6) + integer :: list(6) + !!! CYCLIC PERMUTATIONS + phase(1:3) = 1.d0 + !!! IDENTITY PERMUTATION + list = (/n,l,k,m,j,i/) + idx_list(:,1) = list(:) + !!! FIRST CYCLIC PERMUTATION + list = (/n,l,k,j,i,m/) + idx_list(:,2) = list(:) + !!! FIRST CYCLIC PERMUTATION + list = (/n,l,k,i,m,j/) + idx_list(:,3) = list(:) + + !!! NON CYCLIC PERMUTATIONS + phase(1:3) = -1.d0 + !!! PARTICLE 1 is FIXED + list = (/n,l,k,j,m,i/) + idx_list(:,4) = list(:) + !!! PARTICLE 2 is FIXED + list = (/n,l,k,i,j,m/) + idx_list(:,5) = list(:) + !!! PARTICLE 3 is FIXED + list = (/n,l,k,m,i,j/) + idx_list(:,6) = list(:) + +end + +double precision function sym_3_e_int_from_6_idx_tensor(n,l,k,m,j,i) + implicit none + BEGIN_DOC + ! returns all good combinations of permutations of integrals with the good signs + ! + ! for a given (k^dagger l^dagger n^dagger m j i) when all indices have the same spins + END_DOC + integer, intent(in) :: n,l,k,m,j,i + sym_3_e_int_from_6_idx_tensor = three_body_ints_bi_ort(n,l,k,m,j,i) & ! direct + + three_body_ints_bi_ort(n,l,k,j,i,m) & ! 1st cyclic permutation + + three_body_ints_bi_ort(n,l,k,i,m,j) & ! 2nd cyclic permutation + - three_body_ints_bi_ort(n,l,k,j,m,i) & ! elec 1 is kept fixed + - three_body_ints_bi_ort(n,l,k,i,j,m) & ! elec 2 is kept fixed + - three_body_ints_bi_ort(n,l,k,m,i,j) ! elec 3 is kept fixed + +end + +double precision function direct_sym_3_e_int(n,l,k,m,j,i) + implicit none + BEGIN_DOC + ! returns all good combinations of permutations of integrals with the good signs + ! + ! for a given (k^dagger l^dagger n^dagger m j i) when all indices have the same spins + END_DOC + integer, intent(in) :: n,l,k,m,j,i + double precision :: integral + direct_sym_3_e_int = 0.d0 + call give_integrals_3_body_bi_ort(n,l,k,m,j,i,integral) ! direct + direct_sym_3_e_int += integral + call give_integrals_3_body_bi_ort(n,l,k,j,i,m,integral) ! 1st cyclic permutation + direct_sym_3_e_int += integral + call give_integrals_3_body_bi_ort(n,l,k,i,m,j,integral) ! 2nd cyclic permutation + direct_sym_3_e_int += integral + call give_integrals_3_body_bi_ort(n,l,k,j,m,i,integral) ! elec 1 is kept fixed + direct_sym_3_e_int += -integral + call give_integrals_3_body_bi_ort(n,l,k,i,j,m,integral) ! elec 2 is kept fixed + direct_sym_3_e_int += -integral + call give_integrals_3_body_bi_ort(n,l,k,m,i,j,integral) ! elec 3 is kept fixed + direct_sym_3_e_int += -integral + +end + +double precision function three_e_diag_parrallel_spin(m,j,i) + implicit none + integer, intent(in) :: i,j,m + three_e_diag_parrallel_spin = three_e_3_idx_direct_bi_ort(m,j,i) ! direct + three_e_diag_parrallel_spin += three_e_3_idx_cycle_1_bi_ort(m,j,i) + three_e_3_idx_cycle_2_bi_ort(m,j,i) & ! two cyclic permutations + - three_e_3_idx_exch23_bi_ort(m,j,i) - three_e_3_idx_exch13_bi_ort(m,j,i) & ! two first exchange + - three_e_3_idx_exch12_bi_ort(m,j,i) ! last exchange +end + +double precision function three_e_single_parrallel_spin(m,j,k,i) + implicit none + integer, intent(in) :: i,k,j,m + three_e_single_parrallel_spin = three_e_4_idx_direct_bi_ort(m,j,k,i) ! direct + three_e_single_parrallel_spin += three_e_4_idx_cycle_1_bi_ort(m,j,k,i) + three_e_4_idx_cycle_2_bi_ort(m,j,k,i) & ! two cyclic permutations + - three_e_4_idx_exch23_bi_ort(m,j,k,i) - three_e_4_idx_exch13_bi_ort(m,j,k,i) & ! two first exchange + - three_e_4_idx_exch12_bi_ort(m,j,k,i) ! last exchange +end + +double precision function three_e_double_parrallel_spin(m,l,j,k,i) + implicit none + integer, intent(in) :: i,k,j,m,l + three_e_double_parrallel_spin = three_e_5_idx_direct_bi_ort(m,l,j,k,i) ! direct + three_e_double_parrallel_spin += three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) + three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) & ! two cyclic permutations + - three_e_5_idx_exch23_bi_ort(m,l,j,k,i) - three_e_5_idx_exch13_bi_ort(m,l,j,k,i) & ! two first exchange + - three_e_5_idx_exch12_bi_ort(m,l,j,k,i) ! last exchange +end diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f new file mode 100644 index 00000000..cfa24f3b --- /dev/null +++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f @@ -0,0 +1,61 @@ +program tc_bi_ortho + implicit none + BEGIN_DOC +! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call routine_diag +! call test +end + +subroutine test + implicit none + integer :: i,j + double precision :: hmono,htwoe,hthree,htot + use bitmasks + + print*,'test' +! call htilde_mu_mat_bi_ortho(psi_det(1,1,1), psi_det(1,1,2), N_int, hmono, htwoe, hthree, htot) + call double_htilde_mu_mat_bi_ortho(N_int,psi_det(1,1,1), psi_det(1,1,2), hmono, htwoe, htot) + print*,hmono, htwoe, htot + +end + +subroutine routine_diag + implicit none +! provide eigval_right_tc_bi_orth + provide overlap_bi_ortho +! provide htilde_matrix_elmt_bi_ortho + integer ::i,j + print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) + print*,'e_tc_left_right = ',e_tc_left_right + print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00 + print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth + print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single + print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double + print*,'***' + print*,'e_corr_bi_orth = ',e_corr_bi_orth + print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj + print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth + print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth + print*,'Left/right eigenvectors' + do i = 1,N_det + write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1) + enddo + do j=1,N_states + do i=1,N_det + psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j) + psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j) + enddo + enddo + SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho + call save_tc_bi_ortho_wavefunction +! call routine_save_left_right_bi_ortho +end + diff --git a/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f new file mode 100644 index 00000000..28f122ee --- /dev/null +++ b/src/tc_bi_ortho/tc_bi_ortho_prop.irp.f @@ -0,0 +1,24 @@ +program tc_bi_ortho_prop + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid +! call routine_diag + call test +end + +subroutine test + implicit none + integer :: i + print*,'TC Dipole components' + do i= 1, 3 + print*,tc_bi_ortho_dipole(i,1) + enddo +end diff --git a/src/tc_bi_ortho/tc_cisd_sc2.irp.f b/src/tc_bi_ortho/tc_cisd_sc2.irp.f new file mode 100644 index 00000000..0fb9f524 --- /dev/null +++ b/src/tc_bi_ortho/tc_cisd_sc2.irp.f @@ -0,0 +1,24 @@ +program tc_bi_ortho + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + call test +end + +subroutine test + implicit none +! double precision, allocatable :: dressing_dets(:),e_corr_dets(:) +! allocate(dressing_dets(N_det),e_corr_dets(N_det)) +! e_corr_dets = 0.d0 +! call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets) + provide eigval_tc_cisd_sc2_bi_ortho +end diff --git a/src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f b/src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f new file mode 100644 index 00000000..406ee9e3 --- /dev/null +++ b/src/tc_bi_ortho/tc_cisd_sc2_utils.irp.f @@ -0,0 +1,110 @@ + BEGIN_PROVIDER [ double precision, reigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)] +&BEGIN_PROVIDER [ double precision, leigvec_tc_cisd_sc2_bi_ortho, (N_det,N_states)] +&BEGIN_PROVIDER [ double precision, eigval_tc_cisd_sc2_bi_ortho, (N_states)] + implicit none + integer :: it,n_real,degree,i + double precision :: e_before, e_current,thr, hmono,htwoe,hthree + double precision, allocatable :: e_corr_dets(:),h0j(:), h_sc2(:,:), dressing_dets(:) + double precision, allocatable :: leigvec_tc_bi_orth_tmp(:,:),reigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:) + allocate(leigvec_tc_bi_orth_tmp(N_det,N_det),reigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det)) + allocate(e_corr_dets(N_det),h0j(N_det),h_sc2(N_det,N_det),dressing_dets(N_det)) + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree == 1 .or. degree == 2)then + call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,h0j(i)) + endif + enddo + do i = 1, N_det + e_corr_dets(i) = reigvec_tc_bi_orth(i,1) * h0j(i)/reigvec_tc_bi_orth(1,1) + enddo + print*,'Starting from ',eigval_right_tc_bi_orth(1) + + e_before = 0.d0 + e_current = 10.d0 + thr = 1.d-5 + it = 0 + dressing_dets = 0.d0 + do while (dabs(E_before-E_current).gt.thr) + it += 1 + E_before = E_current + h_sc2 = htilde_matrix_elmt_bi_ortho + call get_cisd_sc2_dressing(psi_det,e_corr_dets,N_det,dressing_dets) + do i = 1, N_det + print*,'dressing_dets(i) = ',dressing_dets(i) + h_sc2(i,i) += dressing_dets(i) + enddo + call non_hrmt_real_diag(N_det,h_sc2,& + leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& + n_real,eigval_right_tmp) + do i = 1, N_det + e_corr_dets(i) = reigvec_tc_bi_orth_tmp(i,1) * h0j(i)/reigvec_tc_bi_orth_tmp(1,1) + enddo + E_current = eigval_right_tmp(1) + print*,'it, E(SC)^2 = ',it,E_current + enddo + eigval_tc_cisd_sc2_bi_ortho(1:N_states) = eigval_right_tmp(1:N_states) + reigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = reigvec_tc_bi_orth_tmp(1:N_det,1:N_states) + leigvec_tc_cisd_sc2_bi_ortho(1:N_det,1:N_states) = leigvec_tc_bi_orth_tmp(1:N_det,1:N_states) + +END_PROVIDER + +subroutine get_cisd_sc2_dressing(dets,e_corr_dets,ndet,dressing_dets) + implicit none + use bitmasks + integer, intent(in) :: ndet + integer(bit_kind), intent(in) :: dets(N_int,2,ndet) + double precision, intent(in) :: e_corr_dets(ndet) + double precision, intent(out) :: dressing_dets(ndet) + integer, allocatable :: degree(:),hole(:,:),part(:,:),spin(:,:) + integer(bit_kind), allocatable :: hole_part(:,:,:) + integer :: i,j,k, exc(0:2,2,2),h1,p1,h2,p2,s1,s2 + integer(bit_kind) :: xorvec(2,N_int) + + double precision :: phase + dressing_dets = 0.d0 + allocate(degree(ndet),hole(2,ndet),part(2,ndet), spin(2,ndet),hole_part(N_int,2,ndet)) + do i = 2, ndet + call get_excitation_degree(HF_bitmask,dets(1,1,i),degree(i),N_int) + do j = 1, N_int + hole_part(j,1,i) = xor( HF_bitmask(j,1), dets(j,1,i)) + hole_part(j,2,i) = xor( HF_bitmask(j,2), dets(j,2,i)) + enddo + if(degree(i) == 1)then + call get_single_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int) + else if(degree(i) == 2)then + call get_double_excitation(HF_bitmask,psi_det(1,1,i),exc,phase,N_int) + endif + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + hole(1,i) = h1 + hole(2,i) = h2 + part(1,i) = p1 + part(2,i) = p2 + spin(1,i) = s1 + spin(2,i) = s2 + enddo + + integer :: same + if(elec_alpha_num+elec_beta_num<3)return + do i = 2, ndet + do j = i+1, ndet + same = 0 + if(degree(i) == degree(j) .and. degree(i)==1)cycle + do k = 1, N_int + xorvec(k,1) = iand(hole_part(k,1,i),hole_part(k,1,j)) + xorvec(k,2) = iand(hole_part(k,2,i),hole_part(k,2,j)) + same += popcnt(xorvec(k,1)) + popcnt(xorvec(k,2)) + enddo +! print*,'i,j',i,j +! call debug_det(dets(1,1,i),N_int) +! call debug_det(hole_part(1,1,i),N_int) +! call debug_det(dets(1,1,j),N_int) +! call debug_det(hole_part(1,1,j),N_int) +! print*,'same = ',same + if(same.eq.0)then + dressing_dets(i) += e_corr_dets(j) + dressing_dets(j) += e_corr_dets(i) + endif + enddo + enddo + +end diff --git a/src/tc_bi_ortho/tc_h_eigvectors.irp.f b/src/tc_bi_ortho/tc_h_eigvectors.irp.f new file mode 100644 index 00000000..651088be --- /dev/null +++ b/src/tc_bi_ortho/tc_h_eigvectors.irp.f @@ -0,0 +1,179 @@ + use bitmasks + + BEGIN_PROVIDER [ integer, index_HF_psi_det] + implicit none + integer :: i,degree + do i = 1, N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int) + if(degree == 0)then + index_HF_psi_det = i + exit + endif + enddo + END_PROVIDER + + + + BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)] +&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)] +&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)] +&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)] +&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ] + + BEGIN_DOC + ! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis + END_DOC + + implicit none + integer :: i, idx_dress, j, istate + logical :: converged, dagger + integer :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l + double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:) + + PROVIDE N_det N_int + + if(n_det.le.N_det_max_full)then + allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det)) + call non_hrmt_real_diag(N_det,htilde_matrix_elmt_bi_ortho,& + leigvec_tc_bi_orth_tmp,reigvec_tc_bi_orth_tmp,& + n_real_tc_bi_orth_eigval_right,eigval_right_tmp) + double precision, allocatable :: coef_hf_r(:),coef_hf_l(:) + integer, allocatable :: iorder(:) + allocate(coef_hf_r(N_det),coef_hf_l(N_det),iorder(N_det)) + do i = 1,N_det + iorder(i) = i + coef_hf_r(i) = -dabs(reigvec_tc_bi_orth_tmp(index_HF_psi_det,i)) + enddo + call dsort(coef_hf_r,iorder,N_det) + igood_r = iorder(1) + print*,'igood_r, coef_hf_r = ',igood_r,coef_hf_r(1) + do i = 1,N_det + iorder(i) = i + coef_hf_l(i) = -dabs(leigvec_tc_bi_orth_tmp(index_HF_psi_det,i)) + enddo + call dsort(coef_hf_l,iorder,N_det) + igood_l = iorder(1) + print*,'igood_l, coef_hf_l = ',igood_l,coef_hf_l(1) + + if(igood_r.ne.igood_l.and.igood_r.ne.1)then + print *,'' + print *,'Warning, the left and right eigenvectors are "not the same" ' + print *,'Warning, the ground state is not dominated by HF...' + print *,'State with largest RIGHT coefficient of HF ',igood_r + print *,'coef of HF in RIGHT eigenvector = ',reigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_r) + print *,'State with largest LEFT coefficient of HF ',igood_l + print *,'coef of HF in LEFT eigenvector = ',leigvec_tc_bi_orth_tmp(index_HF_psi_det,igood_l) + endif + if(state_following_tc)then + print *,'Following the states with the largest coef on HF' + print *,'igood_r,igood_l',igood_r,igood_l + i= igood_r + eigval_right_tc_bi_orth(1) = eigval_right_tmp(i) + do j = 1, N_det + reigvec_tc_bi_orth(j,1) = reigvec_tc_bi_orth_tmp(j,i) +! print*,reigvec_tc_bi_orth(j,1) + enddo + i= igood_l + eigval_left_tc_bi_orth(1) = eigval_right_tmp(i) + do j = 1, N_det + leigvec_tc_bi_orth(j,1) = leigvec_tc_bi_orth_tmp(j,i) + enddo + else + do i = 1, N_states + eigval_right_tc_bi_orth(i) = eigval_right_tmp(i) + eigval_left_tc_bi_orth(i) = eigval_right_tmp(i) + do j = 1, N_det + reigvec_tc_bi_orth(j,i) = reigvec_tc_bi_orth_tmp(j,i) + leigvec_tc_bi_orth(j,i) = leigvec_tc_bi_orth_tmp(j,i) + enddo + enddo + endif + else + double precision, allocatable :: H_jj(:),vec_tmp(:,:) + external htc_bi_ortho_calc_tdav + external htcdag_bi_ortho_calc_tdav + allocate(H_jj(N_det),vec_tmp(N_det,n_states_diag)) + do i = 1, N_det + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, H_jj(i)) + enddo + !!!! Preparing the left-eigenvector + print*,'Computing the left-eigenvector ' + vec_tmp = 0.d0 + do istate = 1, N_states + vec_tmp(:,istate) = psi_l_coef_bi_ortho(:,istate) + enddo + do istate = N_states+1, n_states_diag + vec_tmp(istate,istate) = 1.d0 + enddo + call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_left_tc_bi_orth, N_det, n_states, n_states_diag, converged, htcdag_bi_ortho_calc_tdav) + do istate = 1, N_states + leigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate) + enddo + + print*,'Computing the right-eigenvector ' + !!!! Preparing the right-eigenvector + vec_tmp = 0.d0 + do istate = 1, N_states + vec_tmp(:,istate) = psi_r_coef_bi_ortho(:,istate) + enddo + do istate = N_states+1, n_states_diag + vec_tmp(istate,istate) = 1.d0 + enddo + call davidson_general_ext_rout_nonsym_b1space(vec_tmp, H_jj, eigval_right_tc_bi_orth, N_det, n_states, n_states_diag, converged, htc_bi_ortho_calc_tdav) + do istate = 1, N_states + reigvec_tc_bi_orth(:,istate) = vec_tmp(:,istate) + enddo + + deallocate(H_jj) + endif + call bi_normalize(leigvec_tc_bi_orth,reigvec_tc_bi_orth,N_det,N_det,N_states) + print*,'leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) = ',leigvec_tc_bi_orth(1,1),reigvec_tc_bi_orth(1,1) + norm_ground_left_right_bi_orth = 0.d0 + do j = 1, N_det + norm_ground_left_right_bi_orth += leigvec_tc_bi_orth(j,1) * reigvec_tc_bi_orth(j,1) + enddo + print*,'norm l/r = ',norm_ground_left_right_bi_orth + +END_PROVIDER + + + +subroutine bi_normalize(u_l,u_r,n,ld,nstates) + !!!! Normalization of the scalar product of the left/right eigenvectors + double precision, intent(inout) :: u_l(ld,nstates), u_r(ld,nstates) + integer, intent(in) :: n,ld,nstates + integer :: i + double precision :: accu, tmp + do i = 1, nstates + !!!! Normalization of right eigenvectors |Phi> + accu = 0.d0 + do j = 1, n + accu += u_r(j,i) * u_r(j,i) + enddo + accu = 1.d0/dsqrt(accu) + print*,'accu_r = ',accu + do j = 1, n + u_r(j,i) *= accu + enddo + tmp = u_r(1,i) / dabs(u_r(1,i)) + do j = 1, n + u_r(j,i) *= tmp + enddo + !!!! Adaptation of the norm of the left eigenvector such that = 1 + accu = 0.d0 + do j = 1, n + accu += u_l(j,i) * u_r(j,i) +! print*,j, u_l(j,i) , u_r(j,i) + enddo + if(accu.gt.0.d0)then + accu = 1.d0/dsqrt(accu) + else + accu = 1.d0/dsqrt(-accu) + endif + tmp = (u_l(1,i) * u_r(1,i) )/dabs(u_l(1,i) * u_r(1,i)) + do j = 1, n + u_l(j,i) *= accu * tmp + u_r(j,i) *= accu + enddo + enddo +end diff --git a/src/tc_bi_ortho/tc_hmat.irp.f b/src/tc_bi_ortho/tc_hmat.irp.f new file mode 100644 index 00000000..bf1388e5 --- /dev/null +++ b/src/tc_bi_ortho/tc_hmat.irp.f @@ -0,0 +1,41 @@ + + BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)] + + BEGIN_DOC + ! htilde_matrix_elmt_bi_ortho(j,i) = + ! + ! WARNING !!!!!!!!! IT IS NOT HERMITIAN !!!!!!!!! + END_DOC + + implicit none + integer :: i, j + double precision :: hmono,htwoe,hthree,htot + + PROVIDE N_int + !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hmono, htwoe, hthree, htot) & + !$OMP SHARED (N_det, psi_det, N_int,htilde_matrix_elmt_bi_ortho) + do i = 1, N_det + do j = 1, N_det + ! < J | Htilde | I > + call htilde_mu_mat_bi_ortho(psi_det(1,1,j), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot) + htilde_matrix_elmt_bi_ortho(j,i) = htot + enddo + enddo + !$OMP END PARALLEL DO +! print*,'htilde_matrix_elmt_bi_ortho = ' +! do i = 1, min(100,N_det) +! write(*,'(100(F16.10,X))')htilde_matrix_elmt_bi_ortho(1:min(100,N_det),i) +! enddo + + +END_PROVIDER + + BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)] + implicit none + integer ::i,j + do i = 1, N_det + do j = 1, N_det + htilde_matrix_elmt_bi_ortho_tranp(j,i) = htilde_matrix_elmt_bi_ortho(i,j) + enddo + enddo +END_PROVIDER diff --git a/src/tc_bi_ortho/tc_prop.irp.f b/src/tc_bi_ortho/tc_prop.irp.f new file mode 100644 index 00000000..6db8e80e --- /dev/null +++ b/src/tc_bi_ortho/tc_prop.irp.f @@ -0,0 +1,268 @@ + +BEGIN_PROVIDER [ double precision, tc_transition_matrix, (mo_num, mo_num,N_states,N_states) ] + implicit none + BEGIN_DOC + ! tc_transition_matrix(p,h,istate,jstate) = + ! + ! where are the left/right eigenvectors on a bi-ortho basis + END_DOC + integer :: i,j,istate,jstate,m,n,p,h + double precision :: phase + integer, allocatable :: occ(:,:) + integer :: n_occ_ab(2),degree,exc(0:2,2,2) + allocate(occ(N_int*bit_kind_size,2)) + tc_transition_matrix = 0.d0 + do istate = 1, N_states + do jstate = 1, N_states + do i = 1, N_det + do j = 1, N_det + call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) + if(degree.gt.1)then + cycle + else if (degree == 0)then + call bitstring_to_list_ab(psi_det(1,1,i), occ, n_occ_ab, N_int) + do p = 1, n_occ_ab(1) ! browsing the alpha electrons + m = occ(p,1) + tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + enddo + do p = 1, n_occ_ab(2) ! browsing the beta electrons + m = occ(p,1) + tc_transition_matrix(m,m,istate,jstate)+= psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + enddo + else + call get_single_excitation(psi_det(1,1,j),psi_det(1,1,i),exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Single alpha + h = exc(1,1,1) ! hole in psi_det(1,1,j) + p = exc(1,2,1) ! particle in psi_det(1,1,j) + else + ! Single beta + h = exc(1,1,2) ! hole in psi_det(1,1,j) + p = exc(1,2,2) ! particle in psi_det(1,1,j) + endif + tc_transition_matrix(p,h,istate,jstate)+= phase * psi_l_coef_bi_ortho(i,istate) * psi_r_coef_bi_ortho(j,jstate) + endif + enddo + enddo + enddo + enddo + END_PROVIDER + + + BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_mo, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_mo, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, natorb_tc_eigval, (mo_num)] + implicit none + BEGIN_DOC + ! natorb_tc_reigvec_mo : RIGHT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) + ! natorb_tc_leigvec_mo : LEFT eigenvectors of the ground state transition matrix (equivalent of natural orbitals) + ! natorb_tc_eigval : eigenvalues of the ground state transition matrix (equivalent of the occupation numbers). WARNINING :: can be negative !! + END_DOC + double precision, allocatable :: dm_tmp(:,:) + integer :: i,j,k,n_real + allocate( dm_tmp(mo_num,mo_num)) + dm_tmp(:,:) = -tc_transition_matrix(:,:,1,1) + print*,'dm_tmp' + do i = 1, mo_num + write(*,'(100(F16.10,X))')-dm_tmp(:,i) + enddo +! call non_hrmt_diag_split_degen( mo_num, dm_tmp& + call non_hrmt_fock_mat( mo_num, dm_tmp& +! call non_hrmt_bieig( mo_num, dm_tmp& + , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo& + , n_real, natorb_tc_eigval ) + double precision :: accu + accu = 0.d0 + do i = 1, n_real + print*,'natorb_tc_eigval(i) = ',-natorb_tc_eigval(i) + accu += -natorb_tc_eigval(i) + enddo + print*,'accu = ',accu + dm_tmp = 0.d0 + do i = 1, n_real + accu = 0.d0 + do k = 1, mo_num + accu += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,i) + enddo + accu = 1.d0/dsqrt(dabs(accu)) + natorb_tc_reigvec_mo(:,i) *= accu + natorb_tc_leigvec_mo(:,i) *= accu + do j = 1, n_real + do k = 1, mo_num + dm_tmp(j,i) += natorb_tc_reigvec_mo(k,i) * natorb_tc_leigvec_mo(k,j) + enddo + enddo + enddo + double precision :: accu_d, accu_nd + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, mo_num + accu_d += dm_tmp(i,i) + ! write(*,'(100(F16.10,X))')dm_tmp(:,i) + do j = 1, mo_num + if(i==j)cycle + accu_nd += dabs(dm_tmp(j,i)) + enddo + enddo + print*,'Trace of the overlap between TC natural orbitals ',accu_d + print*,'L1 norm of extra diagonal elements of overlap matrix ',accu_nd + + + END_PROVIDER + + BEGIN_PROVIDER [ double precision, fock_diag_sorted_r_natorb, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, fock_diag_sorted_l_natorb, (mo_num, mo_num)] + &BEGIN_PROVIDER [ double precision, fock_diag_sorted_v_natorb, (mo_num)] + implicit none + integer ::i,j,k + print*,'Diagonal elements of the Fock matrix before ' + do i = 1, mo_num + write(*,*)i,Fock_matrix_tc_mo_tot(i,i) + enddo + double precision, allocatable :: fock_diag(:) + allocate(fock_diag(mo_num)) + fock_diag = 0.d0 + do i = 1, mo_num + fock_diag(i) = 0.d0 + do j = 1, mo_num + do k = 1, mo_num + fock_diag(i) += natorb_tc_leigvec_mo(k,i) * Fock_matrix_tc_mo_tot(k,j) * natorb_tc_reigvec_mo(j,i) + enddo + enddo + enddo + integer, allocatable :: iorder(:) + allocate(iorder(mo_num)) + do i = 1, mo_num + iorder(i) = i + enddo + call dsort(fock_diag,iorder,mo_num) + print*,'Diagonal elements of the Fock matrix after ' + do i = 1, mo_num + write(*,*)i,fock_diag(i) + enddo + do i = 1, mo_num + fock_diag_sorted_v_natorb(i) = natorb_tc_eigval(iorder(i)) + do j = 1, mo_num + fock_diag_sorted_r_natorb(j,i) = natorb_tc_reigvec_mo(j,iorder(i)) + fock_diag_sorted_l_natorb(j,i) = natorb_tc_leigvec_mo(j,iorder(i)) + enddo + enddo + + END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, natorb_tc_reigvec_ao, (ao_num, mo_num)] + &BEGIN_PROVIDER [ double precision, natorb_tc_leigvec_ao, (ao_num, mo_num)] + &BEGIN_PROVIDER [ double precision, overlap_natorb_tc_eigvec_ao, (mo_num, mo_num) ] + + BEGIN_DOC + ! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP + ! + ! THE OVERLAP SHOULD BE THE SAME AS overlap_natorb_tc_eigvec_mo + END_DOC + + implicit none + integer :: i, j, k, q, p + double precision :: accu, accu_d + double precision, allocatable :: tmp(:,:) + + + ! ! MO_R x R + call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & + , mo_r_coef, size(mo_r_coef, 1) & + , fock_diag_sorted_r_natorb, size(fock_diag_sorted_r_natorb, 1) & + , 0.d0, natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) ) + ! + ! MO_L x L + call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & + , mo_l_coef, size(mo_l_coef, 1) & + , fock_diag_sorted_l_natorb, size(fock_diag_sorted_l_natorb, 1) & + , 0.d0, natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1) ) + + + allocate( tmp(mo_num,ao_num) ) + + ! tmp <-- L.T x S_ao + call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 & + , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + ! S <-- tmp x R + call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 & + , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) & + , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) ) + + deallocate( tmp ) + + ! --- + double precision :: norm + do i = 1, mo_num + norm = 1.d0/dsqrt(dabs(overlap_natorb_tc_eigvec_ao(i,i))) + do j = 1, mo_num + natorb_tc_reigvec_ao(j,i) *= norm + natorb_tc_leigvec_ao(j,i) *= norm + enddo + enddo + + allocate( tmp(mo_num,ao_num) ) + + ! tmp <-- L.T x S_ao + call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 & + , natorb_tc_leigvec_ao, size(natorb_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + ! S <-- tmp x R + call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 & + , tmp, size(tmp, 1), natorb_tc_reigvec_ao, size(natorb_tc_reigvec_ao, 1) & + , 0.d0, overlap_natorb_tc_eigvec_ao, size(overlap_natorb_tc_eigvec_ao, 1) ) + + + + deallocate( tmp ) + + accu_d = 0.d0 + accu = 0.d0 + do i = 1, mo_num + accu_d += overlap_natorb_tc_eigvec_ao(i,i) + do j = 1, mo_num + if(i==j)cycle + accu += dabs(overlap_natorb_tc_eigvec_ao(j,i)) + enddo + enddo + print*,'Trace of the overlap_natorb_tc_eigvec_ao = ',accu_d + print*,'mo_num = ',mo_num + print*,'L1 norm of extra diagonal elements of overlap matrix ',accu + accu = accu / dble(mo_num**2) + + END_PROVIDER + + BEGIN_PROVIDER [double precision, tc_bi_ortho_dipole, (3,N_states)] + implicit none + integer :: i,j,istate,m + double precision :: nuclei_part(3) + tc_bi_ortho_dipole = 0.d0 + do istate = 1, N_states + do i = 1, mo_num + do j = 1, mo_num + tc_bi_ortho_dipole(1,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_x(j,i) + tc_bi_ortho_dipole(2,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_y(j,i) + tc_bi_ortho_dipole(3,istate) += -(tc_transition_matrix(j,i,istate,istate)) * mo_bi_orth_bipole_z(j,i) + enddo + enddo + enddo + + nuclei_part = 0.d0 + do m = 1, 3 + do i = 1,nucl_num + nuclei_part(m) += nucl_charge(i) * nucl_coord(i,m) + enddo + enddo +! + do istate = 1, N_states + do m = 1, 3 + tc_bi_ortho_dipole(m,istate) += nuclei_part(m) + enddo + enddo + END_PROVIDER + diff --git a/src/tc_bi_ortho/test_normal_order.irp.f b/src/tc_bi_ortho/test_normal_order.irp.f new file mode 100644 index 00000000..8bdc57ee --- /dev/null +++ b/src/tc_bi_ortho/test_normal_order.irp.f @@ -0,0 +1,73 @@ +program test_normal_order + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + call test +end + +subroutine test + implicit none + use bitmasks ! you need to include the bitmasks_module.f90 features + integer :: h1,h2,p1,p2,s1,s2,i_ok,degree,Ne(2) + integer :: exc(0:2,2,2) + integer(bit_kind), allocatable :: det_i(:,:) + double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal + integer, allocatable :: occ(:,:) + allocate( occ(N_int*bit_kind_size,2) ) + call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int) + allocate(det_i(N_int,2)) + s1 = 1 + s2 = 2 + accu = 0.d0 + do h1 = 1, elec_beta_num + do p1 = elec_beta_num+1, mo_num + do h2 = 1, elec_beta_num + do p2 = elec_beta_num+1, mo_num + det_i = ref_bitmask + call do_single_excitation(det_i,h1,p1,s1,i_ok) + call do_single_excitation(det_i,h2,p2,s2,i_ok) + call htilde_mu_mat_bi_ortho(det_i,HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call get_excitation_degree(ref_bitmask,det_i,degree,N_int) + call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) + hthree *= phase + normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1) + accu += dabs(hthree-normal) + enddo + enddo + enddo + enddo + print*,'accu opposite spin = ',accu + + s1 = 2 + s2 = 2 + accu = 0.d0 + do h1 = 1, elec_beta_num + do p1 = elec_beta_num+1, mo_num + do h2 = h1+1, elec_beta_num + do p2 = elec_beta_num+1, mo_num + det_i = ref_bitmask + call do_single_excitation(det_i,h1,p1,s1,i_ok) + call do_single_excitation(det_i,h2,p2,s2,i_ok) + if(i_ok.ne.1)cycle + call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + call get_excitation_degree(ref_bitmask,det_i,degree,N_int) + call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) + hthree *= phase + normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1) + accu += dabs(hthree-normal) + enddo + enddo + enddo + enddo + print*,'accu same spin = ',accu +end + + diff --git a/src/tc_bi_ortho/test_tc_bi_ortho.irp.f b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f new file mode 100644 index 00000000..2d71b6b2 --- /dev/null +++ b/src/tc_bi_ortho/test_tc_bi_ortho.irp.f @@ -0,0 +1,131 @@ +program tc_bi_ortho + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + ! call routine_2 + call test_rout +end + +subroutine test_rout + implicit none + integer :: i,j,ii,jj + use bitmasks ! you need to include the bitmasks_module.f90 features + integer(bit_kind), allocatable :: det_i(:,:) + allocate(det_i(N_int,2)) + det_i(:,:)= psi_det(:,:,1) + call debug_det(det_i,N_int) + integer, allocatable :: occ(:,:) + integer :: n_occ_ab(2) + allocate(occ(N_int*bit_kind_size,2)) + call bitstring_to_list_ab(det_i, occ, n_occ_ab, N_int) + double precision :: hmono, htwoe, htot + call diag_htilde_mu_mat_bi_ortho(N_int, det_i, hmono, htwoe, htot) + print*,'hmono, htwoe, htot' + print*, hmono, htwoe, htot + print*,'alpha electrons orbital occupancy' + do i = 1, n_occ_ab(1) ! browsing the alpha electrons + j = occ(i,1) + print*,j,mo_bi_ortho_tc_one_e(j,j) + enddo + print*,'beta electrons orbital occupancy' + do i = 1, n_occ_ab(2) ! browsing the beta electrons + j = occ(i,2) + print*,j,mo_bi_ortho_tc_one_e(j,j) + enddo + print*,'alpha beta' + do i = 1, n_occ_ab(1) + ii = occ(i,1) + do j = 1, n_occ_ab(2) + jj = occ(j,2) + print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii) + enddo + enddo + print*,'alpha alpha' + do i = 1, n_occ_ab(1) + ii = occ(i,1) + do j = 1, n_occ_ab(1) + jj = occ(j,1) + print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii) + enddo + enddo + + print*,'beta beta' + do i = 1, n_occ_ab(2) + ii = occ(i,2) + do j = 1, n_occ_ab(2) + jj = occ(j,2) + print*,ii,jj,mo_bi_ortho_tc_two_e(jj,ii,jj,ii), mo_bi_ortho_tc_two_e(ii,jj,jj,ii) + enddo + enddo + + +end + +subroutine routine_2 + implicit none + integer :: i + double precision :: bi_ortho_mo_ints + print*,'H matrix' + do i = 1, N_det + write(*,'(1000(F16.5,X))')htilde_matrix_elmt_bi_ortho(:,i) + enddo + i = 1 + double precision :: phase + integer :: degree,h1, p1, h2, p2, s1, s2, exc(0:2,2,2) + call get_excitation_degree(ref_bitmask, psi_det(1,1,i), degree, N_int) + if(degree==2)then + call get_double_excitation(ref_bitmask, psi_det(1,1,i), exc, phase, N_int) + call decode_exc(exc, 2, h1, p1, h2, p2, s1, s2) + print*,'h1,h2,p1,p2' + print*, h1,h2,p1,p2 + print*,mo_bi_ortho_tc_two_e(p1,p2,h1,h2),mo_bi_ortho_tc_two_e(h1,h2,p1,p2) + endif + + + print*,'coef' + do i = 1, ao_num + print*,i,mo_l_coef(i,8),mo_r_coef(i,8) + enddo +! print*,'mdlqfmlqgmqglj' +! print*,'mo_bi_ortho_tc_two_e()',mo_bi_ortho_tc_two_e(2,2,3,3) +! print*,'bi_ortho_mo_ints ',bi_ortho_mo_ints(2,2,3,3) + print*,'Overlap' + do i = 1, mo_num + write(*,'(100(F16.10,X))')overlap_bi_ortho(:,i) + enddo + +end + +subroutine routine + implicit none + double precision :: hmono,htwoe,hthree,htot + integer(bit_kind), allocatable :: key1(:,:) + integer(bit_kind), allocatable :: key2(:,:) + allocate(key1(N_int,2),key2(N_int,2)) + use bitmasks + key1 = ref_bitmask + call htilde_mu_mat_bi_ortho(key1,key1, N_int, hmono,htwoe,hthree,htot) + key2 = key1 + integer :: h,p,i_ok + h = 1 + p = 8 + call do_single_excitation(key2,h,p,1,i_ok) + call debug_det(key2,N_int) + call htilde_mu_mat_bi_ortho(key2,key1, N_int, hmono,htwoe,hthree,htot) +! print*,'fock_matrix_tc_mo_alpha(p,h) = ',fock_matrix_tc_mo_alpha(p,h) + print*,'htot = ',htot + print*,'hmono = ',hmono + print*,'htwoe = ',htwoe + double precision :: bi_ortho_mo_ints + print*,'bi_ortho_mo_ints(1,p,1,h)',bi_ortho_mo_ints(1,p,1,h) + +end diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f new file mode 100644 index 00000000..badc40b6 --- /dev/null +++ b/src/tc_bi_ortho/test_tc_fock.irp.f @@ -0,0 +1,169 @@ +program test_tc_fock + implicit none + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + print *, 'Hello world' + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + !call routine_1 + !call routine_2 + call routine_3() + +end + +! --- + +subroutine routine_0 + implicit none + use bitmasks ! you need to include the bitmasks_module.f90 features + integer :: i,a,j,m,i_ok + integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2,degree + + integer(bit_kind), allocatable :: det_i(:,:) + double precision :: hmono,htwoe,hthree,htilde_ij,phase + double precision :: same, op, tot, accu + allocate(det_i(N_int,2)) + s1 = 1 + accu = 0.d0 + do i = 1, elec_alpha_num ! occupied + do a = elec_alpha_num+1, mo_num ! virtual + det_i = ref_bitmask + call do_single_excitation(det_i,i,a,s1,i_ok) + if(i_ok == -1)then + print*,'PB !!' + print*,i,a + stop + endif +! call debug_det(det_i,N_int) + call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int) + call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij) + op = fock_3_mat_a_op_sh_bi_orth(a,i) + same = fock_3_mat_a_sa_sh_bi_orth(a,i) +! same = 0.d0 + tot = same + op + if(dabs(tot - phase*hthree).gt.1.d-10)then + print*,'------' + print*,i,a,phase + print*,'hthree = ',phase*hthree + print*,'fock = ',tot + print*,'same,op= ',same,op + print*,dabs(tot - phase*hthree) + stop + endif + accu += dabs(tot - phase*hthree) + enddo + enddo + print*,'accu = ',accu + +end subroutine routine_0 + +! --- + +subroutine routine_1 + + implicit none + integer :: i, a + double precision :: accu + + accu = 0.d0 + do i = 1, mo_num + do a = 1, mo_num + accu += dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) + !if(dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) .gt. 1.d-10)then + print*, i, a + print*, dabs( fock_3_mat_a_op_sh_bi_orth_old(a,i) - fock_3_mat_a_op_sh_bi_orth(a,i) ) & + , fock_3_mat_a_op_sh_bi_orth_old(a,i), fock_3_mat_a_op_sh_bi_orth(a,i) + !endif + enddo + enddo + + print *, 'accu = ', accu + +end subroutine routine_1 + +! --- + +subroutine routine_2 + + implicit none + integer :: i, a + double precision :: accu + + accu = 0.d0 + do i = 1, mo_num + do a = 1, mo_num + accu += dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) + !if(dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) .gt. 1.d-10)then + print*, i, a + print*, dabs( fock_3_mat_a_sa_sh_bi_orth_old(a,i) - fock_3_mat_a_sa_sh_bi_orth(a,i) ) & + , fock_3_mat_a_sa_sh_bi_orth_old(a,i), fock_3_mat_a_sa_sh_bi_orth(a,i) + !endif + enddo + enddo + + print *, 'accu = ', accu + +end subroutine routine_2 + +! --- + +subroutine routine_3() + + use bitmasks ! you need to include the bitmasks_module.f90 features + + implicit none + integer :: i, a, i_ok, s1 + double precision :: hmono, htwoe, hthree, htilde_ij + double precision :: err_ai, err_tot + integer(bit_kind), allocatable :: det_i(:,:) + + allocate(det_i(N_int,2)) + + err_tot = 0.d0 + + s1 = 1 + + det_i = ref_bitmask + call debug_det(det_i, N_int) + print*, ' HF det' + call debug_det(det_i, N_int) + + do i = 1, elec_alpha_num ! occupied + do a = elec_alpha_num+1, mo_num ! virtual + + + det_i = ref_bitmask + call do_single_excitation(det_i, i, a, s1, i_ok) + if(i_ok == -1) then + print*, 'PB !!' + print*, i, a + stop + endif + !print*, ' excited det' + !call debug_det(det_i, N_int) + + call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij) + err_ai = dabs(htilde_ij) + if(err_ai .gt. 1d-7) then + print*, ' warning on', i, a + print*, hmono, htwoe, htilde_ij + endif + err_tot += err_ai + + write(22, *) htilde_ij + enddo + enddo + + print *, ' err_tot = ', err_tot + + deallocate(det_i) + +end subroutine routine_3 + +! ---