diff --git a/plugins/local/bi_ort_ints/bi_ort_ints.irp.f b/plugins/local/bi_ort_ints/bi_ort_ints.irp.f
index b691e47e..0398a18f 100644
--- a/plugins/local/bi_ort_ints/bi_ort_ints.irp.f
+++ b/plugins/local/bi_ort_ints/bi_ort_ints.irp.f
@@ -24,6 +24,7 @@ program bi_ort_ints
   !call test_5idx
 
   call test_mos_in_r()
+  call test_int2_grad1_u12_bimo_t()
 
 end
 
@@ -486,7 +487,7 @@ subroutine test_mos_in_r()
   PROVIDE mos_l_in_r_array_transp_old mos_r_in_r_array_transp_old
   PROVIDE mos_l_in_r_array_transp mos_r_in_r_array_transp
 
-  acc_thr = 1d-12
+  acc_thr = 1d-13
 
   err_tot = 0.d0
   nrm_tot = 0.d0
@@ -527,3 +528,41 @@ end
 
 ! ---
 
+subroutine test_int2_grad1_u12_bimo_t()
+
+  implicit none
+  integer          :: i, j, ipoint, m
+  double precision :: err_tot, nrm_tot, err_loc, acc_thr
+
+  PROVIDE int2_grad1_u12_bimo_t_old
+  PROVIDE int2_grad1_u12_bimo_t
+
+  acc_thr = 1d-13
+
+  err_tot = 0.d0
+  nrm_tot = 0.d0
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do m = 1, 3
+        do ipoint = 1, n_points_final_grid
+          err_loc = dabs(int2_grad1_u12_bimo_t_old(ipoint,m,j,i) - int2_grad1_u12_bimo_t(ipoint,m,j,i))
+          if(err_loc > acc_thr) then
+            print*, " error on", ipoint, m, j, i
+            print*, " old res", int2_grad1_u12_bimo_t_old(ipoint,m,j,i)
+            print*, " new res", int2_grad1_u12_bimo_t    (ipoint,m,j,i)
+            stop
+          endif
+          err_tot = err_tot + err_loc
+          nrm_tot = nrm_tot + dabs(int2_grad1_u12_bimo_t_old(ipoint,m,j,i))
+        enddo
+      enddo
+    enddo
+  enddo
+  print *, ' absolute accuracy on int2_grad1_u12_bimo_t (%) =', 100.d0 * err_tot / nrm_tot
+
+  return
+end
+
+! ---
+
+
diff --git a/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f b/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f
index 77e4cb9b..1fd5d666 100644
--- a/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f
+++ b/plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f
@@ -1,360 +1,54 @@
 
 ! ---
 
-! TODO :: optimization : transform into a DGEMM
-
-BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, n_points_final_grid)]
-
-  BEGIN_DOC
-  !
-  ! mo_v_ki_bi_ortho_erf_rk_cst_mu(k,i,ip) = int dr chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1 )/(2|r - R_ip|) on the BI-ORTHO MO basis 
-  ! 
-  ! where phi_k(r) is a LEFT MOs and phi_i(r) is a RIGHT MO
-  !
-  ! R_ip = the "ip"-th point of the DFT Grid
-  !
-  END_DOC
-
-  implicit none
-  integer :: ipoint
- !$OMP PARALLEL         &
- !$OMP DEFAULT (NONE)   &
- !$OMP PRIVATE (ipoint) & 
- !$OMP SHARED (n_points_final_grid,v_ij_erf_rk_cst_mu,mo_v_ki_bi_ortho_erf_rk_cst_mu)
- !$OMP DO SCHEDULE (dynamic)
-  do ipoint = 1, n_points_final_grid
-    call ao_to_mo_bi_ortho( v_ij_erf_rk_cst_mu            (1,1,ipoint), size(v_ij_erf_rk_cst_mu,             1) &
-                          , mo_v_ki_bi_ortho_erf_rk_cst_mu(1,1,ipoint), size(mo_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
-  enddo
- !$OMP END DO
- !$OMP END PARALLEL
-
-  mo_v_ki_bi_ortho_erf_rk_cst_mu = mo_v_ki_bi_ortho_erf_rk_cst_mu * 0.5d0
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the BI-ORTHO MO basis
-  !
-  END_DOC
-
-  implicit none
-  integer :: ipoint, i, j
-
-  do i = 1, mo_num
-    do j = 1, mo_num
-      do ipoint = 1, n_points_final_grid
-        mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ki_bi_ortho_erf_rk_cst_mu(j,i,ipoint)
-      enddo
-    enddo
-  enddo
-
-  !FREE mo_v_ki_bi_ortho_erf_rk_cst_mu
-
-END_PROVIDER 
-
-! ---
-
-! TODO :: optimization : transform into a DGEMM
-
-BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, 3, n_points_final_grid)]
-
-  BEGIN_DOC
-  !
-  ! mo_x_v_ki_bi_ortho_erf_rk_cst_mu(k,i,m,ip) = int dr x(m) * chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1)/2|r - R_ip| on the BI-ORTHO MO basis 
-  !
-  ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => x(m) = x, m=2 => x(m) = y, m=3 => x(m) = z,
-  !
-  ! R_ip = the "ip"-th point of the DFT Grid
-  !
-  END_DOC
-
-  implicit none
-  integer :: ipoint
-
- !$OMP PARALLEL         &
- !$OMP DEFAULT (NONE)   &
- !$OMP PRIVATE (ipoint) & 
- !$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ki_bi_ortho_erf_rk_cst_mu)
- !$OMP DO SCHEDULE (dynamic)
-  do ipoint = 1, n_points_final_grid
-
-    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,1,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
-                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,1,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
-    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,2,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
-                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,2,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
-    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,3,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
-                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,3,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
-
-  enddo
- !$OMP END DO
- !$OMP END PARALLEL
-
-  mo_x_v_ki_bi_ortho_erf_rk_cst_mu = 0.5d0 * mo_x_v_ki_bi_ortho_erf_rk_cst_mu
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, n_points_final_grid)]
-
-  implicit none
-  integer          :: i, j, ipoint
-  double precision :: wall0, wall1
-
-  !print *, ' providing int2_grad1_u12_ao_transp ...'
-  !call wall_time(wall0)
-
-  if(test_cycle_tc) then
-
-    PROVIDE int2_grad1_u12_ao_test
-
-    do ipoint = 1, n_points_final_grid
-      do i = 1, ao_num
-        do j = 1, ao_num
-          int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1)
-          int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2)
-          int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3)
-        enddo
-      enddo
-    enddo
-
-    FREE int2_grad1_u12_ao_test
-
-  else
-
-    PROVIDE int2_grad1_u12_ao
-
-    do ipoint = 1, n_points_final_grid
-      do i = 1, ao_num
-        do j = 1, ao_num
-          int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1)
-          int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2)
-          int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3)
-        enddo
-      enddo
-    enddo
-
-  endif
-
-  !call wall_time(wall1)
-  !print *, ' wall time for int2_grad1_u12_ao_transp (min) = ', (wall1 - wall0) / 60.d0
-  !call print_memory_usage()
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)]
-
-  implicit none
-  integer          :: ipoint
-  double precision :: wall0, wall1
-
-  PROVIDE mo_l_coef mo_r_coef
-  PROVIDE int2_grad1_u12_ao_transp
-
-  !print *, ' providing int2_grad1_u12_bimo_transp ...'
-  !call wall_time(wall0)
-
-  !$OMP PARALLEL         &
-  !$OMP DEFAULT (NONE)   &
-  !$OMP PRIVATE (ipoint) & 
-  !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp)
-  !$OMP DO SCHEDULE (dynamic)
-  do ipoint = 1, n_points_final_grid
-    call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp  (1,1,1,ipoint), size(int2_grad1_u12_ao_transp  , 1) &
-                          , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
-    call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp  (1,1,2,ipoint), size(int2_grad1_u12_ao_transp  , 1) &
-                          , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
-    call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp  (1,1,3,ipoint), size(int2_grad1_u12_ao_transp  , 1) &
-                          , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
-  enddo
-  !$OMP END DO
-  !$OMP END PARALLEL
-
-  !FREE int2_grad1_u12_ao_transp
-
-  !call wall_time(wall1)
-  !print *, ' wall time for int2_grad1_u12_bimo_transp (min) =', (wall1 - wall0) / 60.d0
-  !call print_memory_usage()
-
-END_PROVIDER 
-
-! ---
-
 BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)]
 
   implicit none
-  integer          :: i, j, ipoint
-  double precision :: wall0, wall1
-
-  !call wall_time(wall0)
-  !print *, ' providing int2_grad1_u12_bimo_t ...'
+  integer                       :: i, j, ipoint
+  double precision              :: tt1, tt2
+  double precision, allocatable :: tmp(:,:,:,:)
 
   PROVIDE mo_l_coef mo_r_coef
-  PROVIDE int2_grad1_u12_bimo_transp
+  PROVIDE int2_grad1_u12_ao
 
+  call wall_time(tt1)
+
+  allocate(tmp(mo_num,mo_num,n_points_final_grid,3))
+
+  !$OMP PARALLEL         &
+  !$OMP DEFAULT (NONE)   &
+  !$OMP PRIVATE (ipoint) & 
+  !$OMP SHARED (ao_num, mo_num, n_points_final_grid, int2_grad1_u12_ao, tmp)
+  !$OMP DO SCHEDULE (dynamic)
+  do ipoint = 1, n_points_final_grid
+    call ao_to_mo_bi_ortho(int2_grad1_u12_ao(1,1,ipoint,1), ao_num, tmp (1,1,ipoint,1), mo_num)
+    call ao_to_mo_bi_ortho(int2_grad1_u12_ao(1,1,ipoint,2), ao_num, tmp (1,1,ipoint,2), mo_num)
+    call ao_to_mo_bi_ortho(int2_grad1_u12_ao(1,1,ipoint,3), ao_num, tmp (1,1,ipoint,3), mo_num)
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  !$OMP PARALLEL               &
+  !$OMP DEFAULT (NONE)         &
+  !$OMP PRIVATE (i, j, ipoint) & 
+  !$OMP SHARED (mo_num, n_points_final_grid, tmp, int2_grad1_u12_bimo_t)
+  !$OMP DO COLLAPSE(2) SCHEDULE (dynamic)
   do ipoint = 1, n_points_final_grid
     do i = 1, mo_num
       do j = 1, mo_num
-        int2_grad1_u12_bimo_t(ipoint,1,j,i) = int2_grad1_u12_bimo_transp(j,i,1,ipoint)
-        int2_grad1_u12_bimo_t(ipoint,2,j,i) = int2_grad1_u12_bimo_transp(j,i,2,ipoint)
-        int2_grad1_u12_bimo_t(ipoint,3,j,i) = int2_grad1_u12_bimo_transp(j,i,3,ipoint)
+        int2_grad1_u12_bimo_t(ipoint,1,j,i) = tmp(j,i,ipoint,1)
+        int2_grad1_u12_bimo_t(ipoint,2,j,i) = tmp(j,i,ipoint,2)
+        int2_grad1_u12_bimo_t(ipoint,3,j,i) = tmp(j,i,ipoint,3)
       enddo                                  
     enddo
   enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
 
-  FREE int2_grad1_u12_bimo_transp
+  deallocate(tmp)
 
-  !call wall_time(wall1)
-  !print *, ' wall time for int2_grad1_u12_bimo_t (min) =', (wall1 - wall0) / 60.d0
-  !call print_memory_usage()
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)]
-
-  implicit none
-  integer          :: i, j, ipoint
-  double precision :: wall0, wall1
-
-  !call wall_time(wall0)
-  !print *, ' providing int2_grad1_u12_ao_t ...'
-
-  PROVIDE int2_grad1_u12_ao
-
-  do ipoint = 1, n_points_final_grid
-    do i = 1, ao_num
-      do j = 1, ao_num
-        int2_grad1_u12_ao_t(ipoint,1,j,i) = int2_grad1_u12_ao(j,i,ipoint,1)
-        int2_grad1_u12_ao_t(ipoint,2,j,i) = int2_grad1_u12_ao(j,i,ipoint,2)
-        int2_grad1_u12_ao_t(ipoint,3,j,i) = int2_grad1_u12_ao(j,i,ipoint,3)
-      enddo                                  
-    enddo
-  enddo
-
-  !call wall_time(wall1)
-  !print *, ' wall time for int2_grad1_u12_ao_t (min) =', (wall1 - wall0) / 60.d0
-  !call print_memory_usage()
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, 3, mo_num, mo_num)]
-
-  implicit none
-  integer :: i, j, ipoint
-
-  do i = 1, mo_num
-    do j = 1, mo_num
-      do ipoint = 1, n_points_final_grid
-        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,1,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,1,ipoint)
-        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,2,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,2,ipoint)
-        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,3,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,3,ipoint)
-      enddo
-    enddo
-  enddo
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, 3, mo_num, mo_num)]
-
-  BEGIN_DOC
-  !
-  ! x_W_ki_bi_ortho_erf_rk(ip,m,k,i) = \int dr chi_k(r) \frac{(1 - erf(mu |r-R_ip|))}{2|r-R_ip|} (x(m)-R_ip(m)) phi_i(r) ON THE BI-ORTHO MO BASIS 
-  !
-  ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z,
-  !
-  ! R_ip = the "ip"-th point of the DFT Grid
-  END_DOC
- 
-  implicit none
-  include 'constants.include.F'
- 
-  integer          :: ipoint, m, i, k
-  double precision :: xyz
-  double precision :: wall0, wall1
- 
-  !print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
-  !call wall_time(wall0)
-
- !$OMP PARALLEL                   &
- !$OMP DEFAULT (NONE)             &
- !$OMP PRIVATE (ipoint,m,i,k,xyz) & 
- !$OMP SHARED (x_W_ki_bi_ortho_erf_rk,n_points_final_grid,mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_num,final_grid_points) 
- !$OMP DO SCHEDULE (dynamic)
-  do i = 1, mo_num
-    do k = 1, mo_num
-      do m = 1, 3
-        do ipoint = 1, n_points_final_grid
-          xyz = final_grid_points(m,ipoint)
-          x_W_ki_bi_ortho_erf_rk(ipoint,m,k,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,m,k,i) - xyz * mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,k,i)
-        enddo
-      enddo
-    enddo
-  enddo
-
- !$OMP END DO
- !$OMP END PARALLEL
-
- ! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu_transp 
- ! FREE mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp
-
-  !call wall_time(wall1)
-  !print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
-
-END_PROVIDER 
-
-! ---
-
-BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_grid, 3, mo_num)]
-  BEGIN_DOC
-  ! x_W_ki_bi_ortho_erf_rk_diag(ip,m,i) = \int dr chi_i(r) (1 - erf(mu |r-R_ip|)) (x(m)-X(m)_ip) phi_i(r) ON THE BI-ORTHO MO BASIS 
-!
-! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z,
-!
-! R_ip = the "ip"-th point of the DFT Grid
-  END_DOC
-
-  implicit none
-  include 'constants.include.F'
- 
-  integer          :: ipoint, m, i
-  double precision :: xyz
-  double precision :: wall0, wall1
- 
-  !print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...'
-  !call wall_time(wall0)
-
- !$OMP PARALLEL                 &
- !$OMP DEFAULT (NONE)           &
- !$OMP PRIVATE (ipoint,m,i,xyz) & 
- !$OMP SHARED (x_W_ki_bi_ortho_erf_rk_diag,n_points_final_grid,mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_num,final_grid_points) 
- !$OMP DO SCHEDULE (dynamic)
-  do i = 1, mo_num
-    do m = 1, 3
-      do ipoint = 1, n_points_final_grid
-        xyz = final_grid_points(m,ipoint)
-        x_W_ki_bi_ortho_erf_rk_diag(ipoint,m,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,m,i,i) - xyz * mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,i,i)
-      enddo
-    enddo
-  enddo
-
- !$OMP END DO
- !$OMP END PARALLEL
-
-  !call wall_time(wall1)
-  !print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0
+  call wall_time(tt2)
+  write(*,"(A,2X,F15.7)") ' wall time for int2_grad1_u12_bimo_t (sec) = ', (tt2 - tt1)
 
 END_PROVIDER 
 
diff --git a/plugins/local/bi_ort_ints/semi_num_ints_mo_old.irp.f b/plugins/local/bi_ort_ints/semi_num_ints_mo_old.irp.f
new file mode 100644
index 00000000..c2b9ad6d
--- /dev/null
+++ b/plugins/local/bi_ort_ints/semi_num_ints_mo_old.irp.f
@@ -0,0 +1,362 @@
+
+! ---
+
+! TODO :: optimization : transform into a DGEMM
+
+BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! mo_v_ki_bi_ortho_erf_rk_cst_mu(k,i,ip) = int dr chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1 )/(2|r - R_ip|) on the BI-ORTHO MO basis 
+  ! 
+  ! where phi_k(r) is a LEFT MOs and phi_i(r) is a RIGHT MO
+  !
+  ! R_ip = the "ip"-th point of the DFT Grid
+  !
+  END_DOC
+
+  implicit none
+  integer :: ipoint
+ !$OMP PARALLEL         &
+ !$OMP DEFAULT (NONE)   &
+ !$OMP PRIVATE (ipoint) & 
+ !$OMP SHARED (n_points_final_grid,v_ij_erf_rk_cst_mu,mo_v_ki_bi_ortho_erf_rk_cst_mu)
+ !$OMP DO SCHEDULE (dynamic)
+  do ipoint = 1, n_points_final_grid
+    call ao_to_mo_bi_ortho( v_ij_erf_rk_cst_mu            (1,1,ipoint), size(v_ij_erf_rk_cst_mu,             1) &
+                          , mo_v_ki_bi_ortho_erf_rk_cst_mu(1,1,ipoint), size(mo_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  mo_v_ki_bi_ortho_erf_rk_cst_mu = mo_v_ki_bi_ortho_erf_rk_cst_mu * 0.5d0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, mo_num, mo_num)]
+
+  BEGIN_DOC
+  !
+  ! int dr phi_i(r) phi_j(r) (erf(mu(R) |r - R|) - 1)/(2|r - R|) on the BI-ORTHO MO basis
+  !
+  END_DOC
+
+  implicit none
+  integer :: ipoint, i, j
+
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do ipoint = 1, n_points_final_grid
+        mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,j,i) = mo_v_ki_bi_ortho_erf_rk_cst_mu(j,i,ipoint)
+      enddo
+    enddo
+  enddo
+
+  !FREE mo_v_ki_bi_ortho_erf_rk_cst_mu
+
+END_PROVIDER 
+
+! ---
+
+! TODO :: optimization : transform into a DGEMM
+
+BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu, (mo_num, mo_num, 3, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! mo_x_v_ki_bi_ortho_erf_rk_cst_mu(k,i,m,ip) = int dr x(m) * chi_k(r) phi_i(r) (erf(mu |r - R_ip|) - 1)/2|r - R_ip| on the BI-ORTHO MO basis 
+  !
+  ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => x(m) = x, m=2 => x(m) = y, m=3 => x(m) = z,
+  !
+  ! R_ip = the "ip"-th point of the DFT Grid
+  !
+  END_DOC
+
+  implicit none
+  integer :: ipoint
+
+ !$OMP PARALLEL         &
+ !$OMP DEFAULT (NONE)   &
+ !$OMP PRIVATE (ipoint) & 
+ !$OMP SHARED (n_points_final_grid,x_v_ij_erf_rk_cst_mu_transp,mo_x_v_ki_bi_ortho_erf_rk_cst_mu)
+ !$OMP DO SCHEDULE (dynamic)
+  do ipoint = 1, n_points_final_grid
+
+    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,1,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
+                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,1,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,2,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
+                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,2,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+    call ao_to_mo_bi_ortho( x_v_ij_erf_rk_cst_mu_transp     (1,1,3,ipoint), size(x_v_ij_erf_rk_cst_mu_transp,      1) &
+                          , mo_x_v_ki_bi_ortho_erf_rk_cst_mu(1,1,3,ipoint), size(mo_x_v_ki_bi_ortho_erf_rk_cst_mu, 1) )
+
+  enddo
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  mo_x_v_ki_bi_ortho_erf_rk_cst_mu = 0.5d0 * mo_x_v_ki_bi_ortho_erf_rk_cst_mu
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3, n_points_final_grid)]
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  !print *, ' providing int2_grad1_u12_ao_transp ...'
+  !call wall_time(wall0)
+
+  if(test_cycle_tc) then
+
+    PROVIDE int2_grad1_u12_ao_test
+
+    do ipoint = 1, n_points_final_grid
+      do i = 1, ao_num
+        do j = 1, ao_num
+          int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,1)
+          int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,2)
+          int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao_test(j,i,ipoint,3)
+        enddo
+      enddo
+    enddo
+
+    FREE int2_grad1_u12_ao_test
+
+  else
+
+    PROVIDE int2_grad1_u12_ao
+
+    do ipoint = 1, n_points_final_grid
+      do i = 1, ao_num
+        do j = 1, ao_num
+          int2_grad1_u12_ao_transp(j,i,1,ipoint) = int2_grad1_u12_ao(j,i,ipoint,1)
+          int2_grad1_u12_ao_transp(j,i,2,ipoint) = int2_grad1_u12_ao(j,i,ipoint,2)
+          int2_grad1_u12_ao_transp(j,i,3,ipoint) = int2_grad1_u12_ao(j,i,ipoint,3)
+        enddo
+      enddo
+    enddo
+
+  endif
+
+  !call wall_time(wall1)
+  !print *, ' wall time for int2_grad1_u12_ao_transp (min) = ', (wall1 - wall0) / 60.d0
+  !call print_memory_usage()
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)]
+
+  implicit none
+  integer          :: ipoint
+  double precision :: wall0, wall1
+
+  PROVIDE mo_l_coef mo_r_coef
+  PROVIDE int2_grad1_u12_ao_transp
+
+  !print *, ' providing int2_grad1_u12_bimo_transp ...'
+  !call wall_time(wall0)
+
+  !$OMP PARALLEL         &
+  !$OMP DEFAULT (NONE)   &
+  !$OMP PRIVATE (ipoint) & 
+  !$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp)
+  !$OMP DO SCHEDULE (dynamic)
+  do ipoint = 1, n_points_final_grid
+    call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp  (1,1,1,ipoint), size(int2_grad1_u12_ao_transp  , 1) &
+                          , int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
+    call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp  (1,1,2,ipoint), size(int2_grad1_u12_ao_transp  , 1) &
+                          , int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
+    call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp  (1,1,3,ipoint), size(int2_grad1_u12_ao_transp  , 1) &
+                          , int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  !FREE int2_grad1_u12_ao_transp
+
+  !call wall_time(wall1)
+  !print *, ' wall time for int2_grad1_u12_bimo_transp (min) =', (wall1 - wall0) / 60.d0
+  !call print_memory_usage()
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_t_old, (n_points_final_grid, 3, mo_num, mo_num)]
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  !call wall_time(wall0)
+  !print *, ' providing int2_grad1_u12_bimo_t_old ...'
+
+  PROVIDE mo_l_coef mo_r_coef
+  PROVIDE int2_grad1_u12_bimo_transp
+
+  do ipoint = 1, n_points_final_grid
+    do i = 1, mo_num
+      do j = 1, mo_num
+        int2_grad1_u12_bimo_t_old(ipoint,1,j,i) = int2_grad1_u12_bimo_transp(j,i,1,ipoint)
+        int2_grad1_u12_bimo_t_old(ipoint,2,j,i) = int2_grad1_u12_bimo_transp(j,i,2,ipoint)
+        int2_grad1_u12_bimo_t_old(ipoint,3,j,i) = int2_grad1_u12_bimo_transp(j,i,3,ipoint)
+      enddo                                  
+    enddo
+  enddo
+
+  FREE int2_grad1_u12_bimo_transp
+
+  !call wall_time(wall1)
+  !print *, ' wall time for int2_grad1_u12_bimo_t_old (min) =', (wall1 - wall0) / 60.d0
+  !call print_memory_usage()
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)]
+
+  implicit none
+  integer          :: i, j, ipoint
+  double precision :: wall0, wall1
+
+  !call wall_time(wall0)
+  !print *, ' providing int2_grad1_u12_ao_t ...'
+
+  PROVIDE int2_grad1_u12_ao
+
+  do ipoint = 1, n_points_final_grid
+    do i = 1, ao_num
+      do j = 1, ao_num
+        int2_grad1_u12_ao_t(ipoint,1,j,i) = int2_grad1_u12_ao(j,i,ipoint,1)
+        int2_grad1_u12_ao_t(ipoint,2,j,i) = int2_grad1_u12_ao(j,i,ipoint,2)
+        int2_grad1_u12_ao_t(ipoint,3,j,i) = int2_grad1_u12_ao(j,i,ipoint,3)
+      enddo                                  
+    enddo
+  enddo
+
+  !call wall_time(wall1)
+  !print *, ' wall time for int2_grad1_u12_ao_t (min) =', (wall1 - wall0) / 60.d0
+  !call print_memory_usage()
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp, (n_points_final_grid, 3, mo_num, mo_num)]
+
+  implicit none
+  integer :: i, j, ipoint
+
+  do i = 1, mo_num
+    do j = 1, mo_num
+      do ipoint = 1, n_points_final_grid
+        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,1,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,1,ipoint)
+        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,2,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,2,ipoint)
+        mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,3,j,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu(j,i,3,ipoint)
+      enddo
+    enddo
+  enddo
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid, 3, mo_num, mo_num)]
+
+  BEGIN_DOC
+  !
+  ! x_W_ki_bi_ortho_erf_rk(ip,m,k,i) = \int dr chi_k(r) \frac{(1 - erf(mu |r-R_ip|))}{2|r-R_ip|} (x(m)-R_ip(m)) phi_i(r) ON THE BI-ORTHO MO BASIS 
+  !
+  ! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z,
+  !
+  ! R_ip = the "ip"-th point of the DFT Grid
+  END_DOC
+ 
+  implicit none
+  include 'constants.include.F'
+ 
+  integer          :: ipoint, m, i, k
+  double precision :: xyz
+  double precision :: wall0, wall1
+ 
+  !print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
+  !call wall_time(wall0)
+
+ !$OMP PARALLEL                   &
+ !$OMP DEFAULT (NONE)             &
+ !$OMP PRIVATE (ipoint,m,i,k,xyz) & 
+ !$OMP SHARED (x_W_ki_bi_ortho_erf_rk,n_points_final_grid,mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_num,final_grid_points) 
+ !$OMP DO SCHEDULE (dynamic)
+  do i = 1, mo_num
+    do k = 1, mo_num
+      do m = 1, 3
+        do ipoint = 1, n_points_final_grid
+          xyz = final_grid_points(m,ipoint)
+          x_W_ki_bi_ortho_erf_rk(ipoint,m,k,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,m,k,i) - xyz * mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,k,i)
+        enddo
+      enddo
+    enddo
+  enddo
+
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+ ! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu_transp 
+ ! FREE mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp
+
+  !call wall_time(wall1)
+  !print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_grid, 3, mo_num)]
+  BEGIN_DOC
+  ! x_W_ki_bi_ortho_erf_rk_diag(ip,m,i) = \int dr chi_i(r) (1 - erf(mu |r-R_ip|)) (x(m)-X(m)_ip) phi_i(r) ON THE BI-ORTHO MO BASIS 
+!
+! where chi_k(r)/phi_i(r) are left/right MOs, m=1 => X(m) = x, m=2 => X(m) = y, m=3 => X(m) = z,
+!
+! R_ip = the "ip"-th point of the DFT Grid
+  END_DOC
+
+  implicit none
+  include 'constants.include.F'
+ 
+  integer          :: ipoint, m, i
+  double precision :: xyz
+  double precision :: wall0, wall1
+ 
+  !print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...'
+  !call wall_time(wall0)
+
+ !$OMP PARALLEL                 &
+ !$OMP DEFAULT (NONE)           &
+ !$OMP PRIVATE (ipoint,m,i,xyz) & 
+ !$OMP SHARED (x_W_ki_bi_ortho_erf_rk_diag,n_points_final_grid,mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_v_ki_bi_ortho_erf_rk_cst_mu_transp,mo_num,final_grid_points) 
+ !$OMP DO SCHEDULE (dynamic)
+  do i = 1, mo_num
+    do m = 1, 3
+      do ipoint = 1, n_points_final_grid
+        xyz = final_grid_points(m,ipoint)
+        x_W_ki_bi_ortho_erf_rk_diag(ipoint,m,i) = mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,m,i,i) - xyz * mo_v_ki_bi_ortho_erf_rk_cst_mu_transp(ipoint,i,i)
+      enddo
+    enddo
+  enddo
+
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+  !call wall_time(wall1)
+  !print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0
+
+END_PROVIDER 
+
+! ---
+
diff --git a/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f b/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f
index 15ed2ce4..b1c2dc87 100644
--- a/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f
+++ b/plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f
@@ -16,7 +16,7 @@
   double precision              :: r(3)
   double precision, allocatable :: aos_r(:,:)
 
-  call cpu_time(tt0)
+  call wall_time(tt0)
 
   allocate(aos_r(ao_num,n_points_final_grid))
 
@@ -27,7 +27,7 @@
   call give_all_aos_at_r(r, aos_r(1,1))
 
 
-  call cpu_time(tt2)
+  call wall_time(tt2)
 
   !$OMP PARALLEL       &
   !$OMP DEFAULT (NONE) &
@@ -43,7 +43,7 @@
   !$OMP END DO
   !$OMP END PARALLEL
 
-  call cpu_time(tt3)
+  call wall_time(tt3)
   write(*,"(A,2X,F15.7)") ' wall time for AOs on r (sec) = ', (tt3 - tt2)
 
 
@@ -63,7 +63,7 @@
 
   deallocate(aos_r)
 
-  call cpu_time(tt1)
+  call wall_time(tt1)
   write(*,"(A,2X,F15.7)") ' wall time for mos_l_in_r_array_transp & mos_r_in_r_array_transp (sec) = ', (tt1 - tt0)
 
 END_PROVIDER