optimized all matrix elements with three body terms

src/tc_bi_ortho/slater_tc_opt.irp.f

@@ -28,19 +28,13 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
   if(degree.gt.2) return
 
   if(degree == 0)then
-    call diag_htilde_mu_mat_fock_bi_ortho(Nint, key_i, hmono, htwoe, hthree, htot)
+    call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
   else if (degree == 1)then
-    call single_htilde_mu_mat_fock_bi_ortho (Nint,key_j, key_i , hmono, htwoe, hthree, htot)
+    call single_htilde_mu_mat_fock_bi_ortho(Nint,key_j, key_i , hmono, htwoe, hthree, htot)
   else if(degree == 2)then
-    call double_htilde_mu_mat_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htot)
+    call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
-    if(three_body_h_tc) then
-     if(.not.double_normal_ord) then
-       call double_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
-     endif
-    endif
   endif
 
-  htot = hmono + htwoe + hthree
   if(degree==0) then
     htot += nuclear_repulsion
   endif

src/tc_bi_ortho/slater_tc_opt_double.irp.f

@@ -1,3 +1,200 @@
+
+subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
+
+  BEGIN_DOC
+  ! <key_j | H_tilde | key_i> 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, hthree, 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
+
+
+  call get_excitation_degree(key_i, key_j, degree, Nint)
+
+  hmono  = 0.d0
+  htwoe  = 0.d0
+  hthree = 0.d0
+  htot   = 0.d0
+
+  if(degree.ne.2)then
+   return
+  endif
+
+  call bitstring_to_list_ab(key_i, occ, Ne, Nint)
+  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(three_body_h_tc)then
+     if(.not.double_normal_ord)then
+      call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
+     elseif(double_normal_ord.and.+Ne(1).gt.2)then
+      htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ???
+     endif
+    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(three_body_h_tc)then
+    if(.not.double_normal_ord)then
+     call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
+    elseif(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
+  endif
+  hthree *= phase
+  htwoe  *= phase
+  htot    =  htwoe + hthree
+
+end
+
+
+
+subroutine three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree)
+ implicit none
+ integer(bit_kind), intent(in) :: key_i(N_int,2)
+ integer, intent(in) :: h1,h2,p1,p2,s1,s2
+ double precision, intent(out) :: hthree
+ integer :: nexc(2),i,ispin,na,nb
+ integer(bit_kind) :: hole(N_int,2)
+ integer(bit_kind) :: particle(N_int,2)
+ integer :: occ_hole(N_int*bit_kind_size,2)
+ integer :: occ_particle(N_int*bit_kind_size,2)
+ integer :: n_occ_ab_hole(2),n_occ_ab_particle(2)
+ integer(bit_kind)              :: det_tmp(N_int,2)
+ integer :: ipart, ihole
+ double precision :: direct_int, exchange_int, three_e_double_parrallel_spin
+
+  nexc(1) = 0
+  nexc(2) = 0
+  !! Get all the holes and particles of key_i with respect to the ROHF determinant
+  do i=1,N_int
+    hole(i,1)     = xor(key_i(i,1),ref_bitmask(i,1))
+    hole(i,2)     = xor(key_i(i,2),ref_bitmask(i,2))
+    particle(i,1) = iand(hole(i,1),key_i(i,1))
+    particle(i,2) = iand(hole(i,2),key_i(i,2))
+    hole(i,1)     = iand(hole(i,1),ref_bitmask(i,1))
+    hole(i,2)     = iand(hole(i,2),ref_bitmask(i,2))
+    nexc(1)       = nexc(1) + popcnt(hole(i,1))
+    nexc(2)       = nexc(2) + popcnt(hole(i,2))
+  enddo
+  integer                        :: tmp(2)
+  !DIR$ FORCEINLINE
+  call bitstring_to_list_ab(particle, occ_particle, tmp, N_int)
+  ASSERT (tmp(1) == nexc(1)) ! Number of particles alpha
+  ASSERT (tmp(2) == nexc(2)) ! Number of particle beta 
+  !DIR$ FORCEINLINE
+  call bitstring_to_list_ab(hole, occ_hole, tmp, N_int)
+  ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
+  ASSERT (tmp(2) == nexc(2)) ! Number of holes beta 
+  if(s1==s2.and.s1==1)then
+   !!!!!!!!!!!!!!!!!!!!!!!!!! alpha/alpha double exc
+   hthree = eff_2_e_from_3_e_aa(p2,p1,h2,h1) 
+   if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant 
+    !!!!!!!! the matrix element is already exact 
+    !!!!!!!! else you need to take care of holes and particles 
+    !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
+    ispin = 1 ! i==alpha ==> pure same spin terms
+    do i = 1, nexc(ispin) ! number of couple of holes/particles 
+     ipart=occ_particle(i,ispin)
+     hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1)
+     ihole=occ_hole(i,ispin)
+     hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1)
+    enddo
+    ispin = 2 ! i==beta ==> alpha/alpha/beta terms
+    do i = 1, nexc(ispin) ! number of couple of holes/particles 
+     ! exchange between (h1,p1) and (h2,p2)
+     ipart=occ_particle(i,ispin)
+     direct_int  = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
+     exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
+     hthree += direct_int - exchange_int
+     ihole=occ_hole(i,ispin)
+     direct_int  = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
+     exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
+     hthree -= direct_int - exchange_int
+    enddo
+   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  elseif(s1==s2.and.s1==2)then 
+   !!!!!!!!!!!!!!!!!!!!!!!!!! beta/beta double exc
+   hthree = eff_2_e_from_3_e_bb(p2,p1,h2,h1)
+   if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant 
+   !!!!!!!! the matrix element is already exact 
+   !!!!!!!! else you need to take care of holes and particles 
+   !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
+   ispin = 2 ! i==beta  ==> pure same spin terms
+   do i = 1, nexc(ispin) ! number of couple of holes/particles 
+    ipart=occ_particle(i,ispin)
+    hthree += three_e_double_parrallel_spin(ipart,p2,h2,p1,h1)
+    ihole=occ_hole(i,ispin)
+    hthree -= three_e_double_parrallel_spin(ihole,p2,h2,p1,h1)
+   enddo
+   ispin = 1 ! i==alpha==> beta/beta/alpha terms
+   do i = 1, nexc(ispin) ! number of couple of holes/particles 
+    ! exchange between (h1,p1) and (h2,p2)
+    ipart=occ_particle(i,ispin)
+    direct_int  = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_exch12_bi_ort(ipart,p2,h2,p1,h1)
+    hthree += direct_int - exchange_int
+    ihole=occ_hole(i,ispin)
+    direct_int  = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_exch12_bi_ort(ihole,p2,h2,p1,h1)
+    hthree -= direct_int - exchange_int
+   enddo
+  else                         ! (h1,p1) == alpha/(h2,p2) == beta 
+   hthree = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
+   if(nexc(1)+nexc(2) ==0)return !! if you're on the reference determinant 
+   !!!!!!!! the matrix element is already exact 
+   !!!!!!!! else you need to take care of holes and particles 
+   !!!!!!!!!!!!! Holes and particles !!!!!!!!!!!!!!!!!!!!!!!
+   ispin = 1 ! i==alpha ==> alpha/beta/alpha terms 
+   do i = 1, nexc(ispin) ! number of couple of holes/particles 
+    ! exchange between (h1,p1) and i
+    ipart=occ_particle(i,ispin)
+    direct_int  = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_exch13_bi_ort(ipart,p2,h2,p1,h1)
+    hthree += direct_int - exchange_int
+    ihole=occ_hole(i,ispin)
+    direct_int  = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_exch13_bi_ort(ihole,p2,h2,p1,h1)
+    hthree -= direct_int - exchange_int
+   enddo
+   ispin = 2 ! i==beta  ==> alpha/beta/beta  terms 
+   do i = 1, nexc(ispin) ! number of couple of holes/particles 
+    ! exchange between (h2,p2) and i
+    ipart=occ_particle(i,ispin)
+    direct_int  = three_e_5_idx_direct_bi_ort(ipart,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_exch23_bi_ort(ipart,p2,h2,p1,h1)
+    hthree += direct_int - exchange_int
+    ihole=occ_hole(i,ispin)
+    direct_int  = three_e_5_idx_direct_bi_ort(ihole,p2,h2,p1,h1)
+    exchange_int = three_e_5_idx_exch23_bi_ort(ihole,p2,h2,p1,h1)
+    hthree -= direct_int - exchange_int
+   enddo
+  endif
+end
+
+
 BEGIN_PROVIDER [ double precision, eff_2_e_from_3_e_ab, (mo_num, mo_num, mo_num, mo_num)]
  implicit none
  BEGIN_DOC

src/tc_bi_ortho/test_normal_order.irp.f

@@ -40,13 +40,15 @@ subroutine test
      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 = eff_2_e_from_3_e_ab(p2,p1,h2,h1)
+     call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal)
+!     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
+stop
 
 !    p2=6
 !    p1=5

src/tc_bi_ortho/test_tc_bi_ortho.irp.f

@@ -50,13 +50,14 @@ subroutine test_slater_tc_opt
      accu += dabs(htot-hnewtot) 
 !    if(dabs(hthree-hnewthree).gt.1.d-8.or.dabs(hthree-hnewthree).gt.dabs(hthree))then
     if(dabs(htot-hnewtot).gt.1.d-8.or.dabs(htot-hnewtot).gt.dabs(htot))then
-     print*,j,i
+     print*,j,i,degree
      call debug_det(psi_det(1,1,i),N_int)
      call debug_det(psi_det(1,1,j),N_int)
      print*,htot,hnewtot,dabs(htot-hnewtot) 
 !     print*,hthree,hnewthree,dabs(hthree-hnewthree) 
      stop
     endif
+    print*,htot,hnewtot,dabs(htot-hnewtot) 
    endif
   enddo
  enddo