beginning to optimize the double excitations for TC

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -7,7 +7,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num,
   !
   ! matrix element of the -L  three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
   !
-  ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mjk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
+  ! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
   !
   ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
   END_DOC
@@ -202,7 +202,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num,
   !
   ! matrix element of the -L  three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
   !
-  ! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
+  ! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|ijm> ::: notice that i is the RIGHT MO and k is the LEFT MO
   !
   ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
   !
@@ -251,7 +251,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num,
   !
   ! matrix element of the -L  three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
   !
-  ! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
+  ! three_e_5_idx_exch12_bi_ort(m,l,j,k,i) = <mlk|-L|mij> ::: notice that i is the RIGHT MO and k is the LEFT MO
   !
   ! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
   !

src/tc_bi_ortho/slater_tc_opt_double.irp.f

@@ -0,0 +1,57 @@
+BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! eff_2_e_from_3_e_ab(p2,p1,h2,h1) = Effective Two-electron operator for alpha/beta double excitations 
+!
+! from contraction with HF density = a^{dagger}_p1_alpha a^{dagger}_p2_beta a_h2_beta a_h1_alpha
+ END_DOC
+ integer :: i,h1,p1,h2,p2
+ integer :: hh1,hh2,pp1,pp2,m,mm
+ integer                        :: Ne(2)
+ integer,           allocatable :: occ(:,:)
+ double precision :: contrib
+ allocate( occ(N_int*bit_kind_size,2) )
+ call bitstring_to_list_ab(ref_bitmask,occ,Ne,N_int)
+ eff_2_e_from_3_e_ab = 0.d0
+  do hh1 = 1, n_act_orb !! alpha 
+    h1 = list_act(hh1) 
+    do hh2 = 1, n_act_orb !! beta 
+      h2 = list_act(hh2) 
+      do pp1 = 1, n_act_orb !! alpha
+        p1 = list_act(pp1)
+        do pp2 = 1, n_act_orb !! beta 
+          p2 = list_act(pp2)
+          call give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
+          eff_2_e_from_3_e_ab(p2,p1,h2,h1) = contrib
+        enddo
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER 
+
+subroutine give_contrib_for_abab(h1,h2,p1,p2,occ,Ne,contrib)
+ implicit none
+ integer, intent(in) :: h1,h2,p1,p2,occ(N_int*bit_kind_size,2),Ne(2)
+ double precision, intent(out) :: contrib
+ integer :: mm,m 
+ double precision :: direct_int, exchange_int
+ !! h1,p1 == alpha 
+ !! h2,p2 == beta
+ contrib = 0.d0
+ do mm = 1, Ne(1) !! alpha 
+   m = occ(m,1)
+   direct_int   = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) 
+   ! exchange between (h1,p1) and m
+   exchange_int = three_e_5_idx_exch13_bi_ort(mm,p2,h2,p1,h1)
+   contrib += direct_int - exchange_int
+ enddo
+
+ do mm = 1, Ne(2) !! beta
+   m = occ(m,2)
+   direct_int   = three_e_5_idx_direct_bi_ort(mm,p2,h2,p1,h1) 
+   ! exchange between (h2,p2) and m
+   exchange_int = three_e_5_idx_exch23_bi_ort(mm,p2,h2,p1,h1)
+   contrib += direct_int - exchange_int
+ enddo
+end

src/tc_bi_ortho/test_normal_order.irp.f

@@ -10,6 +10,7 @@ program test_normal_order
   read_wf = .True.
   touch read_wf
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+  call provide_all_three_ints_bi_ortho
   call test
 end
 
@@ -28,7 +29,7 @@ subroutine test
  s2 = 2
  accu = 0.d0
  do h1 = 1, elec_beta_num
-  do p1 = elec_beta_num+1, mo_num
+  do p1 = elec_alpha_num+1, mo_num
    do h2 = 1, elec_beta_num
     do p2 = elec_beta_num+1, mo_num
      det_i = ref_bitmask
@@ -38,36 +39,37 @@ subroutine test
      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)
+!     normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
+     normal = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
      accu += dabs(hthree-normal)
     enddo
    enddo
   enddo
  enddo
- print*,'accu opposite spin = ',accu
+print*,'accu opposite spin = ',accu
 
- s1 = 2 
+!s1 = 2 
- s2 = 2
+!s2 = 2
- accu = 0.d0
+!accu = 0.d0
- do h1 = 1, elec_beta_num
+!do h1 = 1, elec_beta_num
-  do p1 = elec_beta_num+1, mo_num
+! do p1 = elec_beta_num+1, mo_num
-   do h2 = h1+1, elec_beta_num
+!  do h2 = h1+1, elec_beta_num
-    do p2 = elec_beta_num+1, mo_num
+!   do p2 = elec_beta_num+1, mo_num
-     det_i = ref_bitmask
+!    det_i = ref_bitmask
-     call do_single_excitation(det_i,h1,p1,s1,i_ok)
+!    call do_single_excitation(det_i,h1,p1,s1,i_ok)
-     call do_single_excitation(det_i,h2,p2,s2,i_ok)
+!    call do_single_excitation(det_i,h2,p2,s2,i_ok)
-     if(i_ok.ne.1)cycle
+!    if(i_ok.ne.1)cycle
-     call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+!    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_degree(ref_bitmask,det_i,degree,N_int)
-     call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
+!    call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
-     hthree *= phase
+!    hthree *= phase
-     normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
+!    normal = normal_two_body_bi_orth_aa_bb(p2,h2,p1,h1)
-     accu += dabs(hthree-normal)
+!    accu += dabs(hthree-normal)
-    enddo
+!   enddo
-   enddo
+!  enddo
-  enddo
+! enddo
- enddo
+!enddo
- print*,'accu same spin = ',accu
+!print*,'accu same spin = ',accu
 end