diff --git a/src/becke_numerical_grid/grid_becke_per_atom.irp.f b/src/becke_numerical_grid/grid_becke_per_atom.irp.f
index 6026350b..75610090 100644
--- a/src/becke_numerical_grid/grid_becke_per_atom.irp.f
+++ b/src/becke_numerical_grid/grid_becke_per_atom.irp.f
@@ -46,8 +46,4 @@ END_PROVIDER
       enddo
     enddo
   enddo
-
-
-
-
 END_PROVIDER 
diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f
index 0a8b4d52..2bc6fc30 100644
--- a/src/dft_utils_in_r/mo_in_r.irp.f
+++ b/src/dft_utils_in_r/mo_in_r.irp.f
@@ -181,3 +181,30 @@
  END_PROVIDER
 
 
+ BEGIN_PROVIDER [double precision, mo_abs_dist_per_nuclei, (mo_num,nucl_num)]
+  implicit none
+  BEGIN_DOC
+  ! mo_abs_dist_per_nuclei(j,i_nucl) = <phi_j| |r-R_nucl| |phi_j>
+  END_DOC
+  integer :: ipoint,i_nucl,m,j
+  double precision :: weight, r(3),r_nucl(3),dist
+  mo_abs_dist_per_nuclei = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r(1) = final_grid_points(1,ipoint)
+   r(2) = final_grid_points(2,ipoint)
+   r(3) = final_grid_points(3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   do i_nucl = 1, nucl_num
+    dist = 0.d0
+    do m = 1, 3
+     r_nucl(m) = r(m) - nucl_coord_transp(m,i_nucl) 
+     r_nucl(m) *= r_nucl(m)
+     dist += r_nucl(m)
+    enddo
+    dist = dsqrt(dist)
+    do j = 1, mo_num
+     mo_abs_dist_per_nuclei(j,i_nucl) += weight * mos_in_r_array(j,ipoint)*mos_in_r_array(j,ipoint) * dist
+    enddo
+   enddo
+  enddo
+ END_PROVIDER 
diff --git a/src/tc_bi_ortho/slater_tc_3e.irp.f b/src/tc_bi_ortho/slater_tc_3e.irp.f
index f4899a80..0d5f8542 100644
--- a/src/tc_bi_ortho/slater_tc_3e.irp.f
+++ b/src/tc_bi_ortho/slater_tc_3e.irp.f
@@ -169,8 +169,6 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
        !   is == ispin  in :::   s1 is is  s1 is is      s1 is is s1 is is
        !                       < h1 j  i | p1 j  i > - < h1 j  i | p1 i j >
        !                                                   
-!      direct_int =   three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor 
-!      exchange_int = three_body_ints_bi_ort(jj,ii,p1,ii,jj,h1) ! USES THE 6-IDX tensor 
        direct_int   = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)  
        exchange_int = three_e_4_idx_exch23_bi_ort(jj,ii,p1,h1) 
        hthree += direct_int - exchange_int
@@ -182,14 +180,10 @@ subroutine single_htilde_three_body_ints_bi_ort(Nint, key_j, key_i, hthree)
       ii = occ(i,ispin)  ! other spin 
       do j = 1, Ne(s1)   ! same spin 
        jj = occ(j,s1)    ! same spin 
-!      direct_int = three_body_ints_bi_ort(jj,ii,p1,jj,ii,h1) ! USES THE 6-IDX tensor 
-!      exchange_int = three_body_ints_bi_ort(jj,ii,p1,h1,ii,jj) ! exchange the spin s1 :: 6-IDX tensor
        direct_int   = three_e_4_idx_direct_bi_ort(jj,ii,p1,h1)
        exchange_int = three_e_4_idx_exch13_bi_ort(jj,ii,p1,h1)
        !              < h1 j  i | p1 j i > - < h1 j i | j p1 i >
        hthree += direct_int - exchange_int
-!       print*,'h1,p1,ii,jj = ',h1,p1,ii,jj
-!       print*,direct_int, exchange_int
       enddo
      enddo
 !
diff --git a/src/tc_bi_ortho/test_tc_fock.irp.f b/src/tc_bi_ortho/test_tc_fock.irp.f
index badc40b6..d585ce6c 100644
--- a/src/tc_bi_ortho/test_tc_fock.irp.f
+++ b/src/tc_bi_ortho/test_tc_fock.irp.f
@@ -120,44 +120,54 @@ subroutine routine_3()
   implicit none
   integer                        :: i, a, i_ok, s1
   double precision               :: hmono, htwoe, hthree, htilde_ij
-  double precision               :: err_ai, err_tot
+  double precision               :: err_ai, err_tot, ref, new
   integer(bit_kind), allocatable :: det_i(:,:)
 
   allocate(det_i(N_int,2))
 
   err_tot = 0.d0
  
-  s1 = 1
+  do s1 = 1, 2
 
-  det_i = ref_bitmask
-  call debug_det(det_i, N_int)
-  print*, ' HF det'
-  call debug_det(det_i, N_int)
-
-  do i = 1, elec_alpha_num ! occupied
-    do a = elec_alpha_num+1, mo_num ! virtual 
-
-
-      det_i = ref_bitmask
-      call do_single_excitation(det_i, i, a, s1, i_ok)
-      if(i_ok == -1) then
-       print*, 'PB !!'
-       print*, i, a
-       stop
-      endif
-      !print*, ' excited det'
-      !call debug_det(det_i, N_int)
-
-      call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
-      err_ai = dabs(htilde_ij)
-      if(err_ai .gt. 1d-7) then
-        print*, ' warning on', i, a
-        print*, hmono, htwoe, htilde_ij
-      endif
-      err_tot += err_ai
-
-      write(22, *) htilde_ij
-    enddo
+   det_i = ref_bitmask
+   call debug_det(det_i, N_int)
+   print*, ' HF det'
+   call debug_det(det_i, N_int)
+ 
+   do i = 1, elec_num_tab(s1)
+     do a = elec_num_tab(s1)+1, mo_num ! virtual 
+ 
+ 
+       det_i = ref_bitmask
+       call do_single_excitation(det_i, i, a, s1, i_ok)
+       if(i_ok == -1) then
+        print*, 'PB !!'
+        print*, i, a
+        stop
+       endif
+       !print*, ' excited det'
+       !call debug_det(det_i, N_int)
+ 
+       call htilde_mu_mat_bi_ortho(det_i, ref_bitmask, N_int, hmono, htwoe, hthree, htilde_ij)
+       if(dabs(hthree).lt.1.d-10)cycle
+       ref = hthree 
+       if(s1 == 1)then
+        new = fock_a_tot_3e_bi_orth(a,i)
+       else if(s1 == 2)then
+        new = fock_b_tot_3e_bi_orth(a,i)
+       endif
+       err_ai = dabs(dabs(ref) - dabs(new))
+       if(err_ai .gt. 1d-7) then
+         print*,'s1 = ',s1
+         print*, ' warning on', i, a
+         print*, ref,new,err_ai
+       endif
+       print*, ref,new,err_ai
+       err_tot += err_ai
+ 
+       write(22, *) htilde_ij
+     enddo
+   enddo
   enddo
 
   print *, ' err_tot = ', err_tot
diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f
index 4907f2c8..ca500fbb 100644
--- a/src/tc_scf/fock_tc.irp.f
+++ b/src/tc_scf/fock_tc.irp.f
@@ -76,13 +76,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)]
 END_PROVIDER 
 ! ---
 
-BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
-  implicit none
-  BEGIN_DOC
- ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
-  END_DOC
-  Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta)
-END_PROVIDER
+!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
+!  implicit none
+!  BEGIN_DOC
+! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
+!  END_DOC
+!  Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta)
+!END_PROVIDER
 
 ! ---
 
@@ -94,6 +94,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
   if(bi_ortho)then
    call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
                          , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
+   if(three_body_h_tc)then
+    Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth
+   endif
   else
    call ao_to_mo(  Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
                  , Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
@@ -110,6 +113,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
   if(bi_ortho)then
    call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
                          , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
+   if(three_body_h_tc)then
+    Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth
+   endif
   else
    call ao_to_mo(  Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
                  , Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
@@ -117,17 +123,17 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
 END_PROVIDER
 
 
-BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)]
-  implicit none
-  BEGIN_DOC
- ! Total alpha+beta  TC Fock matrix : h_c + Two-e^TC terms on the MO basis
-  END_DOC
-  Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta)
-  if(three_body_h_tc) then
-    Fock_matrix_tc_mo_tot += fock_3_mat
-  endif
-  !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10)
-END_PROVIDER 
+!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)]
+!  implicit none
+!  BEGIN_DOC
+! ! Total alpha+beta  TC Fock matrix : h_c + Two-e^TC terms on the MO basis
+!  END_DOC
+!  Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta)
+!  if(three_body_h_tc) then
+!    Fock_matrix_tc_mo_tot += fock_3_mat
+!  endif
+!  !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10)
+!END_PROVIDER 
 
 ! ---
 
diff --git a/src/tc_scf/fock_tc_mo_tot.irp.f b/src/tc_scf/fock_tc_mo_tot.irp.f
new file mode 100644
index 00000000..a99c7698
--- /dev/null
+++ b/src/tc_scf/fock_tc_mo_tot.irp.f
@@ -0,0 +1,121 @@
+
+ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ]
+&BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)]
+   implicit none
+   BEGIN_DOC
+   ! Fock matrix on the MO basis.
+   ! For open shells, the ROHF Fock Matrix is ::
+   !
+   !       |   F-K    |  F + K/2  |    F     |
+   !       |---------------------------------|
+   !       | F + K/2  |     F     |  F - K/2 |
+   !       |---------------------------------|
+   !       |    F     |  F - K/2  |  F + K   |
+   !
+   !
+   ! F = 1/2 (Fa + Fb)
+   !
+   ! K = Fb - Fa
+   !
+   END_DOC
+   integer                        :: i,j,n
+   if (elec_alpha_num == elec_beta_num) then
+     Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha
+   else
+
+     do j=1,elec_beta_num
+       ! F-K
+       do i=1,elec_beta_num !CC
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             - (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F+K/2
+       do i=elec_beta_num+1,elec_alpha_num  !CA
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F
+       do i=elec_alpha_num+1, mo_num !CV
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
+       enddo
+     enddo
+
+     do j=elec_beta_num+1,elec_alpha_num
+       ! F+K/2
+       do i=1,elec_beta_num !AC
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F
+       do i=elec_beta_num+1,elec_alpha_num !AA
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
+       enddo
+       ! F-K/2
+       do i=elec_alpha_num+1, mo_num !AV
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+     enddo
+
+     do j=elec_alpha_num+1, mo_num
+       ! F
+       do i=1,elec_beta_num !VC
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))
+       enddo
+       ! F-K/2
+       do i=elec_beta_num+1,elec_alpha_num !VA
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))&
+             - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+       ! F+K
+       do i=elec_alpha_num+1,mo_num !VV
+         Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) &
+             + (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j))
+       enddo
+     enddo
+
+   endif
+
+   do i = 1, mo_num
+     Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i)
+   enddo
+
+
+   if(frozen_orb_scf)then
+     integer                        :: iorb,jorb
+     do i = 1, n_core_orb
+      iorb = list_core(i)
+      do j = 1, n_act_orb
+       jorb = list_act(j)
+       Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+       Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
+      enddo
+     enddo
+   endif
+
+   if(no_oa_or_av_opt)then
+     do i = 1, n_act_orb
+       iorb = list_act(i)
+       do j = 1, n_inact_orb
+         jorb = list_inact(j)
+         Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+         Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_virt_orb
+         jorb = list_virt(j)
+         Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+         Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0
+       enddo
+       do j = 1, n_core_orb
+         jorb = list_core(j)
+         Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0
+         Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0                                                                                                                 
+       enddo
+     enddo
+   endif
+  if(.not.bi_ortho .and. three_body_h_tc)then
+   Fock_matrix_tc_mo_tot += fock_3_mat
+  endif
+
+END_PROVIDER
+
diff --git a/src/tc_scf/fock_three.irp.f b/src/tc_scf/fock_three.irp.f
index e4348892..e9ad4ce5 100644
--- a/src/tc_scf/fock_three.irp.f
+++ b/src/tc_scf/fock_three.irp.f
@@ -16,6 +16,7 @@ BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)]
      fock_3_mat(j,i) = -contrib
     enddo
    enddo
+  else if(bi_ortho.and.three_body_h_tc)then
 !!  !$OMP END DO
 !!  !$OMP END PARALLEL
 !!  do i = 1, mo_num
diff --git a/src/tc_scf/fock_three_bi_ortho_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new.irp.f
new file mode 100644
index 00000000..004a2aa4
--- /dev/null
+++ b/src/tc_scf/fock_three_bi_ortho_new.irp.f
@@ -0,0 +1,177 @@
+BEGIN_PROVIDER [ double precision, fock_a_abb_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_a_abb_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,beta,beta contribution
+ END_DOC
+ fock_a_abb_3e_bi_orth = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_23_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_beta_num
+    do k = j+1, elec_beta_num
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
+      fock_a_abb_3e_bi_orth(a,i) += direct_int - exch_23_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_a_abb_3e_bi_orth = - fock_a_abb_3e_bi_orth
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_a_aba_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,beta contribution
+ END_DOC
+ fock_a_aba_3e_bi_orth = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_alpha_num ! a
+    do k = 1, elec_beta_num ! b
+                                                                       !   a b a   a b a
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 
+      fock_a_aba_3e_bi_orth(a,i) += direct_int - exch_13_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_a_aba_3e_bi_orth = - fock_a_aba_3e_bi_orth
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_a_aaa_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_a_aaa_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution
+ END_DOC
+ fock_a_aaa_3e_bi_orth = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_alpha_num
+    do k = j+1, elec_alpha_num
+      ! positive terms :: cycle contrib 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+      call  give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+      fock_a_aaa_3e_bi_orth(a,i) += direct_int + c_3_int + c_minus_3_int 
+      ! negative terms :: exchange contrib
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+      call  give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+      fock_a_aaa_3e_bi_orth(a,i) += - exch_13_int - exch_23_int  - exch_12_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_a_aaa_3e_bi_orth = - fock_a_aaa_3e_bi_orth
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! fock_a_tot_3e_bi_orth = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions 
+ END_DOC
+ fock_a_tot_3e_bi_orth = fock_a_abb_3e_bi_orth + fock_a_aba_3e_bi_orth + fock_a_aaa_3e_bi_orth 
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_baa_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_b_baa_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,alpha contribution
+ END_DOC
+ fock_b_baa_3e_bi_orth = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_23_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_alpha_num
+    do k = j+1, elec_alpha_num
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23
+      fock_b_baa_3e_bi_orth(a,i) += direct_int - exch_23_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_b_baa_3e_bi_orth = - fock_b_baa_3e_bi_orth
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_b_bab_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,beta contribution
+ END_DOC
+ fock_b_bab_3e_bi_orth = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_beta_num
+    do k = 1, elec_alpha_num
+      !                                                                    b a b   b a b
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13
+      fock_b_bab_3e_bi_orth(a,i) += direct_int - exch_13_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_b_bab_3e_bi_orth = - fock_b_bab_3e_bi_orth
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_bbb_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+! fock_b_bbb_3e_bi_orth(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution
+ END_DOC
+ fock_b_bbb_3e_bi_orth = 0.d0
+ integer :: i,a,j,k
+ double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+ do i = 1, mo_num
+  do a = 1, mo_num
+   
+   do j = 1, elec_beta_num
+    do k = j+1, elec_beta_num
+      ! positive terms :: cycle contrib 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j >
+      call  give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int)      ! < a k j | j i k >
+      call  give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i >
+      fock_b_bbb_3e_bi_orth(a,i) += direct_int + c_3_int + c_minus_3_int 
+      ! negative terms :: exchange contrib
+      call  give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 
+      call  give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23
+      call  give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12
+      fock_b_bbb_3e_bi_orth(a,i) += - exch_13_int - exch_23_int  - exch_12_int 
+    enddo
+   enddo
+
+  enddo
+ enddo
+ fock_b_bbb_3e_bi_orth = - fock_b_bbb_3e_bi_orth
+END_PROVIDER 
+
+BEGIN_PROVIDER [ double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)]
+ implicit none
+ BEGIN_DOC
+ ! fock_b_tot_3e_bi_orth = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions 
+ END_DOC
+ fock_b_tot_3e_bi_orth = fock_b_bbb_3e_bi_orth + fock_b_bab_3e_bi_orth + fock_b_baa_3e_bi_orth
+
+END_PROVIDER 
diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f
index 7e51caaf..3b9f6978 100644
--- a/src/tools/print_wf.irp.f
+++ b/src/tools/print_wf.irp.f
@@ -17,6 +17,7 @@ program print_wf
  ! psi_coef_sorted are the wave function stored in the |EZFIO| directory.
  read_wf = .True.
  touch read_wf
+ call write_wf
  call routine
 end
 
@@ -120,3 +121,19 @@ subroutine routine
  print*,'L2 norm of pert beta  = ',norm_mono_b_pert_2
 
 end
+
+subroutine write_wf
+ implicit none
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen                                                                                 
+ output=trim(ezfio_filename)//'.wf'
+ i_unit_output = getUnitAndOpen(output,'w')
+ integer :: i
+ print*,'Writing the sorted wf'
+ do i = 1, N_det
+  write(i_unit_output,*)i,psi_coef_sorted(i,1)/psi_coef_sorted(1,1)
+ enddo
+
+
+end
+