diff --git a/plugins/local/jastrow/EZFIO.cfg b/plugins/local/jastrow/EZFIO.cfg
index 23dde8ea..8fd2d05a 100644
--- a/plugins/local/jastrow/EZFIO.cfg
+++ b/plugins/local/jastrow/EZFIO.cfg
@@ -1,13 +1,13 @@
 
 [j2e_type]
 type: character*(32)
-doc: type of the 2e-Jastrow: [ None | Mu | Mur | Qmckl ]
+doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ]
 interface: ezfio,provider,ocaml
 default: Mu
 
 [j1e_type]
 type: character*(32)
-doc: type of the 1e-Jastrow: [ None | Gauss | Charge_Harmonizer ]
+doc: type of the 1e-Jastrow: [ None | Gauss | Charge_Harmonizer | Charge_Harmonizer_AO ]
 interface: ezfio,provider,ocaml
 default: None
 
@@ -151,3 +151,45 @@ interface: ezfio,provider,ocaml
 default: 1.0
 ezfio_name: nu_erf
 
+[jBH_size]
+type: integer
+doc: number of terms per atom in Boys-Handy-Jastrow
+interface: ezfio,provider,ocaml
+default: 1
+
+[jBH_c]
+type: double precision
+doc: coefficients of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_m]
+type: integer
+doc: powers of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_n]
+type: integer
+doc: powers of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_o]
+type: integer
+doc: powers of terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (jastrow.jBH_size,nuclei.nucl_num)
+
+[jBH_ee]
+type: double precision
+doc: parameters of e-e terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (nuclei.nucl_num)
+
+[jBH_en]
+type: double precision
+doc: parameters of e-n terms in Boys-Handy-Jastrow
+interface: ezfio
+size: (nuclei.nucl_num)
+
diff --git a/plugins/local/jastrow/bh_param.irp.f b/plugins/local/jastrow/bh_param.irp.f
new file mode 100644
index 00000000..1ed871bc
--- /dev/null
+++ b/plugins/local/jastrow/bh_param.irp.f
@@ -0,0 +1,252 @@
+
+ BEGIN_PROVIDER [double precision, jBH_ee,           (nucl_num)]
+&BEGIN_PROVIDER [double precision, jBH_en,           (nucl_num)]
+&BEGIN_PROVIDER [double precision, jBH_c , (jBH_size, nucl_num)]
+&BEGIN_PROVIDER [integer         , jBH_m , (jBH_size, nucl_num)]
+&BEGIN_PROVIDER [integer         , jBH_n , (jBH_size, nucl_num)]
+&BEGIN_PROVIDER [integer         , jBH_o , (jBH_size, nucl_num)]
+
+  BEGIN_DOC
+  !
+  ! parameters of Boys-Handy-Jastrow
+  !
+  END_DOC
+
+  implicit none
+  logical :: exists
+  integer :: i_nucl, p
+  integer :: ierr
+
+  PROVIDE ezfio_filename
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_ee(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    include 'mpif.h'
+    call MPI_BCAST(jBH_ee, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy e-e param with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_ee ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_ee(jBH_ee)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_ee, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_ee with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_ee = 1.d0
+    call ezfio_set_jastrow_jBH_ee(jBH_ee)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_en(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_en, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy e-n param with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_en ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_en(jBH_en)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_en, (nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if (ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_en with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_en = 1.d0
+    call ezfio_set_jastrow_jBH_en(jBH_en)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_c(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_c, (jBH_size*nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy coeff with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_c ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_c(jBH_c)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_c, (jBH_size*nucl_num), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_c with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_c = 0.d0
+    call ezfio_set_jastrow_jBH_c(jBH_c)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_m(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_m, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy m powers with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_m ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_m(jBH_m)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_m, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_m with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_m = 0
+    call ezfio_set_jastrow_jBH_m(jBH_m)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_n(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_n, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy n powers with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_n ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_n(jBH_n)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_n, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_n with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_n = 0
+    call ezfio_set_jastrow_jBH_n(jBH_n)
+  endif
+
+  ! ---
+
+  if(mpi_master) then
+    call ezfio_has_jastrow_jBH_o(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    call MPI_BCAST(jBH_o, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+    if(ierr /= MPI_SUCCESS) then
+      stop 'Unable to read Boys-Handy o powers with MPI'
+    endif
+  IRP_ENDIF
+
+  if(exists) then
+    if(mpi_master) then
+      write(6,'(A)') '.. >>>>> [ IO READ: jBH_o ] <<<<< ..'
+      call ezfio_get_jastrow_jBH_o(jBH_o)
+      IRP_IF MPI
+        call MPI_BCAST(jBH_o, (jBH_size*nucl_num), MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        if(ierr /= MPI_SUCCESS) then
+          stop 'Unable to read jBH_o with MPI'
+        endif
+      IRP_ENDIF
+    endif
+  else
+
+    jBH_o = 0
+    call ezfio_set_jastrow_jBH_o(jBH_o)
+  endif
+
+  ! ---
+
+  print *, ' parameters for Boys-Handy Jastrow'
+  print *, ' nb of terms per nucleus = ', jBH_size
+
+  do i_nucl = 1, nucl_num
+    print *, ' nucl    = ', nucl_label(i_nucl)
+    print *, ' ee-term = ', jBH_ee(i_nucl)
+    print *, ' en-term = ', jBH_en(i_nucl)
+    print *, '  m     n     o          c'
+    do p = 1, jBH_size
+      write(*,'(3(I4,2x), E15.7)') jBH_m(p,i_nucl), jBH_n(p,i_nucl), jBH_o(p,i_nucl), jBH_c(p,i_nucl)
+    enddo
+  enddo
+
+
+END_PROVIDER
+
+! ---
+
diff --git a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f
index 5777a44a..31ad5756 100644
--- a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f
+++ b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f
@@ -109,6 +109,14 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res)
 
     endif ! env_type
 
+  elseif(j2e_type .eq. "Boys_Handy") then
+
+    PROVIDE jBH_size jBH_en jBH_ee jBH_m jBH_n jBH_o jBH_c
+
+    if(env_type .eq. "None") then
+      call grad1_j12_r1_seq(r1, n_grid2, resx, resy, resz)
+    endif ! env_type
+
   else
 
     print *, ' Error in get_grad1_u12_withsq_r1_seq: Unknown Jastrow'
@@ -157,9 +165,13 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
   double precision, intent(out) :: gradz(n_grid2)
 
   integer                       :: jpoint
+  integer                       :: i_nucl, p, mpA, npA, opA
   double precision              :: r2(3)
   double precision              :: dx, dy, dz, r12, tmp
   double precision              :: mu_val, mu_tmp, mu_der(3)
+  double precision              :: rn(3), f1A, gard1_f1A(3), f2A, gard2_f2A(3), g12, gard1_g12(3)
+  double precision              :: tmp1, tmp2
+
 
   PROVIDE j2e_type
 
@@ -267,6 +279,57 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       gradz(jpoint) = tmp * dz
     enddo
 
+  elseif(j2e_type .eq. "Boys_Handy") then
+
+    do jpoint = 1, n_points_extra_final_grid ! r2
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      gradx(jpoint) = 0.d0 
+      grady(jpoint) = 0.d0 
+      gradz(jpoint) = 0.d0 
+      do i_nucl = 1, nucl_num 
+
+        rn(1) = nucl_coord(i_nucl,1)
+        rn(2) = nucl_coord(i_nucl,2)
+        rn(3) = nucl_coord(i_nucl,3)
+
+        call jBH_elem_fct_grad(jBH_en(i_nucl), r1, rn, f1A, gard1_f1A)
+        call jBH_elem_fct_grad(jBH_en(i_nucl), r2, rn, f2A, gard2_f2A)
+        call jBH_elem_fct_grad(jBH_ee(i_nucl), r1, r2, g12, gard1_g12)
+
+        do p = 1, jBH_size
+          mpA = jBH_m(p,i_nucl)
+          npA = jBH_n(p,i_nucl)
+          opA = jBH_o(p,i_nucl)
+          tmp = jBH_c(p,i_nucl)
+          if(mpA .eq. npA) then
+            tmp = tmp * 0.5d0
+          endif
+
+          tmp1 = 0.d0
+          if(mpA .gt. 0) then
+            tmp1 = tmp1 + dble(mpA) * f1A**dble(mpA-1) * f2A**dble(npA)
+          endif
+          if(npA .gt. 0) then
+            tmp1 = tmp1 + dble(npA) * f1A**dble(npA-1) * f2A**dble(mpA)
+          endif
+          tmp1 = tmp1 * g12**dble(opA)
+
+          tmp2 = 0.d0
+          if(opA .gt. 0) then
+            tmp2 = tmp2 + dble(opA) * g12**dble(opA-1) * (f1A**dble(mpA) * f2A**dble(npA) + f1A**dble(npA) * f2A**dble(mpA))
+          endif
+
+          gradx(jpoint) = gradx(jpoint) + tmp * (tmp1 * gard1_f1A(1) + tmp2 * gard1_g12(1))
+          grady(jpoint) = grady(jpoint) + tmp * (tmp1 * gard1_f1A(2) + tmp2 * gard1_g12(2))
+          gradz(jpoint) = gradz(jpoint) + tmp * (tmp1 * gard1_f1A(3) + tmp2 * gard1_g12(3))
+        enddo ! p
+      enddo ! i_nucl
+    enddo ! jpoint
+
   else
 
     print *, ' Error in grad1_j12_r1_seq: Unknown j2e_type = ', j2e_type
@@ -757,3 +820,34 @@ end
 
 ! ---
 
+subroutine jBH_elem_fct_grad(alpha, r1, r2, fct, gard1_fct)
+
+  implicit none
+  double precision, intent(in)  :: alpha, r1(3), r2(3)
+  double precision, intent(out) :: fct, gard1_fct(3)
+  double precision              :: dist, tmp1, tmp2
+
+  dist = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
+              + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
+              + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
+
+  tmp1 = 1.d0 / (1.d0 + alpha * dist)
+
+  fct = alpha * dist * tmp1
+
+  if(dist .lt. 1d-10) then
+    gard1_fct(1) = 0.d0
+    gard1_fct(2) = 0.d0
+    gard1_fct(3) = 0.d0
+  else
+    tmp2 = alpha * tmp1 * tmp1 / dist
+    gard1_fct(1) = tmp2 * (r1(1) - r2(1))
+    gard1_fct(2) = tmp2 * (r1(2) - r2(2))
+    gard1_fct(3) = tmp2 * (r1(3) - r2(3))
+  endif
+
+  return
+end 
+
+! ---
+
diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f
index 9901a853..59efc943 100644
--- a/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f
+++ b/plugins/local/tc_bi_ortho/slater_tc_opt.irp.f
@@ -8,8 +8,13 @@ subroutine provide_all_three_ints_bi_ortho()
   END_DOC
 
   implicit none
+  double precision :: t1, t2
+
   PROVIDE ao_two_e_integrals_in_map
 
+  print *, ' start provide_all_three_ints_bi_ortho'
+  call wall_time(t1)
+
   if(three_body_h_tc) then
 
     if(three_e_3_idx_term) then
@@ -32,6 +37,9 @@ subroutine provide_all_three_ints_bi_ortho()
 
   endif
 
+  call wall_time(t2)
+  print *, ' end provide_all_three_ints_bi_ortho after (min) = ', (t2-t1)/60.d0
+
   return
 end
 
@@ -83,8 +91,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
   integer,           intent(in) :: Nint
   integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
   double precision, intent(out) :: hmono, htwoe, hthree, htot
+
   integer                       :: degree
 
+  PROVIDE pure_three_body_h_tc
+
   hmono  = 0.d0
   htwoe  = 0.d0
   htot   = 0.d0
@@ -99,7 +110,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
     if(degree == 0) then
       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_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
     endif
@@ -111,7 +122,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
     if(degree == 0) then
       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_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
     else
@@ -149,16 +160,16 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
   double precision, intent(out) :: htot
   integer                       :: degree
 
-  htot   = 0.d0
+  htot = 0.d0
 
   call get_excitation_degree(key_i, key_j, degree, Nint)
   if(degree.gt.2) return
 
-  if(degree == 0)then
+  if(degree == 0) then
     call diag_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_i,htot)
-  else if (degree == 1)then
+  else if (degree == 1) then
     call single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint,key_j, key_i , htot)
-  else if(degree == 2)then
+  else if(degree == 2) then
     call double_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
   endif
 
diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f
index cc1a0603..78f9dc66 100644
--- a/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f
+++ b/plugins/local/tc_bi_ortho/slater_tc_opt_diag.irp.f
@@ -15,6 +15,8 @@
   implicit none
   double precision :: hmono, htwoe, htot, hthree
 
+  PROVIDE N_int
+  PROVIDE HF_bitmask
   PROVIDE mo_l_coef mo_r_coef
 
   call diag_htilde_mu_mat_bi_ortho_slow(N_int, HF_bitmask, hmono, htwoe, htot)
diff --git a/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f b/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f
index 81bf69f4..e57cb05c 100644
--- a/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f
+++ b/plugins/local/tc_bi_ortho/slater_tc_opt_single.irp.f
@@ -19,6 +19,7 @@ subroutine  single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
   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)
@@ -44,27 +45,28 @@ subroutine  single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe,
   call bitstring_to_list_ab(key_i, occ, Ne, Nint)
   call get_single_excitation(key_i, key_j, exc, phase, Nint)
   call decode_exc(exc, 1, h1, p1, h2, p2, s1, s2)
-  call get_single_excitation_from_fock_tc(key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot)
+  call get_single_excitation_from_fock_tc(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, hthree, htot)
 
 end
 
 ! ---
 
-subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot)
+subroutine get_single_excitation_from_fock_tc(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, hthree, htot)
 
   use bitmasks
 
   implicit none
+  integer,           intent(in)  :: Nint
   integer,           intent(in)  :: h, p, spin
   double precision,  intent(in)  :: phase
-  integer(bit_kind), intent(in)  :: key_i(N_int,2), key_j(N_int,2)
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
   double precision,  intent(out) :: hmono, htwoe, hthree, htot
 
-  integer(bit_kind)              :: differences(N_int,2)
-  integer(bit_kind)              :: hole(N_int,2)
-  integer(bit_kind)              :: partcl(N_int,2)
-  integer                        :: occ_hole(N_int*bit_kind_size,2)
-  integer                        :: occ_partcl(N_int*bit_kind_size,2)
+  integer(bit_kind)              :: differences(Nint,2)
+  integer(bit_kind)              :: hole(Nint,2)
+  integer(bit_kind)              :: partcl(Nint,2)
+  integer                        :: occ_hole(Nint*bit_kind_size,2)
+  integer                        :: occ_partcl(Nint*bit_kind_size,2)
   integer                        :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
   integer                        :: i0,i
   double precision               :: buffer_c(mo_num),buffer_x(mo_num)
@@ -74,7 +76,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h
     buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
   enddo
 
-  do i = 1, N_int
+  do i = 1, Nint
     differences(i,1) = xor(key_i(i,1), ref_closed_shell_bitmask(i,1))
     differences(i,2) = xor(key_i(i,2), ref_closed_shell_bitmask(i,2))
     hole       (i,1) = iand(differences(i,1), ref_closed_shell_bitmask(i,1))
@@ -83,8 +85,8 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h
     partcl     (i,2) = iand(differences(i,2), key_i(i,2))
   enddo
 
-  call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
-  call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
+  call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint)
+  call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
   hmono = mo_bi_ortho_tc_one_e(p,h)
   htwoe = fock_op_2_e_tc_closed_shell(p,h)
 
@@ -122,7 +124,7 @@ subroutine get_single_excitation_from_fock_tc(key_i, key_j, h, p, spin, phase, h
 
   hthree = 0.d0
   if (three_body_h_tc .and. elec_num.gt.2 .and. three_e_4_idx_term) then
-    call three_comp_fock_elem(key_i, h, p, spin, hthree)
+    call three_comp_fock_elem(Nint, key_i, h, p, spin, hthree)
   endif
 
   htwoe  = htwoe * phase
@@ -134,24 +136,27 @@ end
 
 ! ---
 
-subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
- implicit none
- integer,intent(in) :: h_fock,p_fock,ispin_fock
- integer(bit_kind), intent(in) :: key_i(N_int,2)
- 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)
+subroutine three_comp_fock_elem(Nint, key_i, h_fock, p_fock, ispin_fock, hthree)
 
+  implicit none
+  integer,           intent(in)  :: Nint
+  integer,           intent(in)  :: h_fock, p_fock, ispin_fock
+  integer(bit_kind), intent(in)  :: key_i(Nint,2)
+  double precision,  intent(out) :: hthree
+
+  integer                        :: nexc(2),i,ispin,na,nb
+  integer(bit_kind)              :: hole(Nint,2)
+  integer(bit_kind)              :: particle(Nint,2)
+  integer                        :: occ_hole(Nint*bit_kind_size,2)
+  integer                        :: occ_particle(Nint*bit_kind_size,2)
+  integer                        :: n_occ_ab_hole(2),n_occ_ab_particle(2)
+  integer(bit_kind)              :: det_tmp(Nint,2)
 
   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
+  do i = 1, Nint
     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))
@@ -161,13 +166,14 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
     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)
+  call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
   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)
+  call bitstring_to_list_ab(hole, occ_hole, tmp, Nint)
   ASSERT (tmp(1) == nexc(1)) ! Number of holes alpha
   ASSERT (tmp(2) == nexc(2)) ! Number of holes beta 
 
@@ -181,15 +187,18 @@ subroutine three_comp_fock_elem(key_i,h_fock,p_fock,ispin_fock,hthree)
   do ispin=1,2
     na = elec_num_tab(ispin)
     nb = elec_num_tab(iand(ispin,1)+1)
-    do i=1,nexc(ispin)
+    do i = 1, nexc(ispin)
       !DIR$ FORCEINLINE
-      call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb)
+      call fock_ac_tc_operator( occ_particle(i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb)
       !DIR$ FORCEINLINE
-      call fock_a_tc_operator ( occ_hole    (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, N_int,na,nb)
+      call fock_a_tc_operator ( occ_hole    (i,ispin), ispin, det_tmp, h_fock,p_fock, ispin_fock, hthree, Nint, na, nb)
     enddo
   enddo
+
 end
 
+! ---
+
 subroutine fock_ac_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,Nint,na,nb)
   use bitmasks
   implicit none
@@ -365,111 +374,118 @@ subroutine fock_a_tc_operator(iorb,ispin,key, h_fock,p_fock, ispin_fock,hthree,N
 
 end
 
+! ---
 
-BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num) ]
- implicit none
- BEGIN_DOC
-! Closed-shell part of the Fock operator for the TC operator
- END_DOC
- integer :: h0,p0,h,p,k0,k,i
- integer :: n_occ_ab(2)
- integer :: occ(N_int*bit_kind_size,2)
- integer :: n_occ_ab_virt(2)
- integer :: occ_virt(N_int*bit_kind_size,2)
- integer(bit_kind) :: key_test(N_int)
- integer(bit_kind) :: key_virt(N_int,2)
- double precision :: accu
+BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)]
 
- fock_op_2_e_tc_closed_shell = -1000.d0
- call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
- do i = 1, N_int
-  key_virt(i,1) = full_ijkl_bitmask(i)
-  key_virt(i,2) = full_ijkl_bitmask(i)
-  key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1))
-  key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2))
- enddo
- call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int)
- ! docc ---> virt single excitations
- do h0 = 1,  n_occ_ab(1)
-  h=occ(h0,1)
-  do p0 = 1, n_occ_ab_virt(1)
-   p = occ_virt(p0,1)
-   accu = 0.d0
-   do k0 = 1, n_occ_ab(1)
-    k = occ(k0,1)
-    accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
-   enddo
-   fock_op_2_e_tc_closed_shell(p,h) = accu 
+  BEGIN_DOC
+  ! Closed-shell part of the Fock operator for the TC operator
+  END_DOC
+ 
+  implicit none
+
+  PROVIDE N_int
+
+  integer           :: h0,p0,h,p,k0,k,i
+  integer           :: n_occ_ab(2)
+  integer           :: occ(N_int*bit_kind_size,2)
+  integer           :: n_occ_ab_virt(2)
+  integer           :: occ_virt(N_int*bit_kind_size,2)
+  integer(bit_kind) :: key_test(N_int)
+  integer(bit_kind) :: key_virt(N_int,2)
+  double precision  :: accu
+ 
+  fock_op_2_e_tc_closed_shell = -1000.d0
+  call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int)
+
+  do i = 1, N_int
+    key_virt(i,1) = full_ijkl_bitmask(i)
+    key_virt(i,2) = full_ijkl_bitmask(i)
+    key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1))
+    key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2))
   enddo
- enddo
-
- do h0 = 1, n_occ_ab_virt(1)
-  h = occ_virt(h0,1)
-  do p0 = 1,  n_occ_ab(1)
-   p=occ(p0,1)
-   accu = 0.d0
-   do k0 = 1, n_occ_ab(1)
-    k = occ(k0,1)
-    accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
-   enddo
-   fock_op_2_e_tc_closed_shell(p,h) = accu 
+  call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int)
+  ! docc ---> virt single excitations
+  do h0 = 1,  n_occ_ab(1)
+    h = occ(h0,1)
+    do p0 = 1, n_occ_ab_virt(1)
+      p = occ_virt(p0,1)
+      accu = 0.d0
+      do k0 = 1, n_occ_ab(1)
+        k = occ(k0,1)
+        accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+      enddo
+      fock_op_2_e_tc_closed_shell(p,h) = accu 
+    enddo
   enddo
- enddo
-
- ! virt ---> virt single excitations
- do h0 = 1,  n_occ_ab_virt(1)
-  h=occ_virt(h0,1)
-  do p0 = 1, n_occ_ab_virt(1)
-   p = occ_virt(p0,1)
-   accu = 0.d0
-   do k0 = 1, n_occ_ab(1)
-    k = occ(k0,1)
-    accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
-   enddo
-   fock_op_2_e_tc_closed_shell(p,h) = accu 
+ 
+  do h0 = 1, n_occ_ab_virt(1)
+    h = occ_virt(h0,1)
+    do p0 = 1,  n_occ_ab(1)
+      p = occ(p0,1)
+      accu = 0.d0
+      do k0 = 1, n_occ_ab(1)
+        k = occ(k0,1)
+        accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+      enddo
+      fock_op_2_e_tc_closed_shell(p,h) = accu 
+    enddo
   enddo
- enddo
-
- do h0 = 1, n_occ_ab_virt(1)
-  h = occ_virt(h0,1)
-  do p0 = 1,  n_occ_ab_virt(1)
-   p=occ_virt(p0,1)
-   accu = 0.d0
-   do k0 = 1, n_occ_ab(1)
-    k = occ(k0,1)
-    accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+ 
+  ! virt ---> virt single excitations
+  do h0 = 1,  n_occ_ab_virt(1)
+   h=occ_virt(h0,1)
+   do p0 = 1, n_occ_ab_virt(1)
+    p = occ_virt(p0,1)
+    accu = 0.d0
+    do k0 = 1, n_occ_ab(1)
+     k = occ(k0,1)
+     accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+    enddo
+    fock_op_2_e_tc_closed_shell(p,h) = accu 
    enddo
-   fock_op_2_e_tc_closed_shell(p,h) = accu 
   enddo
- enddo
-
-
- ! docc ---> docc single excitations
- do h0 = 1,  n_occ_ab(1)
-  h=occ(h0,1)
-  do p0 = 1, n_occ_ab(1)
-   p = occ(p0,1)
-   accu = 0.d0
-   do k0 = 1, n_occ_ab(1)
-    k = occ(k0,1)
-    accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+ 
+  do h0 = 1, n_occ_ab_virt(1)
+   h = occ_virt(h0,1)
+   do p0 = 1,  n_occ_ab_virt(1)
+    p=occ_virt(p0,1)
+    accu = 0.d0
+    do k0 = 1, n_occ_ab(1)
+     k = occ(k0,1)
+     accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+    enddo
+    fock_op_2_e_tc_closed_shell(p,h) = accu 
    enddo
-   fock_op_2_e_tc_closed_shell(p,h) = accu 
   enddo
- enddo
-
- do h0 = 1, n_occ_ab(1)
-  h = occ(h0,1)
-  do p0 = 1,  n_occ_ab(1)
-   p=occ(p0,1)
-   accu = 0.d0
-   do k0 = 1, n_occ_ab(1)
-    k = occ(k0,1)
-    accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+ 
+ 
+  ! docc ---> docc single excitations
+  do h0 = 1,  n_occ_ab(1)
+   h=occ(h0,1)
+   do p0 = 1, n_occ_ab(1)
+    p = occ(p0,1)
+    accu = 0.d0
+    do k0 = 1, n_occ_ab(1)
+     k = occ(k0,1)
+     accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+    enddo
+    fock_op_2_e_tc_closed_shell(p,h) = accu 
+   enddo
+  enddo
+ 
+  do h0 = 1, n_occ_ab(1)
+   h = occ(h0,1)
+   do p0 = 1,  n_occ_ab(1)
+    p=occ(p0,1)
+    accu = 0.d0
+    do k0 = 1, n_occ_ab(1)
+     k = occ(k0,1)
+     accu += 2.d0 * tc_2e_3idx_coulomb_integrals(k,p,h) - tc_2e_3idx_exchange_integrals(k,p,h)
+    enddo
+    fock_op_2_e_tc_closed_shell(p,h) = accu 
    enddo
-   fock_op_2_e_tc_closed_shell(p,h) = accu 
   enddo
- enddo
 
 ! do i = 1, mo_num
 !  write(*,'(100(F10.5,X))')fock_op_2_e_tc_closed_shell(:,i)
@@ -477,8 +493,10 @@ BEGIN_PROVIDER [double precision, fock_op_2_e_tc_closed_shell, (mo_num, mo_num)
 
 END_PROVIDER
 
+! ---
 
 subroutine  single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
+
   BEGIN_DOC
   ! <key_j |H_tilde | key_i> for single excitation ONLY FOR ONE- AND TWO-BODY TERMS 
   !!
@@ -492,8 +510,9 @@ subroutine  single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
   implicit none
   integer,           intent(in) :: Nint
   integer(bit_kind), intent(in) :: key_j(Nint,2), key_i(Nint,2)
-  double precision, intent(out) ::  htot
-  double precision :: hmono, htwoe
+  double precision, intent(out) :: htot
+
+  double precision              :: hmono, htwoe
   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)
@@ -517,75 +536,85 @@ subroutine  single_htilde_mu_mat_fock_bi_ortho_no_3e(Nint, key_j, key_i, htot)
 
   call get_single_excitation(key_i, key_j, exc, phase, Nint)
   call decode_exc(exc,1,h1,p1,h2,p2,s1,s2)
-  call get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h1,p1,s1,phase,hmono,htwoe,htot)
-end
-
-
-subroutine get_single_excitation_from_fock_tc_no_3e(key_i,key_j,h,p,spin,phase,hmono,htwoe,htot)
- use bitmasks
- implicit none
- integer,intent(in) :: h,p,spin
- double precision, intent(in)  :: phase
- integer(bit_kind), intent(in) :: key_i(N_int,2), key_j(N_int,2)
- double precision, intent(out) :: hmono,htwoe,htot
- integer(bit_kind) :: differences(N_int,2)
- integer(bit_kind) :: hole(N_int,2)
- integer(bit_kind) :: partcl(N_int,2)
- integer :: occ_hole(N_int*bit_kind_size,2)
- integer :: occ_partcl(N_int*bit_kind_size,2)
- integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
- integer :: i0,i
- double precision :: buffer_c(mo_num),buffer_x(mo_num)
- do i=1, mo_num
-   buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
-   buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
- enddo
- do i = 1, N_int
-  differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
-  differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
-  hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
-  hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
-  partcl(i,1) = iand(differences(i,1),key_i(i,1))
-  partcl(i,2) = iand(differences(i,2),key_i(i,2))
- enddo
- call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int)
- call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int)
- hmono = mo_bi_ortho_tc_one_e(p,h)
- htwoe = fock_op_2_e_tc_closed_shell(p,h)
- ! holes :: direct terms
- do i0 = 1, n_occ_ab_hole(1)
-  i = occ_hole(i0,1)
-  htwoe -= buffer_c(i)
- enddo
- do i0 = 1, n_occ_ab_hole(2)
-  i = occ_hole(i0,2)
-  htwoe -= buffer_c(i)
- enddo
-
- ! holes :: exchange terms
- do i0 = 1, n_occ_ab_hole(spin)
-  i = occ_hole(i0,spin)
-  htwoe += buffer_x(i)
- enddo
-
- ! particles :: direct terms
- do i0 = 1, n_occ_ab_partcl(1)
-  i = occ_partcl(i0,1)
-  htwoe += buffer_c(i)
- enddo
- do i0 = 1, n_occ_ab_partcl(2)
-  i = occ_partcl(i0,2)
-  htwoe += buffer_c(i)
- enddo
-
- ! particles :: exchange terms
- do i0 = 1, n_occ_ab_partcl(spin)
-  i = occ_partcl(i0,spin)
-  htwoe -= buffer_x(i)
- enddo
- htwoe = htwoe * phase
- hmono = hmono * phase
- htot  = htwoe + hmono 
+  call get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h1, p1, s1, phase, hmono, htwoe, htot)
+
+end
+
+! ---
+
+subroutine get_single_excitation_from_fock_tc_no_3e(Nint, key_i, key_j, h, p, spin, phase, hmono, htwoe, htot)
+
+  use bitmasks
+
+  implicit none
+  integer,           intent(in) :: Nint
+  integer,           intent(in) :: h, p, spin
+  double precision,  intent(in) :: phase
+  integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out) :: hmono,htwoe,htot
+
+  integer(bit_kind)             :: differences(Nint,2)
+  integer(bit_kind)             :: hole(Nint,2)
+  integer(bit_kind)             :: partcl(Nint,2)
+  integer                       :: occ_hole(Nint*bit_kind_size,2)
+  integer                       :: occ_partcl(Nint*bit_kind_size,2)
+  integer                       :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
+  integer                       :: i0,i
+  double precision              :: buffer_c(mo_num), buffer_x(mo_num)
+
+  do i = 1, mo_num
+    buffer_c(i) = tc_2e_3idx_coulomb_integrals(i,p,h)
+    buffer_x(i) = tc_2e_3idx_exchange_integrals(i,p,h)
+  enddo
+
+  do i = 1, Nint
+   differences(i,1) = xor(key_i(i,1),ref_closed_shell_bitmask(i,1))
+   differences(i,2) = xor(key_i(i,2),ref_closed_shell_bitmask(i,2))
+   hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1))
+   hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2))
+   partcl(i,1) = iand(differences(i,1),key_i(i,1))
+   partcl(i,2) = iand(differences(i,2),key_i(i,2))
+  enddo
+
+  call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, Nint)
+  call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, Nint)
+  hmono = mo_bi_ortho_tc_one_e(p,h)
+  htwoe = fock_op_2_e_tc_closed_shell(p,h)
+
+  ! holes :: direct terms
+  do i0 = 1, n_occ_ab_hole(1)
+    i = occ_hole(i0,1)
+    htwoe -= buffer_c(i)
+  enddo
+  do i0 = 1, n_occ_ab_hole(2)
+    i = occ_hole(i0,2)
+    htwoe -= buffer_c(i)
+  enddo
+ 
+  ! holes :: exchange terms
+  do i0 = 1, n_occ_ab_hole(spin)
+    i = occ_hole(i0,spin)
+    htwoe += buffer_x(i)
+  enddo
+ 
+  ! particles :: direct terms
+  do i0 = 1, n_occ_ab_partcl(1)
+    i = occ_partcl(i0,1)
+    htwoe += buffer_c(i)
+  enddo
+  do i0 = 1, n_occ_ab_partcl(2)
+    i = occ_partcl(i0,2)
+    htwoe += buffer_c(i)
+  enddo
+ 
+  ! particles :: exchange terms
+  do i0 = 1, n_occ_ab_partcl(spin)
+    i = occ_partcl(i0,spin)
+    htwoe -= buffer_x(i)
+  enddo
+  htwoe = htwoe * phase
+  hmono = hmono * phase
+  htot  = htwoe + hmono 
 
 end
 
diff --git a/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f b/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f
index 64982ab6..398e96db 100644
--- a/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f
+++ b/plugins/local/tc_bi_ortho/tc_bi_ortho.irp.f
@@ -7,6 +7,10 @@ program tc_bi_ortho
   !
   END_DOC
 
+  implicit none
+
+  PROVIDE N_int
+
   my_grid_becke = .True.
   PROVIDE tc_grid1_a tc_grid1_r
   my_n_pt_r_grid = tc_grid1_r
@@ -66,6 +70,15 @@ subroutine routine_diag()
   ! provide overlap_bi_ortho
   ! provide htilde_matrix_elmt_bi_ortho
 
+  if(noL_standard) then
+    PROVIDE noL_0e
+    PROVIDE noL_1e
+    PROVIDE noL_2e
+  endif
+
+  PROVIDE htilde_matrix_elmt_bi_ortho
+  return
+
   if(N_states .eq. 1) then
 
     print*,'eigval_right_tc_bi_orth   = ',eigval_right_tc_bi_orth(1)
diff --git a/plugins/local/tc_bi_ortho/tc_hmat.irp.f b/plugins/local/tc_bi_ortho/tc_hmat.irp.f
index 88652caa..abec410d 100644
--- a/plugins/local/tc_bi_ortho/tc_hmat.irp.f
+++ b/plugins/local/tc_bi_ortho/tc_hmat.irp.f
@@ -13,16 +13,34 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
  
   implicit none
   integer          :: i, j
+  double precision :: t1, t2
   double precision :: htot
+ 
+  PROVIDE N_int
+  PROVIDE psi_det
+  PROVIDE three_e_3_idx_term
 
-  call provide_all_three_ints_bi_ortho
+  if(noL_standard) then
+    PROVIDE noL_0e
+    PROVIDE noL_1e
+    PROVIDE noL_2e
+  endif
+
+  print *, ' PROVIDING htilde_matrix_elmt_bi_ortho ...'
+  call wall_time(t1)
+
+  call provide_all_three_ints_bi_ortho()
 
   i = 1
   j = 1
   call htilde_mu_mat_opt_bi_ortho_tot(psi_det(1,1,j), psi_det(1,1,i), N_int, htot)
 
-  !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j, htot) &
+
+  !$OMP PARALLEL            &
+  !$OMP DEFAULT(NONE)       &
+  !$OMP PRIVATE(i, j, htot) &
   !$OMP SHARED (N_det, psi_det, N_int, htilde_matrix_elmt_bi_ortho)
+  !$OMP DO 
   do i = 1, N_det
     do j = 1, N_det
       ! < J |Htilde | I >
@@ -31,7 +49,11 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
       htilde_matrix_elmt_bi_ortho(j,i) = htot
     enddo
   enddo
-  !$OMP END PARALLEL DO
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call wall_time(t2)
+  print *, ' wall time for htilde_matrix_elmt_bi_ortho (min) =', (t2-t1)/60.d0 
 
 END_PROVIDER 
 
diff --git a/plugins/local/tc_scf/EZFIO.cfg b/plugins/local/tc_scf/EZFIO.cfg
index 313d6f2b..3dfa9a71 100644
--- a/plugins/local/tc_scf/EZFIO.cfg
+++ b/plugins/local/tc_scf/EZFIO.cfg
@@ -2,3 +2,10 @@
 type: Threshold
 doc: Energy bi-tc HF
 interface: ezfio
+
+[converged_tcscf]
+type: logical
+doc: If |true|, tc-scf has converged
+interface: ezfio,provider,ocaml
+default: False
+
diff --git a/plugins/local/tc_scf/rh_tcscf_diis.irp.f b/plugins/local/tc_scf/rh_tcscf_diis.irp.f
index 12678500..431b6e08 100644
--- a/plugins/local/tc_scf/rh_tcscf_diis.irp.f
+++ b/plugins/local/tc_scf/rh_tcscf_diis.irp.f
@@ -22,6 +22,9 @@ subroutine rh_tcscf_diis()
 
   logical, external             :: qp_stop
 
+  PROVIDE level_shift_TCSCF
+  PROVIDE mo_l_coef mo_r_coef
+
   it          = 0
   e_save      = 0.d0
   dim_DIIS    = 0
@@ -41,19 +44,6 @@ subroutine rh_tcscf_diis()
 
   ! ---
 
-  PROVIDE level_shift_TCSCF
-  PROVIDE mo_l_coef mo_r_coef
-
-  !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
-  !  '====', '================', '================', '================', '================', '================' &
-  !        , '================', '================', '================', '====', '========'
-  !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
-  !  ' it ', '  SCF TC Energy ', '      E(1e)     ', '      E(2e)     ', '      E(3e)     ', '   energy diff  ' &
-  !        , '    gradient    ', '    DIIS error  ', '  level shift   ', 'DIIS', '  WT (m)'
-  !write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')              &
-  !  '====', '================', '================', '================', '================', '================' &
-  !        , '================', '================', '================', '====', '========'
-
   write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)')                      &
     '====', '================', '================', '================', '================', '================' &
           , '================', '================', '====', '========'
@@ -81,8 +71,6 @@ subroutine rh_tcscf_diis()
   er_save = er_DIIS
 
   call wall_time(t1)
-  !write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
-  !  it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
   write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)')  &
     it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
 
@@ -91,6 +79,8 @@ subroutine rh_tcscf_diis()
   PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
 
   converged = .false.
+  call ezfio_set_tc_scf_converged_tcscf(converged)
+
   !do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. dsqrt(thresh_tcscf)))
   do while(.not. converged)
 
@@ -253,8 +243,9 @@ subroutine rh_tcscf_diis()
     endif
     call lock_io
 
-    if (converged) then
+    if(converged) then
       write(json_unit, json_true_fmtx) 'converged'
+      call ezfio_set_tc_scf_converged_tcscf(converged)
     else
       write(json_unit, json_false_fmtx) 'converged'
     endif