diff --git a/plugins/local/jastrow/EZFIO.cfg b/plugins/local/jastrow/EZFIO.cfg
index 8fd2d05a..4fd833cf 100644
--- a/plugins/local/jastrow/EZFIO.cfg
+++ b/plugins/local/jastrow/EZFIO.cfg
@@ -1,3 +1,14 @@
+[mu_of_r_tc]
+type: logical
+doc: If |true|, take the second formula for mu(r)
+interface: ezfio,provider,ocaml
+default: False
+
+[mu_of_r_av]
+type: logical
+doc: If |true|, take the second formula for mu(r)
+interface: ezfio,provider,ocaml
+default: False
 
 [j2e_type]
 type: character*(32)
diff --git a/plugins/local/non_h_ints_mu/NEED b/plugins/local/non_h_ints_mu/NEED
index 5ca1d543..bfc4f311 100644
--- a/plugins/local/non_h_ints_mu/NEED
+++ b/plugins/local/non_h_ints_mu/NEED
@@ -4,3 +4,4 @@ jastrow
 ao_tc_eff_map
 bi_ortho_mos
 trexio
+mu_of_r
diff --git a/plugins/local/non_h_ints_mu/deb_j_mu_of_r b/plugins/local/non_h_ints_mu/deb_j_mu_of_r
new file mode 100755
index 00000000..f9728583
Binary files /dev/null and b/plugins/local/non_h_ints_mu/deb_j_mu_of_r differ
diff --git a/plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f b/plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f
new file mode 100644
index 00000000..480855e9
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/deb_j_mu_of_r.irp.f
@@ -0,0 +1,97 @@
+program test_j_mu_of_r
+ implicit none
+! call routine_test_mu_of_r
+ call routine_test_mu_of_r_tot
+end
+
+subroutine routine_test_mu_of_r_tot
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast,grad_jast_mu_r1(3)
+ double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3)
+
+ dr = 0.00001d0
+ r2 = 0.d0
+ r2(1) =  0.5d0
+ r2(2) = -0.1d0
+ r2(3) =  1.0d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast_mu_r1)
+   double precision :: norm
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call get_j_sum_mu_of_r(r1bis,r2,jast_p)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call get_j_sum_mu_of_r(r1bis,r2,jast_m)
+
+    num_grad_jast_mu_r1(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast_mu_r1(k)*num_grad_jast_mu_r1(k)
+   enddo
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-10)then
+    print*,'/////'
+    print*,grad_jast_mu_r1
+    print*,num_grad_jast_mu_r1
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+
+subroutine routine_test_mu_of_r
+ implicit none
+ integer :: ipoint,k 
+ double precision :: r2(3), weight, dr, r1(3), r1bis(3)
+ double precision :: accu_grad(3)
+ double precision :: jast_mu_r1,grad_jast_mu_r1(3)
+ double precision :: jast_mu_r1_p,jast_mu_r1_m,num_grad_jast_mu_r1(3)
+ double precision :: j12_mu_input
+ double precision :: mu_val_p, mu_val_m, dm_r1, mu_der(3), grad_dm_r1(3)
+
+ dr = 0.0001d0
+ r2 = 0.d0
+ r2(1) =  2.5d0
+ r2(2) =  2.25d0
+ r2(3) = -2.5d0
+ accu_grad = 0.d0
+  do ipoint = 1, n_points_final_grid
+   r1(1:3) = final_grid_points(1:3,ipoint)
+   weight = final_weight_at_r_vector(ipoint)
+   call grad_j_mu_of_r_1(r1,r2,jast_mu_r1,grad_jast_mu_r1, dm_r1, grad_dm_r1)
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call grad_mu_of_r_mean_field(r1bis,mu_val_p, dm_r1, mu_der, grad_dm_r1) 
+    jast_mu_r1_p = j12_mu_input(r1bis, r2, mu_val_p)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call grad_mu_of_r_mean_field(r1bis,mu_val_m, dm_r1, mu_der, grad_dm_r1) 
+    jast_mu_r1_m = j12_mu_input(r1bis, r2, mu_val_m)
+    num_grad_jast_mu_r1(k) = -(jast_mu_r1_p - jast_mu_r1_m)/(2.d0* dr)
+!    print*,jast_mu_r1_p,jast_mu_r1_m
+   enddo
+   print*,'/////'
+   print*,grad_jast_mu_r1
+   print*,num_grad_jast_mu_r1
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
diff --git a/plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f
new file mode 100644
index 00000000..d45c899a
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/jast_deriv_mu_of_r.irp.f
@@ -0,0 +1,164 @@
+subroutine get_j_sum_mu_of_r(r1,r2,jast)
+ implicit none
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: jast
+ double precision :: mu_r1, dm_r1, grad_mu_r1(3), grad_dm_r1(3), j_mu_r1
+ double precision :: mu_r2, dm_r2, grad_mu_r2(3), grad_dm_r2(3), j_mu_r2
+ double precision :: j12_mu_input
+ jast = 0.d0
+ call grad_mu_of_r_mean_field(r1,mu_r1, dm_r1, grad_mu_r1, grad_dm_r1) 
+ call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+ j_mu_r1 = j12_mu_input(r1, r2, mu_r1)
+ j_mu_r2 = j12_mu_input(r1, r2, mu_r2)
+ if(dm_r1 + dm_r2.lt.1.d-7)return
+ jast = (dm_r1 * j_mu_r1 + dm_r2 * j_mu_r2) / (dm_r1 + dm_r2)
+
+end
+
+subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
+ implicit none
+ BEGIN_DOC
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: jast, grad_jast(3)
+ jast = 0.d0
+ grad_jast = 0.d0
+ double precision :: num, denom, grad_num(3), grad_denom(3)
+ double precision :: j_r1, grad_j_r1(3),j_r2, grad_j_r2(3)
+ double precision :: dm_r1, grad_dm_r1(3), grad_jmu_r2(3)
+ double precision :: dm_r2, grad_dm_r2(3),mu_r2, grad_mu_r2(3)
+ double precision :: j12_mu_input
+
+ if(murho_type==0)then
+! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
+! 
+!          = num(r1,r2) / denom(r1,r2)
+!
+! d/dx1 J(r1,r2) = [denom(r1,r2) X d/dx1 num(r1,r2) - num(r1,r2) X d/dx1 denom(r1,r2) ] / denom(r1,r2)^2
+!
+! d/dx1 num(r1,r2) =  j(mu(r1),r12)*d/dx1 rho(r1) + rho(r1) * d/dx1 j(mu(r1),r12) 
+!                   + rho(r2) d/dx1 j(mu(r2),r12)
+! d/dx1 denom(r1,r2) = d/dx1 rho(r1)
+  call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1)
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2)
+  num = dm_r1 * j_r1 + dm_r2 * j_r2
+  denom = dm_r1 + dm_r2
+  if(denom.lt.1.d-7)return
+  jast = num / denom
+ 
+  grad_denom = grad_dm_r1
+ 
+  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2)
+  grad_num =  j_r1 * grad_dm_r1 + dm_r1 * grad_j_r1 + dm_r2 * grad_jmu_r2
+  
+  grad_jast = (grad_num * denom - num * grad_denom)/(denom*denom)
+ else if(murho_type==-1)then
+! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)] 
+!
+! d/dx1 J(r1,r2)  = 0.5 * (d/dx1 j(mu(r1),r12) + d/dx1 j(mu(r2),r12))
+  call grad_j_mu_of_r_1(r1,r2,j_r1, grad_j_r1,dm_r1, grad_dm_r1)
+  call grad_mu_of_r_mean_field(r2,mu_r2, dm_r2, grad_mu_r2, grad_dm_r2) 
+  j_r2 = j12_mu_input(r1, r2, mu_r2) ! j(mu(r2),r1,r2)
+  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2)
+  jast = 0.5d0 * (j_r1 + j_r2)
+  grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2)
+  
+ endif
+
+end
+
+subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+! grad_r1 of j(mu(r1),r12)
+  !
+  ! d/dx1 j(mu(r1),r12) = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 ) d/dx1 mu(r1) 
+  !                     + d/dx1 j(mu(r1),r12) 
+  !
+  ! with 
+  !
+  !           j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
+  !
+  ! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
+  !
+  ! here mu(r1) is obtained by MU MEAN FIELD 
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: jast, grad_jast(3),dm_r1, grad_dm_r1(3)
+ double precision :: dx, dy, dz,  r12, mu_der(3)
+ double precision :: mu_tmp, tmp, grad(3), mu_val
+ jast = 0.d0
+ 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)
+ ! get mu(r1) == mu_val and its gradient d/dx1 mu(r1) == mu_der 
+ call grad_mu_of_r_mean_field(r1,mu_val, dm_r1, mu_der, grad_dm_r1) 
+ mu_tmp  = mu_val * r12
+ ! evalulation of the jastrow j(mu(r1),r12)
+ jast = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_val
+
+ ! tmp = exp(-(mu(r1)*r12)**2) /(2 *sqrt(pi) * mu(r1)**2 )
+ tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
+ ! grad = 
+ grad(1) = tmp * mu_der(1)
+ grad(2) = tmp * mu_der(2)
+ grad(3) = tmp * mu_der(3)
+
+ if(r12 .lt. 1d-10) return
+ tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12) 
+ grad(1) = grad(1) + tmp * dx
+ grad(2) = grad(2) + tmp * dy
+ grad(3) = grad(3) + tmp * dz
+
+ grad_jast = grad 
+end
+
+! ---
+
+double precision function j12_mu_input(r1, r2, mu)
+
+  BEGIN_DOC
+  ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
+  END_DOC
+  include 'constants.include.F'
+
+  implicit none
+  double precision, intent(in) :: r1(3), r2(3), mu
+  double precision             :: mu_tmp, 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_tmp = mu * r12
+
+    j12_mu_input = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu
+
+end
+
+subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu)
+ implicit none
+ BEGIN_DOC
+ ! grad_jmu = 
+ END_DOC
+ double precision, intent(in) :: r1(3), r2(3), mu
+ double precision, intent(out):: grad_jmu(3)
+ double precision             :: mu_tmp, r12, dx, dy, dz, grad(3), tmp
+ grad_jmu = 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
+ mu_tmp  = mu * r12
+ tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12 ! d/dx1 j(mu(r1),r12) 
+ grad(1) = tmp * dx
+ grad(2) = tmp * dy
+ grad(3) = tmp * dz
+
+ grad_jmu = grad 
+end
diff --git a/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f
index 79822508..c3fc91a6 100644
--- a/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f
+++ b/plugins/local/non_h_ints_mu/jast_deriv_utils.irp.f
@@ -3,6 +3,9 @@
 
 double precision function j12_mu(r1, r2)
 
+  BEGIN_DOC
+  ! j(mu,r12) = 1/2 r12 (1 - erf(mu r12)) - 1/2 (sqrt(pi) * mu) e^{-(mu*r12)^2}
+  END_DOC
   include 'constants.include.F'
 
   implicit none
@@ -73,24 +76,41 @@ subroutine grad1_j12_mu(r1, r2, grad)
     grad(3) = tmp * dz
 
   elseif(j2e_type .eq. "Mur") then
+   double precision :: jast
+   call grad_j_sum_mu_of_r(r1,r2,jast,grad)
 
-    dx  = r1(1) - r2(1)
-    dy  = r1(2) - r2(2)
-    dz  = r1(3) - r2(3)
-    r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+!    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
 
-    call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
-    mu_tmp  = mu_val * r12
-    tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
-    grad(1) = tmp * mu_der(1)
-    grad(2) = tmp * mu_der(2)
-    grad(3) = tmp * mu_der(3)
-
-    if(r12 .lt. 1d-10) return
-    tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
-    grad(1) = grad(1) + tmp * dx
-    grad(2) = grad(2) + tmp * dy
-    grad(3) = grad(3) + tmp * dz
+!  elseif(j2e_type .eq. "Mur") then
+!
+!    dx  = r1(1) - r2(1)
+!    dy  = r1(2) - r2(2)
+!    dz  = r1(3) - r2(3)
+!    r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+!
+!    call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
+!    mu_tmp  = mu_val * r12
+!    tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
+!    grad(1) = tmp * mu_der(1)
+!    grad(2) = tmp * mu_der(2)
+!    grad(3) = tmp * mu_der(3)
+!
+!    if(r12 .lt. 1d-10) return
+!    tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
+!    grad(1) = grad(1) + tmp * dx
+!    grad(2) = grad(2) + tmp * dy
+!    grad(3) = grad(3) + tmp * dz
 
   else
 
@@ -369,7 +389,18 @@ end
 ! ---
 
 subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
-
+ BEGIN_DOC
+! various flavours of mu(r1,r2) 
+! depends on essentially the density and other related quantities 
+!
+! change the variable "murho_type" to change type
+!
+!  murho_type == -1 :: mu(r1,r2) = (rho(r1) mu_mf(r1) + rho(r2) mu_mf(r2))/[rho(r1)+rho(r2)]
+!
+!             ==  0 :: mu(r1,r2) = (sqrt(rho(r1)) mu_mf(r1) + sqrt(rho(r2)) mu_mf(r2))/[sqrt(rho(r1))+sqrt(rho(r2))]
+!             
+!             == -2 :: mu(r1,r2) = 0.5(mu_mf(r1) + mu_mf(r2))
+ END_DOC
   implicit none
   double precision, intent(in)  :: r1(3), r2(3)
   double precision, intent(out) :: mu_val, mu_der(3)
@@ -379,11 +410,50 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
   double precision :: rho1, grad_rho1(3),rho2,rho_tot,inv_rho_tot
   double precision :: f_rho1, f_rho2, d_drho_f_rho1
   double precision :: d_dx1_f_rho1(3),d_dx_rho_f_rho(3),nume
+  double precision :: mu_mf_r1, dm_r1, grad_mu_mf_r1(3), grad_dm_r1(3)
+  double precision :: mu_mf_r2, dm_r2, grad_mu_mf_r2(3), grad_dm_r2(3)
+
+  double precision :: num, denom, grad_denom(3), grad_num(3)
+  double precision :: dsqrt_dm_r1
 
   PROVIDE murho_type
   PROVIDE mu_r_ct mu_erf
 
-  if(murho_type .eq. 1) then
+  if(murho_type .eq. 0) then 
+   call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
+   call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
+   dsqrt_dm_r1 = dsqrt(dm_r1)
+   denom = (dsqrt_dm_r1 + dsqrt(dm_r2) )
+   if(denom.lt.1.d-7)then
+    mu_val = 1.d+10
+    mu_der = 0.d0
+    return
+   endif
+   num = (dsqrt(dm_r1) * mu_mf_r1 + dsqrt(dm_r2) * mu_mf_r2) 
+   mu_val = num / denom
+   grad_denom = grad_dm_r1/dsqrt_dm_r1
+   grad_num = dsqrt(dm_r1) * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1
+   mu_der = (grad_num * denom - num * grad_denom)/(denom*denom)
+  else if(murho_type .eq. -1) then
+   call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
+   call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
+   denom = (dm_r1 + dm_r2 )
+   if(denom.lt.1.d-7)then
+    mu_val = 1.d+10
+    mu_der = 0.d0
+    return
+   endif
+   num = (dm_r1 * mu_mf_r1 + dm_r2 * mu_mf_r2) 
+   mu_val = num / denom
+   grad_denom = grad_dm_r1
+   grad_num = dm_r1 * grad_mu_mf_r1 + mu_mf_r1 * grad_dm_r1
+   mu_der = (grad_num * denom - num * grad_denom)/(denom*denom)
+  else if(murho_type .eq. -2) then
+   call grad_mu_of_r_mean_field(r1,mu_mf_r1, dm_r1, grad_mu_mf_r1, grad_dm_r1)
+   call grad_mu_of_r_mean_field(r2,mu_mf_r2, dm_r2, grad_mu_mf_r2, grad_dm_r2)
+   mu_val = 0.5d0 * (mu_mf_r1 + mu_mf_r2) 
+   mu_der = 0.5d0 * grad_mu_mf_r1
+  else if(murho_type .eq. 1) then
 
     !
     ! r = 0.5 (r1 + r2)
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 2c41b535..a1f5c504 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
@@ -216,31 +216,35 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       r2(1) = final_grid_points_extra(1,jpoint)
       r2(2) = final_grid_points_extra(2,jpoint)
       r2(3) = final_grid_points_extra(3,jpoint)
-
-      dx  = r1(1) - r2(1)
-      dy  = r1(2) - r2(2)
-      dz  = r1(3) - r2(3)
-      r12 = dsqrt(dx * dx + dy * dy + dz * dz)
-
-      call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
-      mu_tmp  = mu_val * r12
-      tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
-      gradx(jpoint) = tmp * mu_der(1)
-      grady(jpoint) = tmp * mu_der(2)
-      gradz(jpoint) = tmp * mu_der(3)
-
-      if(r12 .lt. 1d-10) then
-        gradx(jpoint) = 0.d0
-        grady(jpoint) = 0.d0
-        gradz(jpoint) = 0.d0
-        cycle
-      endif
-
-      tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
-
-      gradx(jpoint) = gradx(jpoint) + tmp * dx
-      grady(jpoint) = grady(jpoint) + tmp * dy
-      gradz(jpoint) = gradz(jpoint) + tmp * dz
+      double precision :: jast, grad_jast(3)
+      call grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
+      gradx(jpoint) = grad_jast(1)
+      grady(jpoint) = grad_jast(2)
+      gradz(jpoint) = grad_jast(3)
+!      dx  = r1(1) - r2(1)
+!      dy  = r1(2) - r2(2)
+!      dz  = r1(3) - r2(3)
+!      r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+!
+!      call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
+!      mu_tmp  = mu_val * r12
+!      tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
+!      gradx(jpoint) = tmp * mu_der(1)
+!      grady(jpoint) = tmp * mu_der(2)
+!      gradz(jpoint) = tmp * mu_der(3)
+!
+!      if(r12 .lt. 1d-10) then
+!        gradx(jpoint) = 0.d0
+!        grady(jpoint) = 0.d0
+!        gradz(jpoint) = 0.d0
+!        cycle
+!      endif
+!
+!      tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
+!
+!      gradx(jpoint) = gradx(jpoint) + tmp * dx
+!      grady(jpoint) = grady(jpoint) + tmp * dy
+!      gradz(jpoint) = gradz(jpoint) + tmp * dz
     enddo
 
   elseif(j2e_type .eq. "Boys") then
diff --git a/plugins/local/non_h_ints_mu/tc_integ.irp.f b/plugins/local/non_h_ints_mu/tc_integ.irp.f
index 58e3db48..23aae519 100644
--- a/plugins/local/non_h_ints_mu/tc_integ.irp.f
+++ b/plugins/local/non_h_ints_mu/tc_integ.irp.f
@@ -127,8 +127,8 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
         ! TODO combine 1shot & int2_grad1_u12_ao_num
         PROVIDE int2_grad1_u12_ao_num
         int2_grad1_u12_ao = int2_grad1_u12_ao_num
-        !PROVIDE int2_grad1_u12_ao_num_1shot
-        !int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
+!         PROVIDE int2_grad1_u12_ao_num_1shot
+!         int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
       endif
 
     elseif(tc_integ_type .eq. "semi-analytic") then
diff --git a/plugins/local/tc_keywords/EZFIO.cfg b/plugins/local/tc_keywords/EZFIO.cfg
index f3bd75c8..b858fa5b 100644
--- a/plugins/local/tc_keywords/EZFIO.cfg
+++ b/plugins/local/tc_keywords/EZFIO.cfg
@@ -230,7 +230,7 @@ default: 70
 type: character*(32)
 doc: approach used to evaluate TC integrals [ analytic | numeric | semi-analytic ]
 interface: ezfio,ocaml,provider
-default: semi-analytic
+default: numeric
 
 [minimize_lr_angles]
 type: logical
diff --git a/src/mu_of_r/mu_of_r_mean_field.irp.f b/src/mu_of_r/mu_of_r_mean_field.irp.f
index 6abc7e4f..b55e4f6f 100644
--- a/src/mu_of_r/mu_of_r_mean_field.irp.f
+++ b/src/mu_of_r/mu_of_r_mean_field.irp.f
@@ -146,26 +146,18 @@ subroutine mu_of_r_mean_field(r,mu_mf, dm)
  endif
 end
 
-subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
+subroutine mu_of_r_mean_field_tc(r,mu_mf, dm)
  implicit none
- include 'constants.include.F'
+  include 'constants.include.F'
  double precision, intent(in) :: r(3)
- double precision, intent(out):: grad_mu_mf(3), grad_dm(3)
  double precision, intent(out):: mu_mf, dm
- double precision :: grad_f_mf_ab(3), grad_two_bod_dens(3),grad_dm_a(3), grad_dm_b(3)
  double precision :: f_mf_ab,two_bod_dens, dm_a, dm_b
- call get_grad_f_mf_ab(r,grad_f_mf_ab, grad_two_bod_dens,f_mf_ab,two_bod_dens, dm_a, dm_b,grad_dm_a, grad_dm_b)
-
+ call get_f_mf_ab(r,f_mf_ab,two_bod_dens, dm_a, dm_b)
  dm = dm_a + dm_b
- grad_dm(1:3) = grad_dm_a(1:3) + grad_dm_b(1:3)
-
  if(dabs(two_bod_dens).lt.1.d-10)then
   mu_mf = 1.d+10
-  grad_mu_mf = 0.d0
  else
-  mu_mf = 0.5d0 * sqpi * f_mf_ab/two_bod_dens
-  grad_mu_mf(1:3) = 0.5d0 * sqpi * (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& 
-                                   /(two_bod_dens*two_bod_dens)
- endif 
-
+  mu_mf = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0)
+ endif
 end
+
diff --git a/src/mu_of_r/test_proj_op.irp.f b/src/mu_of_r/test_proj_op.irp.f
index cf53c772..b03d958c 100644
--- a/src/mu_of_r/test_proj_op.irp.f
+++ b/src/mu_of_r/test_proj_op.irp.f
@@ -20,7 +20,8 @@ program projected_operators
 ! call test
 !  call test_f_mean_field
 ! call test_grad_f_mean_field
- call test_grad_mu_mf
+! call test_grad_mu_mf
+ call plot_mu_of_r_mf
 end
 
 
@@ -174,3 +175,33 @@ subroutine test_grad_mu_mf
  print*, accu_grad_mu_mf
 
 end
+
+subroutine plot_mu_of_r_mf
+ implicit none
+  include 'constants.include.F'
+ integer :: ipoint,npoint
+ double precision :: dx,r(3),xmax,xmin
+ double precision :: accu_mu,accu_nelec,mu_mf, dm,mu_mf_tc
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen
+ output=trim(ezfio_filename)//'.mu_mf'
+ i_unit_output = getUnitAndOpen(output,'w')
+ xmax = 5.D0
+ xmin = 0.d0
+ npoint = 10000
+ dx = (xmax - xmin)/dble(npoint)
+ r = 0.d0
+ r(1) = xmin
+ accu_mu = 0.d0
+ accu_nelec = 0.d0
+ do ipoint = 1, npoint
+  call mu_of_r_mean_field(r,mu_mf, dm)
+  call mu_of_r_mean_field_tc(r,mu_mf_tc, dm)
+  write(i_unit_output,'(100(F16.10,X))')r(1),mu_mf,mu_mf_tc,dm
+  accu_mu    += mu_mf * dm * r(1)**2*dx*4.D0*pi
+  accu_nelec +=         dm * r(1)**2*dx*4.D0*pi
+  r(1) += dx
+ enddo
+ print*,'nelec      = ',accu_nelec
+ print*,'mu average = ',accu_mu/accu_nelec
+end