added some comments for normal ordering old

src/tc_bi_ortho/normal_ordered_old.irp.f

@@ -120,6 +120,13 @@ 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
+  BEGIN_DOC
+!  give the contribution for a double excitation of opposite spin BUT averaged over spin
+!
+! it is the average of <p1_down p2_up |h1_down h2_up> and <p1_up p2_down |h1_up h2_down> 
+!
+! because the orbitals h1,h2,p1,p2 are spatial orbitals and therefore can be of different spins
+  END_DOC
 
   implicit none
   integer, intent(in)           :: Nint, h1, h2, p1, p2
@@ -158,7 +165,8 @@ subroutine give_aba_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
     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)
+    hthree += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12) ! spin average
+!    hthree += 1.d0 * int_direct - 1.0d0 * (int_exc_13 + int_exc_12)
   enddo
 
   return

src/tc_bi_ortho/test_normal_order.irp.f

@@ -20,7 +20,7 @@ subroutine test
  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
+ double precision :: hmono,htwoe,hthree,htilde_ij,accu,phase,normal,hthree_tmp
  integer,           allocatable :: occ(:,:)
  allocate( occ(N_int*bit_kind_size,2) )
  call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
@@ -32,15 +32,44 @@ subroutine test
   do p1 = elec_alpha_num+1, mo_num
    do h2 = 1, elec_beta_num
     do p2 = elec_beta_num+1, mo_num
+     hthree = 0.d0
+
      det_i = ref_bitmask
+ s1 = 1 
+ s2 = 2
      call do_single_excitation(det_i,h1,p1,s1,i_ok)
+    if(i_ok.ne.1)cycle
      call do_single_excitation(det_i,h2,p2,s2,i_ok)
-     call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree,htilde_ij)
+    if(i_ok.ne.1)cycle
+     call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,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
+     hthree_tmp *= phase
-!    !normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
+     hthree += 0.5d0 * hthree_tmp
-     call three_comp_two_e_elem(det_i,h1,h2,p1,p2,s1,s2,normal)
+     det_i = ref_bitmask
+ s1 = 2 
+ s2 = 1
+     call do_single_excitation(det_i,h1,p1,s1,i_ok)
+    if(i_ok.ne.1)cycle
+     call do_single_excitation(det_i,h2,p2,s2,i_ok)
+    if(i_ok.ne.1)cycle
+     call htilde_mu_mat_bi_ortho_slow(det_i,HF_bitmask,N_int,hmono,htwoe,hthree_tmp,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_tmp *= phase
+     hthree += 0.5d0 * hthree_tmp
+
+
+!     normal = normal_two_body_bi_orth_ab(p2,h2,p1,h1)
+     call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, normal)
+     if(dabs(hthree).lt.1.d-10)cycle
+     if(dabs(hthree-normal).gt.1.d-10)then
+!      print*,pp2,pp1,hh2,hh1
+      print*,p2,p1,h2,h1
+      print*,hthree,normal,dabs(hthree-normal)
+      stop
+     endif
+!     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
@@ -73,8 +102,8 @@ do h1 = 1, elec_alpha_num
     integer :: hh1, pp1, hh2, pp2, ss1, ss2
     call decode_exc(exc, 2, hh1, pp1, hh2, pp2, ss1, ss2)
     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)
-     normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
+!     normal = eff_2_e_from_3_e_aa(p2,p1,h2,h1)
     if(dabs(hthree).lt.1.d-10)cycle
     if(dabs(hthree-normal).gt.1.d-10)then
      print*,pp2,pp1,hh2,hh1