diff --git a/src/bi_ort_ints/total_twoe_pot.irp.f b/src/bi_ort_ints/total_twoe_pot.irp.f
index f5f5959a..bdebe890 100644
--- a/src/bi_ort_ints/total_twoe_pot.irp.f
+++ b/src/bi_ort_ints/total_twoe_pot.irp.f
@@ -40,20 +40,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
 
   provide j1b_type
 
-  if(j1b_type .eq. 3) then
-
-    do j = 1, ao_num
-      do l = 1, ao_num
-        do i = 1, ao_num
-          do k = 1, ao_num
-            ao_two_e_tc_tot(k,i,l,j) = ao_tc_int_chemist(k,i,l,j)
-            !write(222,*) ao_two_e_tc_tot(k,i,l,j) 
-          enddo
-        enddo
-      enddo
-    enddo
-
-  else
+  if(j1b_type .eq. 0) then
 
     PROVIDE ao_tc_sym_two_e_pot_in_map
 
@@ -77,6 +64,19 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
       enddo
     enddo
 
+  else
+
+    do j = 1, ao_num
+      do l = 1, ao_num
+        do i = 1, ao_num
+          do k = 1, ao_num
+            ao_two_e_tc_tot(k,i,l,j) = ao_tc_int_chemist(k,i,l,j)
+            !write(222,*) ao_two_e_tc_tot(k,i,l,j) 
+          enddo
+        enddo
+      enddo
+    enddo
+
   endif
 
 END_PROVIDER 
diff --git a/src/non_h_ints_mu/debug_fit.irp.f b/src/non_h_ints_mu/debug_fit.irp.f
index af441335..5995bffa 100644
--- a/src/non_h_ints_mu/debug_fit.irp.f
+++ b/src/non_h_ints_mu/debug_fit.irp.f
@@ -82,9 +82,9 @@ subroutine test_grad_j1b_nucl()
   integer                    :: ipoint
   double precision           :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz
   double precision           :: r(3)
-  double precision, external :: grad_x_j1b_nucl
-  double precision, external :: grad_y_j1b_nucl
-  double precision, external :: grad_z_j1b_nucl
+  double precision, external :: grad_x_j1b_nucl_num
+  double precision, external :: grad_y_j1b_nucl_num
+  double precision, external :: grad_z_j1b_nucl_num
 
   print*, ' test_grad_j1b_nucl ...'
 
@@ -101,7 +101,7 @@ subroutine test_grad_j1b_nucl()
     r(3) = final_grid_points(3,ipoint)
 
     i_exc  = v_1b_grad(1,ipoint) 
-    i_num  = grad_x_j1b_nucl(r)
+    i_num  = grad_x_j1b_nucl_num(r)
     acc_ij = dabs(i_exc - i_num)
     if(acc_ij .gt. eps_ij) then
       print *, ' problem in x of v_1b_grad on', ipoint
@@ -111,7 +111,7 @@ subroutine test_grad_j1b_nucl()
     endif
 
     i_exc  = v_1b_grad(2,ipoint) 
-    i_num  = grad_y_j1b_nucl(r)
+    i_num  = grad_y_j1b_nucl_num(r)
     acc_ij = dabs(i_exc - i_num)
     if(acc_ij .gt. eps_ij) then
       print *, ' problem in y of v_1b_grad on', ipoint
@@ -121,7 +121,7 @@ subroutine test_grad_j1b_nucl()
     endif
 
     i_exc  = v_1b_grad(3,ipoint) 
-    i_num  = grad_z_j1b_nucl(r)
+    i_num  = grad_z_j1b_nucl_num(r)
     acc_ij = dabs(i_exc - i_num)
     if(acc_ij .gt. eps_ij) then
       print *, ' problem in z of v_1b_grad on', ipoint
@@ -317,7 +317,7 @@ subroutine test_fit_ugradu()
       i_fit = i_fit / dsqrt(x2)
   
       tmp = j12_mu(r1, r2) 
-      call grad1_j12_mu_exc(r1, r2, grad)
+      call grad1_j12_mu(r1, r2, grad)
   
       ! ---
   
diff --git a/src/non_h_ints_mu/grad_squared.irp.f b/src/non_h_ints_mu/grad_squared.irp.f
index 1fd39f6a..8b801f9d 100644
--- a/src/non_h_ints_mu/grad_squared.irp.f
+++ b/src/non_h_ints_mu/grad_squared.irp.f
@@ -232,37 +232,33 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
 
   PROVIDE j1b_type
 
-  if(j1b_type .eq. 3) then
-
-    do ipoint = 1, n_points_final_grid
-      tmp1 = v_1b(ipoint)
-      tmp1 = tmp1 * tmp1
-      do j = 1, ao_num
-        do i = 1, ao_num
-          grad12_j12(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
-        enddo
+  do ipoint = 1, n_points_final_grid
+    tmp1 = v_1b(ipoint)
+    tmp1 = tmp1 * tmp1
+    do j = 1, ao_num
+      do i = 1, ao_num
+        grad12_j12(i,j,ipoint) = tmp1 * int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
       enddo
     enddo
+  enddo
 
-  else
-
-    grad12_j12 = 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)
-      do j = 1, ao_num
-        do i = 1, ao_num
-          do igauss = 1, n_max_fit_slat
-            delta = expo_gauss_1_erf_x_2(igauss)
-            coef  = coef_gauss_1_erf_x_2(igauss)
-            grad12_j12(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j)
-          enddo
-        enddo
-      enddo
-    enddo
-
-  endif
+  !if(j1b_type .eq. 0) then 
+  !  grad12_j12 = 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)
+  !    do j = 1, ao_num
+  !      do i = 1, ao_num
+  !        do igauss = 1, n_max_fit_slat
+  !          delta = expo_gauss_1_erf_x_2(igauss)
+  !          coef  = coef_gauss_1_erf_x_2(igauss)
+  !          grad12_j12(i,j,ipoint) += -0.25d0 * coef * overlap_gauss_r12_ao(r, delta, i, j)
+  !        enddo
+  !      enddo
+  !    enddo
+  !  enddo
+  !endif
 
   call wall_time(time1)
   print*, ' Wall time for grad12_j12 = ', time1 - time0
@@ -351,6 +347,7 @@ END_PROVIDER
 
 ! ---
 
+
 BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao_num)]
 
   BEGIN_DOC
@@ -376,7 +373,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   else
 
-    allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
+    allocate(b_mat(n_points_final_grid,ao_num,ao_num))
 
     b_mat = 0.d0
    !$OMP PARALLEL               &
@@ -392,29 +389,13 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
       enddo
     enddo
    !$OMP END DO
-   !$OMP END PARALLEL
-
-    tmp = 0.d0
-   !$OMP PARALLEL               &
-   !$OMP DEFAULT (NONE)         &
-   !$OMP PRIVATE (j, l, ipoint) &
-   !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
-   !$OMP DO SCHEDULE (static)
-    do ipoint = 1, n_points_final_grid
-      do j = 1, ao_num
-        do l = 1, ao_num
-          tmp(l,j,ipoint) = u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) + 0.5d0 * grad12_j12(l,j,ipoint)
-        enddo
-      enddo
-    enddo
-   !$OMP END DO
    !$OMP END PARALLEL
 
     tc_grad_square_ao = 0.d0
-    call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
-              , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid      &
-              , 1.d0, tc_grad_square_ao, ao_num*ao_num)
-    deallocate(tmp, b_mat)
+    call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0                  &
+              , int2_grad1_u12_squared_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
+              , 0.d0, tc_grad_square_ao, ao_num*ao_num)
+    deallocate(b_mat)
 
     call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)
 
@@ -450,3 +431,4 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 END_PROVIDER
 
 ! ---
+
diff --git a/src/non_h_ints_mu/grad_tc_int.irp.f b/src/non_h_ints_mu/grad_tc_int.irp.f
index cb3b71a3..f4eb02e2 100644
--- a/src/non_h_ints_mu/grad_tc_int.irp.f
+++ b/src/non_h_ints_mu/grad_tc_int.irp.f
@@ -16,9 +16,11 @@ BEGIN_PROVIDER [double precision, ao_non_hermit_term_chemist, (ao_num, ao_num, a
   double precision              :: wall1, wall0
   double precision, allocatable :: b_mat(:,:,:,:), ac_mat(:,:,:,:)
 
+  print*, ' providing ao_non_hermit_term_chemist ...'
+  call wall_time(wall0)
+
   provide v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
 
-  call wall_time(wall0)
   allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
 
  !$OMP PARALLEL                         &
@@ -102,7 +104,7 @@ BEGIN_PROVIDER [double precision, ao_non_hermit_term_chemist, (ao_num, ao_num, a
  !$OMP END PARALLEL
 
   call wall_time(wall1)
-  print *, ' wall time dgemm ', wall1 - wall0
+  print *, ' wall time for ao_non_hermit_term_chemist ', wall1 - wall0
 
 END_PROVIDER 
 
diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f
index a515e0b8..c9f62b18 100644
--- a/src/non_h_ints_mu/j12_nucl_utils.irp.f
+++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f
@@ -204,39 +204,6 @@ END_PROVIDER
 
 ! ---
 
-double precision function jmu_modif(r1, r2)
-
-  implicit none
-  double precision, intent(in) :: r1(3), r2(3)
-  double precision, external   :: j12_mu, j12_nucl
-
-  jmu_modif = j12_mu(r1, r2) * j12_nucl(r1, r2)
-
-  return
-end function jmu_modif
-
-! ---
-
-double precision function j12_mu(r1, r2)
-
-  include 'constants.include.F'
-
-  implicit none
-  double precision, intent(in) :: r1(3), r2(3)
-  double precision             :: mu_r12, r12
-
-  r12 = 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)) )
-  mu_r12 = mu_erf * r12
-
-  j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf
-
-  return
-end function j12_mu
-
-! ---
-
 double precision function j12_mu_r12(r12)
 
   include 'constants.include.F'
@@ -254,6 +221,19 @@ end function j12_mu_r12
 
 ! ---
 
+double precision function jmu_modif(r1, r2)
+
+  implicit none
+  double precision, intent(in) :: r1(3), r2(3)
+  double precision, external   :: j12_mu, j12_nucl
+
+  jmu_modif = j12_mu(r1, r2) * j12_nucl(r1, r2)
+
+  return
+end function jmu_modif
+
+! ---
+
 double precision function j12_mu_gauss(r1, r2)
 
   implicit none
@@ -278,30 +258,6 @@ end function j12_mu_gauss
 
 ! ---
 
-double precision function j1b_nucl(r)
-
-  implicit none
-  double precision, intent(in) :: r(3)
-  integer                      :: i
-  double precision             :: a, d, e
-
-  j1b_nucl = 1.d0
-
-  do i = 1, nucl_num
-    a = j1b_pen(i)
-    d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
-        + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
-        + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
-    e = 1.d0 - exp(-a*d)
-
-    j1b_nucl = j1b_nucl * e
-  enddo
-
-  return
-end function j1b_nucl
-
-! ---
-
 double precision function j12_nucl(r1, r2)
 
   implicit none
@@ -317,7 +273,7 @@ end function j12_nucl
 
 ! ---------------------------------------------------------------------------------------
 
-double precision function grad_x_j1b_nucl(r)
+double precision function grad_x_j1b_nucl_num(r)
 
   implicit none
   double precision, intent(in) :: r(3)
@@ -333,12 +289,12 @@ double precision function grad_x_j1b_nucl(r)
   r_eps(1) = r_eps(1) - 2.d0 * delta
   fm       = j1b_nucl(r_eps)
 
-  grad_x_j1b_nucl = 0.5d0 * (fp - fm) / delta
+  grad_x_j1b_nucl_num = 0.5d0 * (fp - fm) / delta
 
   return
-end function grad_x_j1b_nucl
+end function grad_x_j1b_nucl_num
 
-double precision function grad_y_j1b_nucl(r)
+double precision function grad_y_j1b_nucl_num(r)
 
   implicit none
   double precision, intent(in) :: r(3)
@@ -354,12 +310,12 @@ double precision function grad_y_j1b_nucl(r)
   r_eps(2) = r_eps(2) - 2.d0 * delta
   fm       = j1b_nucl(r_eps)
 
-  grad_y_j1b_nucl = 0.5d0 * (fp - fm) / delta
+  grad_y_j1b_nucl_num = 0.5d0 * (fp - fm) / delta
 
   return
-end function grad_y_j1b_nucl
+end function grad_y_j1b_nucl_num
 
-double precision function grad_z_j1b_nucl(r)
+double precision function grad_z_j1b_nucl_num(r)
 
   implicit none
   double precision, intent(in) :: r(3)
@@ -375,10 +331,10 @@ double precision function grad_z_j1b_nucl(r)
   r_eps(3) = r_eps(3) - 2.d0 * delta
   fm       = j1b_nucl(r_eps)
 
-  grad_z_j1b_nucl = 0.5d0 * (fp - fm) / delta
+  grad_z_j1b_nucl_num = 0.5d0 * (fp - fm) / delta
 
   return
-end function grad_z_j1b_nucl
+end function grad_z_j1b_nucl_num
 
 ! ---------------------------------------------------------------------------------------
 
@@ -389,9 +345,9 @@ double precision function lapl_j1b_nucl(r)
   implicit none
   double precision, intent(in) :: r(3)
   double precision             :: r_eps(3), eps, fp, fm, delta
-  double precision, external   :: grad_x_j1b_nucl
-  double precision, external   :: grad_y_j1b_nucl
-  double precision, external   :: grad_z_j1b_nucl
+  double precision, external   :: grad_x_j1b_nucl_num
+  double precision, external   :: grad_y_j1b_nucl_num
+  double precision, external   :: grad_z_j1b_nucl_num
 
   eps   = 1d-5
   r_eps = r
@@ -402,9 +358,9 @@ double precision function lapl_j1b_nucl(r)
 
   delta    = max(eps, dabs(eps*r(1)))
   r_eps(1) = r_eps(1) + delta
-  fp       = grad_x_j1b_nucl(r_eps)
+  fp       = grad_x_j1b_nucl_num(r_eps)
   r_eps(1) = r_eps(1) - 2.d0 * delta
-  fm       = grad_x_j1b_nucl(r_eps)
+  fm       = grad_x_j1b_nucl_num(r_eps)
   r_eps(1) = r_eps(1) + delta
 
   lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta
@@ -413,9 +369,9 @@ double precision function lapl_j1b_nucl(r)
 
   delta    = max(eps, dabs(eps*r(2)))
   r_eps(2) = r_eps(2) + delta
-  fp       = grad_y_j1b_nucl(r_eps)
+  fp       = grad_y_j1b_nucl_num(r_eps)
   r_eps(2) = r_eps(2) - 2.d0 * delta
-  fm       = grad_y_j1b_nucl(r_eps)
+  fm       = grad_y_j1b_nucl_num(r_eps)
   r_eps(2) = r_eps(2) + delta
 
   lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta
@@ -424,9 +380,9 @@ double precision function lapl_j1b_nucl(r)
 
   delta    = max(eps, dabs(eps*r(3)))
   r_eps(3) = r_eps(3) + delta
-  fp       = grad_z_j1b_nucl(r_eps)
+  fp       = grad_z_j1b_nucl_num(r_eps)
   r_eps(3) = r_eps(3) - 2.d0 * delta
-  fm       = grad_z_j1b_nucl(r_eps)
+  fm       = grad_z_j1b_nucl_num(r_eps)
   r_eps(3) = r_eps(3) + delta
 
   lapl_j1b_nucl += 0.5d0 * (fp - fm) / delta
@@ -574,35 +530,6 @@ end function grad1_z_j12_mu_num
 
 ! ---------------------------------------------------------------------------------------
 
-! ---
-
-subroutine grad1_j12_mu_exc(r1, r2, grad)
-
-  implicit none
-  double precision, intent(in)  :: r1(3), r2(3)
-  double precision, intent(out) :: grad(3)
-  double precision              :: dx, dy, dz, r12, tmp
-
-  grad = 0.d0
-
-  dx = r1(1) - r2(1)
-  dy = r1(2) - r2(2)
-  dz = r1(3) - r2(3)
-
-  r12 = dsqrt( dx * dx + dy * dy + dz * dz )
-  if(r12 .lt. 1d-10) return
-
-  tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
-
-  grad(1) = tmp * dx 
-  grad(2) = tmp * dy 
-  grad(3) = tmp * dz 
-
-  return
-end subroutine grad1_j12_mu_exc
-
-! ---
-
 subroutine grad1_jmu_modif_num(r1, r2, grad)
 
   implicit none
@@ -614,11 +541,11 @@ subroutine grad1_jmu_modif_num(r1, r2, grad)
 
   double precision, external    :: j12_mu
   double precision, external    :: j1b_nucl
-  double precision, external    :: grad_x_j1b_nucl
-  double precision, external    :: grad_y_j1b_nucl
-  double precision, external    :: grad_z_j1b_nucl
+  double precision, external    :: grad_x_j1b_nucl_num
+  double precision, external    :: grad_y_j1b_nucl_num
+  double precision, external    :: grad_z_j1b_nucl_num
 
-  call grad1_j12_mu_exc(r1, r2, grad_u12)
+  call grad1_j12_mu(r1, r2, grad_u12)
 
   tmp0 = j1b_nucl(r1) 
   tmp1 = j1b_nucl(r2)
@@ -626,9 +553,9 @@ subroutine grad1_jmu_modif_num(r1, r2, grad)
   tmp3 = tmp0 * tmp1
   tmp4 = tmp2 * tmp1
 
-  grad(1) = tmp3 * grad_u12(1) + tmp4 * grad_x_j1b_nucl(r1)
-  grad(2) = tmp3 * grad_u12(2) + tmp4 * grad_y_j1b_nucl(r1)
-  grad(3) = tmp3 * grad_u12(3) + tmp4 * grad_z_j1b_nucl(r1)
+  grad(1) = tmp3 * grad_u12(1) + tmp4 * grad_x_j1b_nucl_num(r1)
+  grad(2) = tmp3 * grad_u12(2) + tmp4 * grad_y_j1b_nucl_num(r1)
+  grad(3) = tmp3 * grad_u12(3) + tmp4 * grad_z_j1b_nucl_num(r1)
 
   return
 end subroutine grad1_jmu_modif_num
diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f
new file mode 100644
index 00000000..b2772f92
--- /dev/null
+++ b/src/non_h_ints_mu/jast_deriv.irp.f
@@ -0,0 +1,245 @@
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, grad1_u12_num,         (n_points_extra_final_grid, n_points_final_grid, 3)]
+&BEGIN_PROVIDER [ double precision, grad1_u12_squared_num, (n_points_extra_final_grid, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! grad_1 u(r1,r2)
+  !
+  ! this will be integrated numerically over r2:
+  !   we use grid for r1 and extra_grid for r2
+  !
+  ! for 99 < j1b_type < 199
+  !
+  !   u(r1,r2)        = j12_mu(r12) x v(r1) x v(r2)
+  !   grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2)
+  !
+  END_DOC
+
+  implicit none
+  integer          :: ipoint, jpoint
+  double precision :: r1(3), r2(3)
+
+  PROVIDE j1b_type
+  PROVIDE final_grid_points_extra
+
+  grad1_u12_num         = 0.d0
+  grad1_u12_squared_num = 0.d0
+
+  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+
+    double precision           :: v1b_r1, v1b_r2, u2b_r12
+    double precision           :: grad1_v1b(3), grad1_u2b(3)
+    double precision           :: dx, dy, dz
+    double precision, external :: j12_mu, j1b_nucl
+
+    !$OMP PARALLEL                                                                                    &
+    !$OMP DEFAULT (NONE)                                                                              &
+    !$OMP PRIVATE (ipoint, jpoint, r1, r2, v1b_r1, v1b_r2, u2b_r12, grad1_v1b, grad1_u2b, dx, dy, dz) &
+    !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points,                  &
+    !$OMP         final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
+    !$OMP DO SCHEDULE (static)
+    do ipoint = 1, n_points_final_grid  ! r1
+
+      r1(1) = final_grid_points(1,ipoint)
+      r1(2) = final_grid_points(2,ipoint)
+      r1(3) = final_grid_points(3,ipoint)
+
+      v1b_r1 = j1b_nucl(r1)
+      call grad1_j1b_nucl(r1, grad1_v1b)
+
+      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)
+
+        v1b_r2  = j1b_nucl(r2)
+        u2b_r12 = j12_mu(r1, r2)
+        call grad1_j12_mu(r1, r2, grad1_u2b)
+
+        dx = (grad1_u2b(1) * v1b_r1 + u2b_r12 * grad1_v1b(1)) * v1b_r2
+        dy = (grad1_u2b(2) * v1b_r1 + u2b_r12 * grad1_v1b(2)) * v1b_r2
+        dz = (grad1_u2b(3) * v1b_r1 + u2b_r12 * grad1_v1b(3)) * v1b_r2
+
+        grad1_u12_num(jpoint,ipoint,1) = dx 
+        grad1_u12_num(jpoint,ipoint,2) = dy 
+        grad1_u12_num(jpoint,ipoint,3) = dz 
+
+        grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+
+  endif
+
+END_PROVIDER
+
+! ---
+
+double precision function j12_mu(r1, r2)
+
+  include 'constants.include.F'
+
+  implicit none
+  double precision, intent(in) :: r1(3), r2(3)
+  double precision             :: mu_r12, r12
+
+  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+
+    r12 = 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)) )
+    mu_r12 = mu_erf * r12
+
+    j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+
+  endif
+
+  return
+end function j12_mu
+
+! ---
+
+subroutine grad1_j12_mu(r1, r2, grad)
+
+  implicit none
+  double precision, intent(in)  :: r1(3), r2(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: dx, dy, dz, r12, tmp
+
+  grad = 0.d0
+
+  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+
+    dx = r1(1) - r2(1)
+    dy = r1(2) - r2(2)
+    dz = r1(3) - r2(3)
+
+    r12 = dsqrt( dx * dx + dy * dy + dz * dz )
+    if(r12 .lt. 1d-10) return
+
+    tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
+
+    grad(1) = tmp * dx
+    grad(2) = tmp * dy
+    grad(3) = tmp * dz
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+
+  endif
+
+  return
+end subroutine grad1_j12_mu
+
+! ---
+
+double precision function j1b_nucl(r)
+
+  implicit none
+  double precision, intent(in) :: r(3)
+  integer                      :: i
+  double precision             :: a, d, e
+
+  if(j1b_type .eq. 103) then
+
+    j1b_nucl = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      e = 1.d0 - dexp(-a*d)
+      j1b_nucl = j1b_nucl * e
+    enddo
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+
+  endif
+
+  return
+end function j1b_nucl
+
+! ---
+
+subroutine grad1_j1b_nucl(r, grad)
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: grad(3)
+  integer                       :: ipoint, i, j, phase
+  double precision              :: x, y, z, dx, dy, dz
+  double precision              :: a, d, e
+  double precision              :: fact_x, fact_y, fact_z
+  double precision              :: ax_der, ay_der, az_der, a_expo
+
+  if(j1b_type .eq. 103) then
+
+    x = r(1)
+    y = r(2)
+    z = r(3)
+
+    fact_x = 0.d0
+    fact_y = 0.d0
+    fact_z = 0.d0
+    do i = 1, List_all_comb_b2_size
+
+      phase  = 0
+      a_expo = 0.d0
+      ax_der = 0.d0
+      ay_der = 0.d0
+      az_der = 0.d0
+      do j = 1, nucl_num
+        a  = dble(List_all_comb_b2(j,i)) * j1b_pen(j)
+        dx = x - nucl_coord(j,1)
+        dy = y - nucl_coord(j,2)
+        dz = z - nucl_coord(j,3)
+
+        phase  += List_all_comb_b2(j,i)
+        a_expo += a * (dx*dx + dy*dy + dz*dz)
+        ax_der += a * dx
+        ay_der += a * dy
+        az_der += a * dz
+      enddo
+      e = -2.d0 * (-1.d0)**dble(phase) * dexp(-a_expo)
+
+      fact_x += e * ax_der
+      fact_y += e * ay_der
+      fact_z += e * az_der
+    enddo
+
+    grad(1) = fact_x
+    grad(2) = fact_y
+    grad(3) = fact_z
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+
+  endif
+
+  return
+end subroutine grad1_j1b_nucl
+
+! ---
+
diff --git a/src/non_h_ints_mu/new_grad_tc.irp.f b/src/non_h_ints_mu/new_grad_tc.irp.f
index a6673252..064eb9f1 100644
--- a/src/non_h_ints_mu/new_grad_tc.irp.f
+++ b/src/non_h_ints_mu/new_grad_tc.irp.f
@@ -1,101 +1,5 @@
 ! ---
 
-BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)]
-
-  BEGIN_DOC
-  !
-  ! int2_grad1_u12_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) 
-  !
-  ! where r1 = r(ipoint)
-  !
-  ! if J(r1,r2) = u12:
-  !
-  ! int2_grad1_u12_ao(i,j,ipoint,:) = 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r2) \phi_j(r2)
-  !                                 = 0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ]
-  !
-  ! if J(r1,r2) = u12 x v1 x v2
-  !
-  ! int2_grad1_u12_ao(i,j,ipoint,:) =      v1    x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
-  !                                 - \grad_1 v1 x [       \int dr2                  u12              v2 \phi_i(r2) \phi_j(r2) ] 
-  !                                 =    0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) 
-  !                                 -    0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) 
-  !                                 - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
-  !
-  !
-  END_DOC
-
-  implicit none
-  integer          :: ipoint, i, j, m
-  double precision :: time0, time1
-  double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
-
-  print*, ' providing int2_grad1_u12_ao ...'
-  call wall_time(time0)
-
-  PROVIDE j1b_type
-
-  if(read_tc_integ) then
-
-    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read")
-    read(11) int2_grad1_u12_ao
-    close(11)
-
-  else
-
-    if(j1b_type .eq. 3) then
-      do ipoint = 1, n_points_final_grid
-        x = final_grid_points(1,ipoint)
-        y = final_grid_points(2,ipoint)
-        z = final_grid_points(3,ipoint)
-        tmp0  = 0.5d0 * v_1b(ipoint)
-        tmp_x =  v_1b_grad(1,ipoint)
-        tmp_y =  v_1b_grad(2,ipoint)
-        tmp_z =  v_1b_grad(3,ipoint)
-        do j = 1, ao_num
-          do i = 1, ao_num
-            tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
-            tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
-            int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
-            int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
-            int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
-          enddo
-        enddo
-      enddo
-    else
-      do ipoint = 1, n_points_final_grid
-        x = final_grid_points(1,ipoint)
-        y = final_grid_points(2,ipoint)
-        z = final_grid_points(3,ipoint)
-        do j = 1, ao_num
-          do i = 1, ao_num
-            tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
-
-            int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,1)
-            int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,2)
-            int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,3)
-          enddo
-        enddo
-      enddo
-      int2_grad1_u12_ao *= 0.5d0
-    endif
-
-  endif
-
-  if(write_tc_integ.and.mpi_master) then
-    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write")
-    call ezfio_set_work_empty(.False.)
-    write(11) int2_grad1_u12_ao
-    close(11)
-    call ezfio_set_tc_keywords_io_tc_integ('Read')
-  endif
-
-  call wall_time(time1)
-  print*, ' Wall time for int2_grad1_u12_ao = ', time1 - time0
-
-END_PROVIDER 
-
-! ---
-
 BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)]
 
   BEGIN_DOC
@@ -303,6 +207,8 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
   print*, ' providing tc_grad_and_lapl_ao ...'
   call wall_time(time0)
 
+  PROVIDE int2_grad1_u12_ao
+
   if(read_tc_integ) then
 
     open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao', action="read")
@@ -341,8 +247,7 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
     do m = 1, 3
       call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0              &
                 , int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
-                , 1.d0, tc_grad_and_lapl_ao, ao_num*ao_num) 
-  
+                , 0.d0, tc_grad_and_lapl_ao, ao_num*ao_num) 
     enddo
     deallocate(b_mat)
   
diff --git a/src/non_h_ints_mu/new_grad_tc_manu.irp.f b/src/non_h_ints_mu/new_grad_tc_manu.irp.f
index 901e3048..7ab5b327 100644
--- a/src/non_h_ints_mu/new_grad_tc_manu.irp.f
+++ b/src/non_h_ints_mu/new_grad_tc_manu.irp.f
@@ -39,7 +39,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
     read(11) int2_grad1_u12_ao_test
     close(11)
 
-
   else
 
     if(j1b_type .eq. 3) then
diff --git a/src/non_h_ints_mu/numerical_integ.irp.f b/src/non_h_ints_mu/numerical_integ.irp.f
index dcd7a52a..f9457247 100644
--- a/src/non_h_ints_mu/numerical_integ.irp.f
+++ b/src/non_h_ints_mu/numerical_integ.irp.f
@@ -322,9 +322,9 @@ double precision function num_gradu_squared_u_ij_mu(i, j, ipoint)
   double precision, external :: ao_value
   double precision, external :: j1b_nucl
   double precision, external :: j12_mu
-  double precision, external :: grad_x_j1b_nucl
-  double precision, external :: grad_y_j1b_nucl
-  double precision, external :: grad_z_j1b_nucl
+  double precision, external :: grad_x_j1b_nucl_num
+  double precision, external :: grad_y_j1b_nucl_num
+  double precision, external :: grad_z_j1b_nucl_num
 
   r1(1) = final_grid_points(1,ipoint)
   r1(2) = final_grid_points(2,ipoint)
@@ -342,11 +342,11 @@ double precision function num_gradu_squared_u_ij_mu(i, j, ipoint)
     tmp_z = r1(3) - r2(3)
     r12   = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
 
-    dx1_v1 = grad_x_j1b_nucl(r1)
-    dy1_v1 = grad_y_j1b_nucl(r1)
-    dz1_v1 = grad_z_j1b_nucl(r1)
+    dx1_v1 = grad_x_j1b_nucl_num(r1)
+    dy1_v1 = grad_y_j1b_nucl_num(r1)
+    dz1_v1 = grad_z_j1b_nucl_num(r1)
 
-    call grad1_j12_mu_exc(r1, r2, grad_u12)
+    call grad1_j12_mu(r1, r2, grad_u12)
 
     tmp1    = 1.d0 - derf(mu_erf * r12)
     v1_tmp  = j1b_nucl(r1)
@@ -390,9 +390,9 @@ double precision function num_grad12_j12(i, j, ipoint)
   double precision, external :: ao_value
   double precision, external :: j1b_nucl
   double precision, external :: j12_mu
-  double precision, external :: grad_x_j1b_nucl
-  double precision, external :: grad_y_j1b_nucl
-  double precision, external :: grad_z_j1b_nucl
+  double precision, external :: grad_x_j1b_nucl_num
+  double precision, external :: grad_y_j1b_nucl_num
+  double precision, external :: grad_z_j1b_nucl_num
 
   r1(1) = final_grid_points(1,ipoint)
   r1(2) = final_grid_points(2,ipoint)
@@ -410,11 +410,11 @@ double precision function num_grad12_j12(i, j, ipoint)
     tmp_z = r1(3) - r2(3)
     r12   = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
 
-    dx1_v1 = grad_x_j1b_nucl(r1)
-    dy1_v1 = grad_y_j1b_nucl(r1)
-    dz1_v1 = grad_z_j1b_nucl(r1)
+    dx1_v1 = grad_x_j1b_nucl_num(r1)
+    dy1_v1 = grad_y_j1b_nucl_num(r1)
+    dz1_v1 = grad_z_j1b_nucl_num(r1)
 
-    call grad1_j12_mu_exc(r1, r2, grad_u12)
+    call grad1_j12_mu(r1, r2, grad_u12)
 
     tmp1    = 1.d0 - derf(mu_erf * r12)
     v1_tmp  = j1b_nucl(r1)
@@ -456,9 +456,9 @@ double precision function num_u12sq_j1bsq(i, j, ipoint)
   double precision, external :: ao_value
   double precision, external :: j1b_nucl
   double precision, external :: j12_mu
-  double precision, external :: grad_x_j1b_nucl
-  double precision, external :: grad_y_j1b_nucl
-  double precision, external :: grad_z_j1b_nucl
+  double precision, external :: grad_x_j1b_nucl_num
+  double precision, external :: grad_y_j1b_nucl_num
+  double precision, external :: grad_z_j1b_nucl_num
 
   r1(1) = final_grid_points(1,ipoint)
   r1(2) = final_grid_points(2,ipoint)
@@ -476,11 +476,11 @@ double precision function num_u12sq_j1bsq(i, j, ipoint)
     tmp_z = r1(3) - r2(3)
     r12   = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
 
-    dx1_v1 = grad_x_j1b_nucl(r1)
-    dy1_v1 = grad_y_j1b_nucl(r1)
-    dz1_v1 = grad_z_j1b_nucl(r1)
+    dx1_v1 = grad_x_j1b_nucl_num(r1)
+    dy1_v1 = grad_y_j1b_nucl_num(r1)
+    dz1_v1 = grad_z_j1b_nucl_num(r1)
 
-    call grad1_j12_mu_exc(r1, r2, grad_u12)
+    call grad1_j12_mu(r1, r2, grad_u12)
 
     tmp1    = 1.d0 - derf(mu_erf * r12)
     v1_tmp  = j1b_nucl(r1)
@@ -522,9 +522,9 @@ double precision function num_u12_grad1_u12_j1b_grad1_j1b(i, j, ipoint)
   double precision, external :: ao_value
   double precision, external :: j1b_nucl
   double precision, external :: j12_mu
-  double precision, external :: grad_x_j1b_nucl
-  double precision, external :: grad_y_j1b_nucl
-  double precision, external :: grad_z_j1b_nucl
+  double precision, external :: grad_x_j1b_nucl_num
+  double precision, external :: grad_y_j1b_nucl_num
+  double precision, external :: grad_z_j1b_nucl_num
 
   r1(1) = final_grid_points(1,ipoint)
   r1(2) = final_grid_points(2,ipoint)
@@ -542,11 +542,11 @@ double precision function num_u12_grad1_u12_j1b_grad1_j1b(i, j, ipoint)
     tmp_z = r1(3) - r2(3)
     r12   = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
 
-    dx1_v1 = grad_x_j1b_nucl(r1)
-    dy1_v1 = grad_y_j1b_nucl(r1)
-    dz1_v1 = grad_z_j1b_nucl(r1)
+    dx1_v1 = grad_x_j1b_nucl_num(r1)
+    dy1_v1 = grad_y_j1b_nucl_num(r1)
+    dz1_v1 = grad_z_j1b_nucl_num(r1)
 
-    call grad1_j12_mu_exc(r1, r2, grad_u12)
+    call grad1_j12_mu(r1, r2, grad_u12)
 
     tmp1    = 1.d0 - derf(mu_erf * r12)
     v1_tmp  = j1b_nucl(r1)
diff --git a/src/non_h_ints_mu/tc_integ.irp.f b/src/non_h_ints_mu/tc_integ.irp.f
new file mode 100644
index 00000000..865522e0
--- /dev/null
+++ b/src/non_h_ints_mu/tc_integ.irp.f
@@ -0,0 +1,256 @@
+
+! ---
+
+ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao,         (ao_num, ao_num, n_points_final_grid, 3)]
+&BEGIN_PROVIDER [double precision, int2_grad1_u12_squared_ao, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! int2_grad1_u12_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) 
+  !
+  ! where r1 = r(ipoint)
+  !
+  ! if J(r1,r2) = u12 (j1b_type .eq. 1)
+  !
+  ! int2_grad1_u12_ao(i,j,ipoint,:) = 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r2) \phi_j(r2)
+  !                                 = 0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ]
+  !
+  ! if J(r1,r2) = u12 x v1 x v2 (j1b_type .eq. 3)
+  !
+  ! int2_grad1_u12_ao(i,j,ipoint,:) =      v1    x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
+  !                                 - \grad_1 v1 x [       \int dr2                  u12              v2 \phi_i(r2) \phi_j(r2) ] 
+  !                                 =    0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) 
+  !                                 -    0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) 
+  !                                 - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
+  !
+  !
+  !
+  ! int2_grad1_u12_squared_ao = int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2
+  !
+  END_DOC
+
+  implicit none
+  integer          :: ipoint, i, j, m, jpoint
+  double precision :: time0, time1
+  double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
+
+  print*, ' providing int2_grad1_u12_ao & int2_grad1_u12_squared_ao ...'
+
+  PROVIDE j1b_type
+
+  if(read_tc_integ) then
+    call wall_time(time0)
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read")
+    read(11) int2_grad1_u12_ao
+    call wall_time(time1)
+    print*, ' wall time for int2_grad1_u12_ao =', time1-time0 
+  endif
+
+  if(j1b_type .eq. 3) then
+
+    ! ---
+
+    if(.not.read_tc_integ) then
+      call wall_time(time0)
+      PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
+      int2_grad1_u12_ao = 0.d0
+      ! TODO OPENMP
+      do ipoint = 1, n_points_final_grid
+        x = final_grid_points(1,ipoint)
+        y = final_grid_points(2,ipoint)
+        z = final_grid_points(3,ipoint)
+        tmp0  = 0.5d0 * v_1b(ipoint)
+        tmp_x =  v_1b_grad(1,ipoint)
+        tmp_y =  v_1b_grad(2,ipoint)
+        tmp_z =  v_1b_grad(3,ipoint)
+        do j = 1, ao_num
+          do i = 1, ao_num
+            tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
+            tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
+            int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
+            int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
+            int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
+          enddo
+        enddo
+      enddo
+      call wall_time(time1)
+      print*, ' wall time for int2_grad1_u12_ao =', time1-time0 
+    endif
+
+    ! ---
+
+    call wall_time(time0)
+    PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+    int2_grad1_u12_squared_ao = 0.d0
+    !$OMP PARALLEL               &
+    !$OMP DEFAULT (NONE)         &
+    !$OMP PRIVATE (i, j, ipoint) &
+    !$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
+    !$OMP DO SCHEDULE (static)
+    do ipoint = 1, n_points_final_grid
+      do j = 1, ao_num
+        do i = 1, ao_num
+          int2_grad1_u12_squared_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+    call wall_time(time1)
+    print*, ' wall time for int2_grad1_u12_squared_ao =', time1-time0 
+
+    ! ---
+
+  elseif(j1b_type .ge. 100) then
+
+    ! ---
+
+    call wall_time(time0)
+    PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
+    PROVIDE grad1_u12_squared_num
+    PROVIDE grad1_u12_num
+    double precision, allocatable :: tmp(:,:,:)
+    allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
+    tmp = 0.d0
+    !$OMP PARALLEL               &
+    !$OMP DEFAULT (NONE)         &
+    !$OMP PRIVATE (j, i, jpoint) &
+    !$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
+    !$OMP DO SCHEDULE (static)
+    do j = 1, ao_num
+      do i = 1, ao_num
+        do jpoint = 1, n_points_extra_final_grid
+          tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+    if(.not.read_tc_integ) then
+      int2_grad1_u12_ao = 0.d0
+      do m = 1, 3
+        call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, 1.d0          &
+                  , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid &
+                  , 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num) 
+      enddo
+
+      !! DEBUG
+      !PROVIDE v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
+      !int2_grad1_u12_ao = 0.d0
+      !do ipoint = 1, n_points_final_grid
+      !  x = final_grid_points(1,ipoint)
+      !  y = final_grid_points(2,ipoint)
+      !  z = final_grid_points(3,ipoint)
+      !  tmp0  = 0.5d0 * v_1b(ipoint)
+      !  tmp_x =  v_1b_grad(1,ipoint)
+      !  tmp_y =  v_1b_grad(2,ipoint)
+      !  tmp_z =  v_1b_grad(3,ipoint)
+      !  do j = 1, ao_num
+      !    do i = 1, ao_num
+      !      tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
+      !      tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
+      !      int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
+      !      int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
+      !      int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
+      !    enddo
+      !  enddo
+      !enddo
+
+      ! these dgemm are equivalen to
+      !!$OMP PARALLEL                                                           &
+      !!$OMP DEFAULT (NONE)                                                     &
+      !!$OMP PRIVATE (j, i, ipoint, jpoint, w)                                  &
+      !!$OMP SHARED (int2_grad1_u12_ao, ao_num, n_points_final_grid,            &
+      !!$OMP         n_points_extra_final_grid, final_weight_at_r_vector_extra, &
+      !!$OMP         aos_in_r_array_extra_transp, grad1_u12_num, tmp)
+      !!$OMP DO SCHEDULE (static)
+      !do ipoint = 1, n_points_final_grid
+      !  do j = 1, ao_num
+      !    do i = 1, ao_num
+      !      do jpoint = 1, n_points_extra_final_grid
+      !        w = tmp(jpoint,i,j)
+      !        int2_grad1_u12_ao(i,j,ipoint,1) += w * grad1_u12_num(jpoint,ipoint,1)
+      !        int2_grad1_u12_ao(i,j,ipoint,2) += w * grad1_u12_num(jpoint,ipoint,2)
+      !        int2_grad1_u12_ao(i,j,ipoint,3) += w * grad1_u12_num(jpoint,ipoint,3)
+      !      enddo
+      !    enddo
+      !  enddo
+      !enddo
+      !!$OMP END DO
+      !!$OMP END PARALLEL
+      call wall_time(time1)
+      print*, ' wall time for int2_grad1_u12_ao =', time1-time0 
+    endif
+
+    ! ---
+
+    call wall_time(time0)
+    int2_grad1_u12_squared_ao = 0.d0
+    call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0              &
+              , tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid &
+              , 0.d0, int2_grad1_u12_squared_ao(1,1,1), ao_num*ao_num) 
+
+    !! DEBUG
+    !PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+    !int2_grad1_u12_squared_ao = 0.d0
+    !!$OMP PARALLEL               &
+    !!$OMP DEFAULT (NONE)         &
+    !!$OMP PRIVATE (i, j, ipoint) &
+    !!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
+    !!$OMP DO SCHEDULE (static)
+    !do ipoint = 1, n_points_final_grid
+    !  do j = 1, ao_num
+    !    do i = 1, ao_num
+    !      int2_grad1_u12_squared_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+    !    enddo
+    !  enddo
+    !enddo
+    !!$OMP END DO
+    !!$OMP END PARALLEL
+    !call wall_time(time1)
+
+    !! this dgemm is equivalen to
+    !!$OMP PARALLEL                                                           &
+    !!$OMP DEFAULT (NONE)                                                     &
+    !!$OMP PRIVATE (i, j, ipoint, jpoint, w)                                  &
+    !!$OMP SHARED (int2_grad1_u12_squared_ao, ao_num, n_points_final_grid,    &
+    !!$OMP         n_points_extra_final_grid, final_weight_at_r_vector_extra, &
+    !!$OMP         aos_in_r_array_extra_transp, grad1_u12_squared_num, tmp)
+    !!$OMP DO SCHEDULE (static)
+    !do ipoint = 1, n_points_final_grid
+    !  do j = 1, ao_num
+    !    do i = 1, ao_num
+    !      do jpoint = 1, n_points_extra_final_grid
+    !        w = 0.5d0 * tmp(jpoint,i,j)
+    !        int2_grad1_u12_squared_ao(i,j,ipoint) += w * grad1_u12_squared_num(jpoint,ipoint)
+    !      enddo
+    !    enddo
+    !  enddo
+    !enddo
+    !!$OMP END DO
+    !!$OMP END PARALLEL
+    call wall_time(time1)
+    print*, ' wall time for int2_grad1_u12_squared_ao =', time1-time0 
+
+    ! ---
+
+    deallocate(tmp)
+  
+  else
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+  endif
+
+  if(write_tc_integ.and.mpi_master) then
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write")
+    call ezfio_set_work_empty(.False.)
+    write(11) int2_grad1_u12_ao
+    close(11)
+    call ezfio_set_tc_keywords_io_tc_integ('Read')
+  endif
+
+END_PROVIDER
+
+! ---
+