diagonal matrix elements work with 3-e a la fock

src/bi_ort_ints/three_body_ijm.irp.f

@@ -60,7 +60,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num
   !
   ! matrix element of the -L  three-body operator ON A BI ORTHONORMAL BASIS for the first cyclic permutation 
   !
-  ! three_e_3_idx_direct_bi_ort(m,j,i) = <mji|-L|jim>
+  ! three_e_3_idx_cycle_1_bi_ort(m,j,i) = <mji|-L|jim>
   !
   ! 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_3e.irp.f

@@ -49,8 +49,6 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
 
   if(Ne(1)+Ne(2).ge.3)then
 !!  ! alpha/alpha/beta three-body
-  double precision :: accu 
-  accu = 0.d0
    do i = 1, Ne(1)
     ii = occ(i,1) 
     do j = i+1, Ne(1)
@@ -62,14 +60,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
       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
-      accu += direct_int - exchange_int
      enddo
     enddo
    enddo
-   !print*,'aab = ',accu
   
    ! beta/beta/alpha three-body
-  accu = 0.d0
    do i = 1, Ne(2)
     ii = occ(i,2) 
     do j = i+1, Ne(2)
@@ -79,14 +74,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
       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
-      accu += direct_int - exchange_int
      enddo
     enddo
    enddo
-   !print*,'abb = ',accu
 
    ! alpha/alpha/alpha three-body
-   accu = 0.d0
    do i = 1, Ne(1)
     ii = occ(i,1) ! 1
     do j = i+1, Ne(1)
@@ -95,14 +87,11 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
       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
-      accu += three_e_diag_parrallel_spin(mm,jj,ii) 
      enddo
     enddo
    enddo
-   !print*,'aaa = ',accu
 
    ! beta/beta/beta three-body
-   accu = 0.d0
    do i = 1, Ne(2)
     ii = occ(i,2) ! 1
     do j = i+1, Ne(2)
@@ -111,11 +100,9 @@ subroutine diag_htilde_three_body_ints_bi_ort(Nint, key_i, hthree)
       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
-      accu += three_e_diag_parrallel_spin(mm,jj,ii) 
      enddo
     enddo
    enddo
-   !print*,'bbb = ',accu
   endif
 
 end

src/tc_bi_ortho/slater_tc_opt.irp.f

@@ -55,6 +55,9 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree,
   enddo
 
   if (nexc(1)+nexc(2) == 0) then
+    hmono = ref_tc_energy_1e
+    htwoe = ref_tc_energy_2e
+    hthree= ref_tc_energy_3e
     htot = ref_tc_energy_tot
     return
   endif
@@ -75,7 +78,6 @@ subroutine diag_htilde_mu_mat_fock_bi_ortho(Nint, det_in, hmono, htwoe, hthree,
   hmono = ref_tc_energy_1e
   htwoe = ref_tc_energy_2e 
   hthree= ref_tc_energy_3e
-
   do ispin=1,2
     na = elec_num_tab(ispin)
     nb = elec_num_tab(iand(ispin,1)+1)
@@ -110,7 +112,9 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
 
   integer                        :: occ(Nint*bit_kind_size,2)
   integer                        :: other_spin
-  integer                        :: k,l,i
+  integer                        :: k,l,i,jj,mm,j,m
+  double precision :: three_e_diag_parrallel_spin, direct_int, exchange_int
+  
 
   if (iorb < 1) then
     print *,  irp_here, ': iorb < 1'
@@ -151,6 +155,39 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
   do i=1,nb
     htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
   enddo
+
+  if(three_body_h_tc)then
+   !!!!! 3-e part 
+   !! same-spin/same-spin
+   do j = 1, na
+    jj = occ(j,ispin)
+    do m = j+1, na
+     mm = occ(m,ispin)
+     hthree += three_e_diag_parrallel_spin(mm,jj,iorb)
+    enddo
+   enddo
+   !! same-spin/oposite-spin
+   do j = 1, na
+    jj = occ(j,ispin)
+    do m = 1, nb
+     mm = occ(m,other_spin)
+     direct_int   = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+     exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+     hthree += direct_int - exchange_int
+    enddo
+   enddo
+   !! oposite-spin/opposite-spin
+    do j = 1, nb
+     jj = occ(j,other_spin) 
+     do m = j+1, nb 
+      mm = occ(m,other_spin) 
+      direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+      exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+      hthree += direct_int - exchange_int
+     enddo
+    enddo
+  endif
+
   na = na+1
 end
 
@@ -173,9 +210,10 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
   integer(bit_kind), intent(inout) :: key(Nint,2)
   double precision, intent(inout) :: hmono,htwoe,hthree
   
+  double precision  :: direct_int, exchange_int, three_e_diag_parrallel_spin
   integer                        :: occ(Nint*bit_kind_size,2)
   integer                        :: other_spin
-  integer                        :: k,l,i
+  integer                        :: k,l,i,jj,mm,j,m
   integer                        :: tmp(2)
 
   ASSERT (iorb > 0)
@@ -205,5 +243,37 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
     htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
   enddo
 
+  if(three_body_h_tc)then
+   !!!!! 3-e part 
+   !! same-spin/same-spin
+   do j = 1, na
+    jj = occ(j,ispin)
+    do m = j+1, na
+     mm = occ(m,ispin)
+     hthree -= three_e_diag_parrallel_spin(mm,jj,iorb)
+    enddo
+   enddo
+   !! same-spin/oposite-spin
+   do j = 1, na
+    jj = occ(j,ispin)
+    do m = 1, nb
+     mm = occ(m,other_spin)
+     direct_int   = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+     exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+     hthree -= (direct_int - exchange_int)
+    enddo
+   enddo
+   !! oposite-spin/opposite-spin
+    do j = 1, nb
+     jj = occ(j,other_spin) 
+     do m = j+1, nb 
+      mm = occ(m,other_spin) 
+      direct_int   = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+      exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+      hthree -= (direct_int - exchange_int)
+     enddo
+    enddo
+  endif
+
 end
 

src/tc_bi_ortho/test_tc_bi_ortho.irp.f

@@ -19,28 +19,36 @@ subroutine test_slater_tc_opt
  integer :: i,j
  double precision :: hmono, htwoe, htot, hthree 
  double precision :: hnewmono, hnewtwoe, hnewthnewree, hnewtot
- double precision :: accu ,i_count
+ double precision :: accu_d ,i_count, accu
  accu = 0.d0
+ accu_d = 0.d0
  i_count = 0.d0
  do i = 1, N_det
 ! do i = 14,14
   call diag_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,i), hmono, htwoe, htot)
+  call htilde_mu_mat_bi_ortho(psi_det(1,1,i), psi_det(1,1,i), N_int, hmono, htwoe, hthree, htot)
   call diag_htilde_mu_mat_fock_bi_ortho(N_int, psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot)
+!  print*,hthree,hnewthnewree
+!  print*,htot,hnewtot,dabs(hnewtot-htot)
+  accu_d += dabs(htot-hnewtot) 
+!  if(dabs(htot-hnewtot).gt.1.d-8)then
+   print*,i
+   print*,htot,hnewtot,dabs(htot-hnewtot)
+!  endif
   do j = 1, N_det
-!  do j = 1, 1
    if(i==j)cycle
    call single_htilde_mu_mat_bi_ortho(N_int, psi_det(1,1,j), psi_det(1,1,i), hmono, htwoe, htot)
    call single_htilde_mu_mat_fock_bi_ortho (N_int, psi_det(1,1,j), psi_det(1,1,i), hnewmono, hnewtwoe, hnewthnewree, hnewtot)
    if(dabs(htot).gt.1.d-10)then
-!    if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then
+    if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then
      print*,j,i
      i_count += 1.D0
      print*,htot,hnewtot,dabs(htot-hnewtot) 
      accu += dabs(htot-hnewtot) 
-!    endif
+    endif
    endif
   enddo
  enddo
- print*,'accu = ',accu/i_count
+ print*,'accu_d = ',accu_d/N_det
 
 end