added the 3body contribution in the energy for TCSCF

src/bi_ort_ints/semi_num_ints_mo.irp.f

@@ -127,6 +127,9 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
   implicit none
   integer :: ipoint
 
+ print*,'providing int2_grad1_u12_bimo_transp'
+ double precision :: wall0, wall1
+ call wall_time(wall0)
  !$OMP PARALLEL         &
  !$OMP DEFAULT (NONE)   &
  !$OMP PRIVATE (ipoint) & 
@@ -142,6 +145,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
   enddo
  !$OMP END DO
  !$OMP END PARALLEL
+ call wall_time(wall1)
+ print*,'Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
 
 END_PROVIDER 
 

src/tc_bi_ortho/slater_tc_3e.irp.f

@@ -49,6 +49,8 @@ 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)
@@ -60,11 +62,14 @@ 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)
@@ -74,11 +79,14 @@ 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)
@@ -87,11 +95,14 @@ 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)
@@ -100,9 +111,11 @@ 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/test_tc_fock.irp.f

@@ -13,105 +13,44 @@ program test_tc_fock
 
   !call routine_1
   !call routine_2
-  call routine_3()
+!  call routine_3()
 
+ call test_3e
 end
 
 ! ---
 
-subroutine routine_0
+subroutine test_3e
  implicit none
-  use bitmasks ! you need to include the bitmasks_module.f90 features
+ double precision :: integral_aaa,integral_aab,integral_abb,integral_bbb,accu
- integer :: i,a,j,m,i_ok
+ double precision ::  hmono, htwoe, hthree, htot
- integer :: exc(0:2,2,2),h1,p1,s1,h2,p2,s2,degree
+ call htilde_mu_mat_bi_ortho(ref_bitmask, ref_bitmask, N_int, hmono, htwoe, hthree, htot)
-
+! call diag_htilde_three_body_ints_bi_ort(N_int, ref_bitmask, hthree)
- integer(bit_kind), allocatable :: det_i(:,:)
+ print*,'hmono = ',hmono
- double precision :: hmono,htwoe,hthree,htilde_ij,phase
+ print*,'htwoe = ',htwoe
- double precision :: same, op, tot, accu
+ print*,'hthree= ',hthree
- allocate(det_i(N_int,2))
+ print*,'htot  = ',htot
- s1 = 1
+ print*,''
- accu = 0.d0
+ print*,''
- do i = 1, elec_alpha_num ! occupied
+ print*,'TC_one= ',TC_HF_one_electron_energy
-  do a = elec_alpha_num+1, mo_num ! virtual 
+ print*,'TC_two= ',TC_HF_two_e_energy
-    det_i = ref_bitmask
+ print*,'TC_3e = ',diag_three_elem_hf
-    call do_single_excitation(det_i,i,a,s1,i_ok)
+ print*,'TC_tot= ',TC_HF_energy
-    if(i_ok == -1)then
+ print*,''
-     print*,'PB !!'
+ print*,''
-     print*,i,a
+ call give_aaa_contrib(integral_aaa)
-     stop
+ print*,'integral_aaa = ',integral_aaa
-    endif
+ call give_aab_contrib(integral_aab)
-!    call debug_det(det_i,N_int)
+ print*,'integral_aab = ',integral_aab
-    call get_excitation(ref_bitmask,det_i,exc,degree,phase,N_int)
+ call give_abb_contrib(integral_abb)
-    call htilde_mu_mat_bi_ortho(det_i,ref_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+ print*,'integral_abb = ',integral_abb
-    op   = fock_3_mat_a_op_sh_bi_orth(a,i)
+ call give_bbb_contrib(integral_bbb)
-    same = fock_3_mat_a_sa_sh_bi_orth(a,i)
+ print*,'integral_bbb = ',integral_bbb
-!    same = 0.d0
+ accu = integral_aaa + integral_aab + integral_abb + integral_bbb
-    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
+ print*,'delta = ',hthree - accu
 
-end subroutine routine_0
+end
-
-! ---
-
-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()
 

src/tc_scf/fock_three.irp.f

@@ -74,8 +74,11 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf]
  implicit none
  integer :: i,j,k,ipoint,mm
  double precision :: contrib,weight,four_third,one_third,two_third,exchange_int_231
+ print*,'providing diag_three_elem_hf'
+ if(.not.three_body_h_tc)then
+  diag_three_elem_hf = 0.d0
+ else
   if(.not.bi_ortho)then
-  if(three_body_h_tc)then
    one_third = 1.d0/3.d0
    two_third = 2.d0/3.d0
    four_third = 4.d0/3.d0
@@ -102,10 +105,14 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf]
    enddo
    diag_three_elem_hf = - diag_three_elem_hf
   else
-   diag_three_elem_hf = 0.D0
+   double precision :: integral_aaa,hthree, integral_aab,integral_abb,integral_bbb
+   provide mo_l_coef mo_r_coef
+   call give_aaa_contrib(integral_aaa)
+   call give_aab_contrib(integral_aab)
+   call give_abb_contrib(integral_abb)
+   call give_bbb_contrib(integral_bbb)
+   diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb
   endif
- else
-   diag_three_elem_hf = 0.D0
  endif
 END_PROVIDER 
 

src/tc_scf/fock_three_bi_ortho_new_new.irp.f

@@ -154,6 +154,7 @@ BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)]
 END_PROVIDER 
 
 subroutine contrib_3e_sss(a,i,j,k,integral)
+ implicit none
  integer, intent(in) :: a,i,j,k
  BEGIN_DOC
  ! returns the pure same spin contribution to F(a,i) from two orbitals j,k
@@ -173,6 +174,7 @@ subroutine contrib_3e_sss(a,i,j,k,integral)
 end
 
 subroutine contrib_3e_soo(a,i,j,k,integral)
+ implicit none
  integer, intent(in) :: a,i,j,k
  BEGIN_DOC
  ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k

src/tc_scf/tc_scf.irp.f

@@ -81,7 +81,7 @@ subroutine routine_scf()
   print*,'TC HF total energy = ', TC_HF_energy
   print*,'TC HF 1 e   energy = ', TC_HF_one_electron_energy
   print*,'TC HF 2 e   energy = ', TC_HF_two_e_energy
-  if(.not. bi_ortho)then
+  if(three_body_h_tc)then
    print*,'TC HF 3 body       = ', diag_three_elem_hf
   endif
   print*,'***'
@@ -124,6 +124,9 @@ subroutine routine_scf()
       print*,'TC HF total energy = ', TC_HF_energy
       print*,'TC HF 1 e   energy = ', TC_HF_one_electron_energy
       print*,'TC HF 2 non hermit = ', TC_HF_two_e_energy
+      if(three_body_h_tc)then
+       print*,'TC HF 3 body       = ', diag_three_elem_hf
+      endif
       print*,'***'
       e_delta = dabs( TC_HF_energy - e_save )
       print*, 'it, delta E = ', it, e_delta

src/tc_scf/three_e_energy_bi_ortho.irp.f

@@ -0,0 +1,174 @@
+
+subroutine contrib_3e_diag_sss(i,j,k,integral)
+ implicit none
+ integer, intent(in) :: i,j,k
+ BEGIN_DOC
+ ! returns the pure same spin contribution to diagonal matrix element of 3e term
+ END_DOC
+ double precision, intent(out) :: integral
+ double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+ call  give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j >
+ call  give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int)      ! < i k j | j i k >
+ call  give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i >
+ integral = direct_int + c_3_int + c_minus_3_int 
+ ! negative terms :: exchange contrib
+ call  give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 
+ call  give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23
+ call  give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12
+ integral += - exch_13_int - exch_23_int  - exch_12_int 
+ integral = -integral
+end
+
+subroutine contrib_3e_diag_soo(i,j,k,integral)
+ implicit none
+ integer, intent(in) :: i,j,k
+ BEGIN_DOC
+ ! returns the pure same spin contribution to diagonal matrix element of 3e term
+ END_DOC
+ double precision, intent(out) :: integral
+ double precision :: direct_int, exch_23_int
+ call  give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int) ! < i k j | i k j >
+ call  give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)! < i k j | i j k > : E_23
+ integral = direct_int - exch_23_int 
+ integral = -integral
+end
+
+
+subroutine give_aaa_contrib_bis(integral_aaa)
+ implicit none
+ double precision, intent(out) :: integral_aaa
+ double precision :: integral
+ integer :: i,j,k
+ integral_aaa = 0.d0
+ do i = 1, elec_alpha_num
+  do j = i+1, elec_alpha_num
+   do k = j+1, elec_alpha_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_aaa += integral
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine give_aaa_contrib(integral_aaa)
+ implicit none
+ double precision, intent(out) :: integral_aaa
+ double precision :: integral
+ integer :: i,j,k
+ integral_aaa = 0.d0
+ do i = 1, elec_alpha_num
+  do j = 1, elec_alpha_num
+   do k = 1, elec_alpha_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_aaa += integral
+   enddo
+  enddo
+ enddo
+ integral_aaa *= 1.d0/6.d0 
+end
+
+
+subroutine give_aab_contrib(integral_aab)
+ implicit none
+ double precision, intent(out) :: integral_aab
+ double precision :: integral
+ integer :: i,j,k
+ integral_aab = 0.d0
+ do i = 1, elec_beta_num
+  do j = 1, elec_alpha_num
+   do k = 1, elec_alpha_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_aab += integral
+   enddo
+  enddo
+ enddo
+ integral_aab *= 0.5d0
+end
+
+
+subroutine give_aab_contrib_bis(integral_aab)
+ implicit none
+ double precision, intent(out) :: integral_aab
+ double precision :: integral
+ integer :: i,j,k
+ integral_aab = 0.d0
+ do i = 1, elec_beta_num
+  do j = 1, elec_alpha_num
+   do k = j+1, elec_alpha_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_aab += integral
+   enddo
+  enddo
+ enddo
+end
+
+
+subroutine give_abb_contrib(integral_abb)
+ implicit none
+ double precision, intent(out) :: integral_abb
+ double precision :: integral
+ integer :: i,j,k
+ integral_abb = 0.d0
+ do i = 1, elec_alpha_num
+  do j = 1, elec_beta_num
+   do k = 1, elec_beta_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_abb += integral
+   enddo
+  enddo
+ enddo
+ integral_abb *= 0.5d0
+end
+
+subroutine give_abb_contrib_bis(integral_abb)
+ implicit none
+ double precision, intent(out) :: integral_abb
+ double precision :: integral
+ integer :: i,j,k
+ integral_abb = 0.d0
+ do i = 1, elec_alpha_num
+  do j = 1, elec_beta_num
+   do k = j+1, elec_beta_num
+    call contrib_3e_diag_soo(i,j,k,integral)
+    integral_abb += integral
+   enddo
+  enddo
+ enddo
+end
+
+subroutine give_bbb_contrib_bis(integral_bbb)
+ implicit none
+ double precision, intent(out) :: integral_bbb
+ double precision :: integral
+ integer :: i,j,k
+ integral_bbb = 0.d0
+ do i = 1, elec_beta_num
+  do j = i+1, elec_beta_num
+   do k = j+1, elec_beta_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_bbb += integral
+   enddo
+  enddo
+ enddo
+
+end
+
+subroutine give_bbb_contrib(integral_bbb)
+ implicit none
+ double precision, intent(out) :: integral_bbb
+ double precision :: integral
+ integer :: i,j,k
+ integral_bbb = 0.d0
+ do i = 1, elec_beta_num
+  do j = 1, elec_beta_num
+   do k = 1, elec_beta_num
+    call contrib_3e_diag_sss(i,j,k,integral)
+    integral_bbb += integral
+   enddo
+  enddo
+ enddo
+ integral_bbb *= 1.d0/6.d0 
+end
+
+