diff --git a/plugins/local/jastrow/EZFIO.cfg b/plugins/local/jastrow/EZFIO.cfg
index 4fd833cf..1b616f02 100644
--- a/plugins/local/jastrow/EZFIO.cfg
+++ b/plugins/local/jastrow/EZFIO.cfg
@@ -1,8 +1,8 @@
 [mu_of_r_tc]
-type: logical
-doc: If |true|, take the second formula for mu(r)
+type: character*(32)
+doc: type of the mu(r): [ Standard | Erfmu | Erfmugauss ]
 interface: ezfio,provider,ocaml
-default: False
+default: Standard
 
 [mu_of_r_av]
 type: logical
@@ -12,7 +12,7 @@ default: False
 
 [j2e_type]
 type: character*(32)
-doc: type of the 2e-Jastrow: [ None | Mu | Mu_Nu | Mur | Boys | Boys_Handy | Qmckl ]
+doc: type of the 2e-Jastrow: [ None | Mu | Mugauss | Mu_Nu | Mur | Murgauss | Bump | Boys | Boys_Handy | Qmckl ]
 interface: ezfio,provider,ocaml
 default: Mu
 
diff --git a/plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f b/plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f
new file mode 100644
index 00000000..1c8d7198
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/deb_deriv_mu.irp.f
@@ -0,0 +1,28 @@
+program test_j_mu_of_r
+ implicit none
+ double precision :: x,mu_min,dmu,mu_max, mu, mu_p, mu_m
+ double precision :: j_simple,j_p, j_m,numeric_d_mu,d_dx_mu
+ double precision :: accu
+ integer :: npt,i
+ npt = 1000
+ mu_min = 0.3d0
+ mu_max = 10.d0
+ dmu = (mu_max - mu_min)/dble(npt)
+ x = 0.7d0
+ mu = mu_min
+ do i = 1, npt
+  call get_deriv_mu_j12(x,mu,d_dx_mu)
+  mu_p = mu + dmu
+  mu_m = mu - dmu
+  j_p = j_simple(x,mu_p)
+  j_m = j_simple(x,mu_m)
+  numeric_d_mu = 0.5d0 * (j_p - j_m)/dmu
+  print*,mu
+  print*,numeric_d_mu,d_dx_mu,dabs(d_dx_mu-numeric_d_mu)
+  accu += dabs(d_dx_mu-numeric_d_mu)
+  mu += dmu
+ enddo
+ accu *= dmu
+ print*,'accu = ',accu
+end
+
diff --git a/plugins/local/non_h_ints_mu/deb_j_bump.irp.f b/plugins/local/non_h_ints_mu/deb_j_bump.irp.f
new file mode 100644
index 00000000..82872357
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/deb_j_bump.irp.f
@@ -0,0 +1,98 @@
+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),j_bump
+ 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)
+   call get_grad_j_bump_mu_of_r(r1,r2,grad_jast_mu_r1)
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    jast_p = j_bump(r1bis,r2,a_boys)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    jast_m = j_bump(r1bis,r2,a_boys)
+
+    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
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.dr)then
+     print*,'/////'
+     print*,error,norm
+     print*,grad_jast_mu_r1
+     print*,num_grad_jast_mu_r1
+    endif
+   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 :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3)
+ double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3)
+ double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp
+ double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm
+
+ dr = 0.0001d0
+ 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_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_r1, grad_dm_r1) 
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1) 
+
+    r1bis= r1
+    r1bis(k) -= dr
+    call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1) 
+    num_grad_mu_r1(k) = (mu_rp - mu_rm)/(2.d0* dr)
+!    print*,jast_mu_r1_p,jast_mu_r1_m
+   enddo
+   print*,'/////'
+   print*,mu_der_r1
+   print*,num_grad_mu_r1
+   do k = 1,3
+    accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
diff --git a/plugins/local/non_h_ints_mu/deb_j_gauss.irp.f b/plugins/local/non_h_ints_mu/deb_j_gauss.irp.f
new file mode 100644
index 00000000..264a6f04
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/deb_j_gauss.irp.f
@@ -0,0 +1,62 @@
+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(3),j_bump,j12_mu
+ double precision :: jast_p,jast_m,num_grad_jast(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 grad1_j12_mu(r1, r2, grad_jast)
+   grad_jast = - grad_jast
+   double precision :: norm,error
+   norm = 0.D0
+   do k = 1, 3
+    r1bis= r1
+    r1bis(k) += dr
+    jast_p = j12_mu(r1bis, r2)
+
+    r1bis= r1
+    r1bis(k) -= dr
+    jast_m = j12_mu(r1bis, r2)
+
+    num_grad_jast(k) = (jast_p - jast_m)/(2.d0* dr)
+    norm += num_grad_jast(k)*num_grad_jast(k)
+   enddo
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+   error *= 0.33333333d0
+   norm = dsqrt(norm)
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.dr)then
+     print*,'/////'
+     print*,error,norm
+     print*,grad_jast
+     print*,num_grad_jast
+    endif
+   endif
+   do k = 1,3
+    accu_grad(k) += weight * dabs(grad_jast(k) - num_grad_jast(k))
+   enddo
+  enddo
+ print*,'accu_grad = '
+ print*, accu_grad
+
+end
+
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
deleted file mode 100755
index f9728583..00000000
Binary files a/plugins/local/non_h_ints_mu/deb_j_mu_of_r and /dev/null 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
index 480855e9..fc29bf50 100644
--- 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
@@ -12,7 +12,7 @@ subroutine routine_test_mu_of_r_tot
  double precision :: jast,grad_jast_mu_r1(3)
  double precision :: jast_p,jast_m,num_grad_jast_mu_r1(3)
 
- dr = 0.00001d0
+ dr = 0.000001d0
  r2 = 0.d0
  r2(1) =  0.5d0
  r2(2) = -0.1d0
@@ -22,7 +22,7 @@ subroutine routine_test_mu_of_r_tot
    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
+   double precision :: norm,error
    norm = 0.D0
    do k = 1, 3
     r1bis= r1
@@ -36,11 +36,19 @@ subroutine routine_test_mu_of_r_tot
     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
+   error = 0.d0
+   do k = 1, 3
+    error += dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+   enddo
+   error *= 0.33333333d0
    norm = dsqrt(norm)
-   if(norm.gt.1.d-10)then
-    print*,'/////'
-    print*,grad_jast_mu_r1
-    print*,num_grad_jast_mu_r1
+   if(norm.gt.1.d-05)then
+    if(dabs(error/norm).gt.10.d0*dr)then
+     print*,'/////'
+     print*,error,norm,dabs(error/norm)
+     print*,grad_jast_mu_r1
+     print*,num_grad_jast_mu_r1
+    endif
    endif
    do k = 1,3
     accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
@@ -54,41 +62,33 @@ 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)
+ double precision :: weight, dr, r1(3), r1bis(3),accu_grad(3),num_grad_mu_r1(3)
+ double precision :: mu_r1,dm_r1, mu_der_r1(3), grad_dm_r1(3)
+ double precision :: mu_der_rp(3), grad_dm_rp(3),mu_rp
+ double precision :: mu_der_rm(3), grad_dm_rm(3),mu_rm
 
  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)
+   call grad_mu_of_r_mean_field(r1,mu_r1,dm_r1, mu_der_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)
+    call grad_mu_of_r_mean_field(r1bis,mu_rp, dm_r1, mu_der_rp, grad_dm_r1) 
 
     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)
+    call grad_mu_of_r_mean_field(r1bis,mu_rm, dm_r1, mu_der_rm, grad_dm_r1) 
+    num_grad_mu_r1(k) = (mu_rp - mu_rm)/(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
+   print*,mu_der_r1
+   print*,num_grad_mu_r1
    do k = 1,3
-    accu_grad(k) += weight * dabs(grad_jast_mu_r1(k) - num_grad_jast_mu_r1(k))
+    accu_grad(k) += weight * dabs(mu_der_r1(k) - num_grad_mu_r1(k))
    enddo
   enddo
  print*,'accu_grad = '
diff --git a/plugins/local/non_h_ints_mu/j_bump.irp.f b/plugins/local/non_h_ints_mu/j_bump.irp.f
new file mode 100644
index 00000000..1731bb72
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/j_bump.irp.f
@@ -0,0 +1,90 @@
+double precision function wigner_radius(rho)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: rho
+ wigner_radius = 4.d0 * pi * rho * 0.333333333333d0
+ wigner_radius = wigner_radius**(-0.3333333d0)
+end
+
+double precision function j_bump(r1,r2,a)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: r1(3),r2(3),a
+ double precision :: inv_a,factor,x_scaled,scalar
+ double precision :: r12
+ r12   = (r1(1) - r2(1))*(r1(1) - r2(1))
+ r12  += (r1(2) - r2(2))*(r1(2) - r2(2))
+ r12  += (r1(3) - r2(3))*(r1(3) - r2(3))
+ r12 = dsqrt(r12)
+ inv_a = 1.d0/a
+ x_scaled = r12*inv_a*inv_sq_pi
+ x_scaled*= x_scaled
+ j_bump = 0.5d0 * (r12-a) * dexp(-x_scaled)
+end
+
+subroutine get_grad_j_bump(x,a,grad)
+ implicit none
+ BEGIN_DOC
+ ! gradient of the Jastrow with a bump
+ !
+ ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
+ !
+ ! d/dx j(x,a) = 1/(2 pi a^2) * exp[-(x/(a*sqrt(pi)))^2] * (pi a^2 + 2 a x - 2x^2)
+ END_DOC
+  include 'constants.include.F'
+ double precision, intent(in)  :: x,a
+ double precision, intent(out) :: grad
+ double precision :: inv_a,factor,x_scaled,scalar
+ inv_a = 1.d0/a
+ factor = 0.5d0*inv_pi*inv_a*inv_a 
+ x_scaled = x*inv_a*inv_sq_pi
+ x_scaled*= x_scaled
+ grad = factor * dexp(-x_scaled) * (pi*a*a + 2.d0 * a*x - 2.d0*x*x)
+end
+
+subroutine get_d_da_j_bump(x,a,d_da)
+ implicit none
+ BEGIN_DOC
+ ! Derivative with respect by to the parameter "a" of the Jastrow with a bump
+ !
+ ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
+ !
+ ! d/da j(x,a) = - 1/(pi*a^3) * exp[-(x/(a*sqrt(pi)))^2] * (-2 x^3 + 2 a x^2 + pi a^x3)
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in)  :: x,a
+ double precision, intent(out) :: d_da
+ double precision :: factor, inv_a,x_scaled,scalar
+ inv_a = 1.d0/a
+ factor = inv_a*inv_a*inv_a*inv_pi
+ x_scaled = x*inv_a*inv_sq_pi
+ x_scaled*= x_scaled
+ d_da = factor * dexp(-x_scaled) * (-2.d0 * x*x*x + 2.d0*x*x*a+pi*a*a*a)
+end
+
+subroutine get_grad_j_bump_mu_of_r(r1,r2,grad_j_bump)
+ implicit none
+ BEGIN_DOC
+ ! d/dx1 j(x,a(r1,r2)) where j(x,a) is the Jastrow with a bump
+ !
+ ! j(x,a) = 1/2 * (x-a)* exp[-(x/(a*sqrt(pi)))^2]
+ !
+ ! a(r1,r2) = [rho(r1) a(r1) + rho(r2) a(r2)]/[rho(r1) + rho(r2)]
+ !
+ ! d/dx1 j(x,a) = d/dx1 j(x,a(r1,r2))
+ END_DOC
+ double precision, intent(in) :: r1(3),r2(3)
+ double precision, intent(out):: grad_j_bump(3)
+ double precision :: r12,r12_vec(3),grad_scal,inv_r12
+ r12_vec = r1 - r2
+ r12   = (r1(1) - r2(1))*(r1(1) - r2(1))
+ r12  += (r1(2) - r2(2))*(r1(2) - r2(2))
+ r12  += (r1(3) - r2(3))*(r1(3) - r2(3))
+ r12 = dsqrt(r12)
+ call get_grad_j_bump(r12,a_boys,grad_scal)
+ if(r12.lt.1.d-10)then
+  grad_j_bump = 0.d0
+ else 
+  grad_j_bump = grad_scal * r12_vec/r12
+ endif
+end
diff --git a/plugins/local/non_h_ints_mu/jast_deriv.irp.f b/plugins/local/non_h_ints_mu/jast_deriv.irp.f
index 9a430135..1f97c18a 100644
--- a/plugins/local/non_h_ints_mu/jast_deriv.irp.f
+++ b/plugins/local/non_h_ints_mu/jast_deriv.irp.f
@@ -31,7 +31,7 @@
   grad1_u12_squared_num = 0.d0
 
   if( ((j2e_type .eq. "Mu") .and. (env_type .eq. "None")) .or. &
-       (j2e_type .eq. "Mur") ) then
+       (j2e_type .eq. "Mur").or.(j2e_type .eq. "Mugauss") .or. (j2e_type .eq. "Murgauss")) then
 
     !$OMP PARALLEL                                                                                    &
     !$OMP DEFAULT (NONE)                                                                              &
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
index d45c899a..8f5aee0c 100644
--- 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
@@ -4,19 +4,48 @@ subroutine get_j_sum_mu_of_r(r1,r2,jast)
  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
+ double precision :: j12_mu_input,mu_tot,r12,j_simple
  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)
+ if(murho_type==0)then
+! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
+  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)
+ else if(murho_type==1)then
+! J(r1,r2) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2))
+  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) 
+  mu_tot = 0.5d0 * (mu_r1 + mu_r2)
+  jast = j12_mu_input(r1, r2, mu_tot)
+ else if(murho_type==2)then
+! MU(r1,r2) = (rho(1) * mu(r1)+ rho(2) * mu(r2))/(rho(1)+rho(2))
+! J(r1,r2) = j(MU(r1,r2),r12)
+  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) 
+  double precision :: mu_tmp, dm_tot, dm_tot_inv
+  dm_tot = dm_r1**a_boys + dm_r2**a_boys  ! rho(1)**alpha+rho(2)**alpha
+  if(dm_tot.lt.1.d-12)then
+   dm_tot_inv = 1.d+12
+  else
+   dm_tot_inv = 1.d0/dm_tot
+  endif
+  mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2)
+  mu_tot = nu_erf * mu_tmp*dm_tot_inv ! 
+  r12  = (r1(1) - r2(1)) * (r1(1) - r2(1))
+  r12 += (r1(2) - r2(2)) * (r1(2) - r2(2))
+  r12 += (r1(3) - r2(3)) * (r1(3) - r2(3))
+  r12 = dsqrt(r12)
+  jast = j_simple(r12,mu_tot)
+ endif
 
 end
 
 subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
  implicit none
+  include 'constants.include.F'
  BEGIN_DOC
  END_DOC
  double precision, intent(in) :: r1(3),r2(3)
@@ -26,8 +55,25 @@ subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
  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
+ double precision :: dm_r2, grad_dm_r2(3),mu_r2, grad_mu_r2(3),mu_r1
+ double precision :: j12_mu_input,r12,grad_mu_r1(3),grad_jmu_r1(3)
+ double precision :: mu_tot,dx,dy,dz,r12_vec(3),d_dmu_j,d_dr12_j
+
+ 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.gt.1.d-10)then
+  r12_vec(1) = dx
+  r12_vec(2) = dy
+  r12_vec(3) = dz
+  r12_vec *= 1.d0/r12
+  ! r12_vec = grad_r1 (r12)
+ else
+  r12 = 1.d-10
+  r12_vec = 0.d0
+ endif
 
  if(murho_type==0)then
 ! J(r1,r2) = [rho(r1) * j(mu(r1),r12) + rho(r2) * j(mu(r2),r12)] / [rho(r1) + rho(r2)]
@@ -48,11 +94,65 @@ subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
   jast = num / denom
  
   grad_denom = grad_dm_r1
- 
-  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2)
+  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12)
   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) = j(0.5 * (mu(r1)+mu(r2)),r12), MU(r1,r2) = 0.5 *(mu(r1)+mu(r2))
+!
+! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst 
+!                + d/dMU [j(MU,r12)]
+!                x d/d(mu(r1)) MU(r1,r2)
+!                x d/dx1 mu(r1)
+!                = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1
+!                + e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2) 
+!                x 0.5 
+!                x d/dx1 mu(r1)
+ 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)
+ mu_tot = 0.5d0 * (mu_r1 + mu_r2)
+ call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12)
+ grad_jast = grad_jmu_r1
+ grad_jast+= dexp(-r12*mu_tot*r12*mu_tot) * inv_sq_pi_2 /(mu_tot* mu_tot) * 0.5d0 * grad_mu_r1
+ else if(murho_type==2)then
+! MU(r1,r2) = beta * (rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2))/(rho(1)**alpha+rho(2)**alpha)
+! J(r1,r2) = j(MU(r1,r2),r12)
+!
+! d/dx1 J(r1,r2) = d/dx1 j(MU(r1,r2),r12)|MU=cst 
+!                + d/dMU [j(MU,r12)] 
+!                x d/d(mu(r1)) MU(r1,r2)
+!                x d/dx1 mu(r1)
+!                = 0.5 * (1 - erf(MU(r1,r2) *r12))/r12 * (x1 - x2) == grad_jmu_r1
+!                + 0.5 e^{-(r12*MU(r1,r2))^2}/(2 sqrt(pi) * MU(r1,r2)^2) 
+!                x d/dx1 MU(r1,r2)
+! with d/dx1 MU(r1,r2) = beta * {[mu(1) d/dx1 [rho(1)**alpha] + rho(1)**alpha * d/dx1 mu(1)](rho(1)**alpha+rho(2)**alpha) 
+!                       - MU(1,2) d/dx1 [rho(1)]**alpha}/(rho(1)**alpha+rho(2)**alpha)^2
+! d/dx1 [rho(1)]**alpha = alpha [rho(1)]**(alpha-1) d/dx1 rho(1)
+!                        
+ 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)
+ double precision :: dm_tot,dm_tot_inv,grad_mu_tot(3),mu_tmp,grad_dm_r1_alpha(3),d_dx_j
+ dm_tot = dm_r1**a_boys + dm_r2**a_boys  ! rho(1)**alpha+rho(2)**alpha
+ grad_dm_r1_alpha = a_boys * dm_r1**(a_boys-1) * grad_dm_r1
+ if(dm_tot.lt.1.d-12)then
+  dm_tot_inv = 1.d+12
+ else
+  dm_tot_inv = 1.d0/dm_tot
+ endif
+ mu_tmp = dm_r1**a_boys * mu_r1 + dm_r2**a_boys * mu_r2 !rho(1)**alpha * mu(r1)+ rho(2)**alpha * mu(r2)
+ mu_tot = nu_erf * mu_tmp*dm_tot_inv ! 
+ grad_mu_tot = ( mu_r1 * grad_dm_r1_alpha + dm_r1**a_boys * grad_mu_r1 ) * dm_tot & 
+              -  mu_tmp * grad_dm_r1_alpha
+ grad_mu_tot *= dm_tot_inv * dm_tot_inv * nu_erf
+ call get_deriv_r12_j12(r12,mu_tot,d_dr12_j) ! d/dr12 j(MU(r1,r2,r12)
+ ! d/dx1 j(MU(r1,r2),r12) | MU(r1,r2) = cst
+ ! d/dr12 j(MU(r1,r2,r12) x d/dx1 r12
+ grad_jmu_r1 = d_dr12_j * r12_vec 
+! call grad_j12_mu_input(r1, r2, mu_tot, grad_jmu_r1,r12)
+ grad_jast = grad_jmu_r1
+ ! d/dMU j(MU(r1,r2),r12)
+ call get_deriv_mu_j12(r12,mu_tot,d_dmu_j)
+ grad_jast+= d_dmu_j * grad_mu_tot
  else if(murho_type==-1)then
 ! J(r1,r2) = 0.5 * [j(mu(r1),r12) + j(mu(r2),r12)] 
 !
@@ -60,7 +160,7 @@ subroutine grad_j_sum_mu_of_r(r1,r2,jast,grad_jast)
   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)
+  call grad_j12_mu_input(r1, r2, mu_r2, grad_jmu_r2,r12)
   jast = 0.5d0 * (j_r1 + j_r2)
   grad_jast = 0.5d0 * (grad_j_r1 + grad_jmu_r2)
   
@@ -73,6 +173,7 @@ subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
   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) 
@@ -83,6 +184,8 @@ subroutine grad_j_mu_of_r_1(r1,r2,jast, grad_jast, dm_r1, grad_dm_r1)
   !
   ! and d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
   !
+  !     d/d mu j(mu,r12) = e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2)
+  !
   ! here mu(r1) is obtained by MU MEAN FIELD 
  END_DOC
  double precision, intent(in) :: r1(3),r2(3)
@@ -140,14 +243,16 @@ double precision function j12_mu_input(r1, r2, mu)
 
 end
 
-subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu)
+subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu,r12)
  implicit none
  BEGIN_DOC
- ! grad_jmu = 
+ ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1)
+ !
+ !          = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)]
  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
+ double precision, intent(out):: grad_jmu(3),r12
+ double precision             :: mu_tmp, dx, dy, dz, grad(3), tmp
  grad_jmu = 0.d0
  dx  = r1(1) - r2(1)
  dy  = r1(2) - r2(2)
@@ -162,3 +267,40 @@ subroutine grad_j12_mu_input(r1, r2, mu, grad_jmu)
 
  grad_jmu = grad 
 end
+
+subroutine j12_and_grad_j12_mu_input(r1, r2, mu, jmu, grad_jmu)
+ implicit none
+ include 'constants.include.F'
+ BEGIN_DOC
+ ! jmu = j(mu,r12) 
+ ! grad_jmu = d/dx1 j(mu,r12) assuming mu=cst(r1)
+ !
+ !          = 0.5/r_12 * (x_1 - x_2) * [1 - erf(mu*r12)]
+ END_DOC
+ double precision, intent(in) :: r1(3), r2(3), mu
+ double precision, intent(out):: grad_jmu(3), jmu
+ double precision             :: mu_tmp, r12, dx, dy, dz, grad(3), tmp
+ double precision :: erfc_mur12,inv_mu
+ inv_mu = 1.d0/mu
+
+ grad_jmu = 0.d0 ! initialization when r12 --> 0
+ jmu = - inv_sq_pi_2 * inv_mu ! initialization when r12 --> 0
+
+ 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
+ erfc_mur12 = (1.d0 - derf(mu_tmp))
+ mu_tmp  = mu * r12
+ tmp     = 0.5d0 * erfc_mur12  / r12 ! d/dx1 j(mu(r1),r12) 
+ grad(1) = tmp * dx
+ grad(2) = tmp * dy
+ grad(3) = tmp * dz
+
+ grad_jmu = grad 
+
+ jmu= 0.5d0 * r12 * erfc_mur12 - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) * inv_mu
+
+
+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 c3fc91a6..f13990de 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
@@ -1,3 +1,65 @@
+subroutine get_deriv_r12_j12(x,mu,d_dx_j)
+ implicit none
+  include 'constants.include.F'
+ BEGIN_DOC
+ ! d/dr12 j(mu,r12)
+ END_DOC
+ double precision, intent(in)  :: x,mu
+ double precision, intent(out) :: d_dx_j
+ 
+ d_dx_j = 0.d0
+ if(x .lt. 1d-10) return
+ if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then
+  d_dx_j = 0.5d0 * (1.d0 - derf(mu * x))
+ else if(j2e_type .eq. "Mugauss" .or. j2e_type .eq. "Murgauss" ) then
+  double precision :: x_tmp
+  x_tmp = mu * x
+  ! gradient of j(mu,x)
+  d_dx_j = 0.5d0 * (1.d0 - derf(x_tmp))
+
+  ! gradient of gaussian additional term
+  x_tmp *= alpha_mu_gauss
+  x_tmp *= x_tmp 
+  d_dx_j += -0.5d0 * mu * c_mu_gauss * x * dexp(-x_tmp)
+ else
+   print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type
+   stop
+ endif
+end
+
+
+subroutine get_deriv_mu_j12(x,mu,d_d_mu)
+ implicit none
+ BEGIN_DOC
+ ! d/dmu j(mu,r12)
+ END_DOC
+  include 'constants.include.F'
+ double precision, intent(in)  :: x,mu
+ double precision, intent(out) :: d_d_mu
+ double precision :: x_tmp,inv_mu_2,inv_alpha_2
+ 
+ d_d_mu = 0.d0
+ if(x .lt. 1d-10) return
+ x_tmp = x*mu
+ if(mu.lt.1.d-10) return
+ inv_mu_2 = mu*mu
+ inv_mu_2 = 1.d0/inv_mu_2
+ if(j2e_type .eq. "Mu" .or. j2e_type .eq. "Mur") then
+  ! e^{-(r12*mu)^2}/(2 sqrt(pi) * mu^2)
+  d_d_mu = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2
+ else if(j2e_type .eq. "Mugauss"  .or. j2e_type .eq. "Murgauss" ) then
+  d_d_mu  = dexp(-x_tmp*x_tmp) * inv_sq_pi_2 * inv_mu_2
+  inv_alpha_2 = 1.d0/alpha_mu_gauss
+  inv_alpha_2 *= inv_alpha_2
+  x_tmp *= alpha_mu_gauss
+  x_tmp *= x_tmp
+  d_d_mu += -0.25d0 * c_mu_gauss*inv_alpha_2*dexp(-x_tmp) * (1.d0 + 2.d0 * x_tmp) * inv_mu_2
+ else
+   print *, ' Error in get_deriv_r12_j12: Unknown j2e_type = ', j2e_type
+   stop
+ endif
+end
+
 
 ! ---
 
@@ -21,6 +83,18 @@ double precision function j12_mu(r1, r2)
 
     j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf
 
+  else if(j2e_type .eq. "Mugauss") 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)) )
+    double precision :: r12_tmp
+    r12_tmp = mu_erf * r12
+
+    j12_mu = 0.5d0 * r12 * (1.d0 - derf(r12_tmp)) - inv_sq_pi_2 * dexp(-r12_tmp*r12_tmp) / mu_erf
+    r12_tmp *= alpha_mu_gauss
+    j12_mu += 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu_erf) * dexp(-r12_tmp*r12_tmp)
+
   else
 
     print *, ' Error in j12_mu: Unknown j2e_type = ', j2e_type
@@ -60,7 +134,7 @@ subroutine grad1_j12_mu(r1, r2, grad)
 
   grad = 0.d0
 
-  if(j2e_type .eq. "Mu") then
+  if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mugauss") then
 
     dx = r1(1) - r2(1)
     dy = r1(2) - r2(2)
@@ -69,48 +143,42 @@ subroutine grad1_j12_mu(r1, r2, grad)
     r12 = dsqrt(dx * dx + dy * dy + dz * dz)
     if(r12 .lt. 1d-10) return
 
-    tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
+    call get_deriv_r12_j12(r12,mu_erf,tmp)
+!    tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
 
     grad(1) = tmp * dx
     grad(2) = tmp * dy
     grad(3) = tmp * dz
+    grad *= 1.d0/r12
 
-  elseif(j2e_type .eq. "Mur") then
+  elseif(j2e_type .eq. "Mur" .or. j2e_type .eq. "Murgauss") then
    double precision :: jast
    call grad_j_sum_mu_of_r(r1,r2,jast,grad)
+ 
+  elseif(j2e_type .eq. "Bump") then
+   double precision ::grad_jast(3)
+   call get_grad_j_bump_mu_of_r(r1,r2,grad_jast)
+     dx = r1(1) - r2(1)
+     dy = r1(2) - r2(2)
+     dz = r1(3) - r2(3)
 
-!    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
+     r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+     if(r12 .lt. 1d-10) then
+       grad(1) = 0.d0
+       grad(2) = 0.d0
+       grad(3) = 0.d0
+       return
+     endif                                                                                                                                                                            
+
+     tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
+
+     grad(1) = 0.5d0 * tmp * dx
+     grad(2) = 0.5d0 * tmp * dy
+     grad(3) = 0.5d0 * tmp * dz
+     grad(1) += 0.5d0 * grad_jast(1)
+     grad(2) += 0.5d0 * grad_jast(2)
+     grad(3) += 0.5d0 * grad_jast(3)
 
-!  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
 
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 a1f5c504..2ce9fd07 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
@@ -33,8 +33,11 @@ subroutine get_grad1_u12_withsq_r1_seq(ipoint, n_grid2, resx, resy, resz, res)
   r1(2) = final_grid_points(2,ipoint)
   r1(3) = final_grid_points(3,ipoint)
 
-  if( (j2e_type .eq. "Mu")  .or. &
-      (j2e_type .eq. "Mur") .or. &
+  if( (j2e_type .eq. "Mu")        .or. &
+      (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Mugauss")   .or. &
+      (j2e_type .eq. "Murgauss")  .or. &
+      (j2e_type .eq. "Bump")      .or. &
       (j2e_type .eq. "Boys") ) then
 
     if(env_type .eq. "None") then
@@ -206,7 +209,43 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       gradz(jpoint) = tmp * dz
     enddo
 
-  elseif(j2e_type .eq. "Mur") then
+  else if(j2e_type .eq. "Mugauss") then
+
+    !  d/dx1 j(mu,r12) = 0.5 * [(1 - erf(mu * r12)) / r12 - mu*c*r12*exp(-(mu*alpha*r12)^2] * (x1 - x2)
+
+    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)
+
+      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) then
+        gradx(jpoint) = 0.d0 
+        grady(jpoint) = 0.d0 
+        gradz(jpoint) = 0.d0 
+        cycle
+      endif
+
+      double precision :: r12_tmp
+      r12_tmp = mu_erf * r12
+      ! gradient of j(mu,r12)
+      tmp = 0.5d0 * (1.d0 - derf(r12_tmp)) / r12
+      ! gradient of gaussian additional term
+      r12_tmp *= alpha_mu_gauss
+      r12_tmp *= r12_tmp 
+      tmp += -0.5d0 * mu_erf * c_mu_gauss * r12 * dexp(-r12_tmp)/r12
+
+      gradx(jpoint) = tmp * dx
+      grady(jpoint) = tmp * dy
+      gradz(jpoint) = tmp * dz
+    enddo
+
+  elseif(j2e_type .eq. "Mur".or.j2e_type .eq. "Murgauss") then
 
     ! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2) 
     !                        + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
@@ -221,30 +260,41 @@ subroutine grad1_j12_r1_seq(r1, n_grid2, gradx, grady, gradz)
       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. "Bump") then
+
+    ! d/dx1 jbump(r1,r2)
+
+    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)
+      call get_grad_j_bump_mu_of_r(r1,r2,grad_jast)
+
+      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) then
+        gradx(jpoint) = 0.d0 
+        grady(jpoint) = 0.d0 
+        gradz(jpoint) = 0.d0 
+        cycle
+      endif
+
+      tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
+
+      gradx(jpoint) = 0.5d0 * tmp * dx
+      grady(jpoint) = 0.5d0 * tmp * dy
+      gradz(jpoint) = 0.5d0 * tmp * dz
+      gradx(jpoint) += 0.5d0 * grad_jast(1)
+      grady(jpoint) += 0.5d0 * grad_jast(2)
+      gradz(jpoint) += 0.5d0 * grad_jast(3)
+!      gradx(jpoint) =  grad_jast(1)
+!      grady(jpoint) =  grad_jast(2)
+!      gradz(jpoint) =  grad_jast(3)
     enddo
 
   elseif(j2e_type .eq. "Boys") then
@@ -673,8 +723,11 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz)
   r1(2) = final_grid_points(2,ipoint)
   r1(3) = final_grid_points(3,ipoint)
 
-  if( (j2e_type .eq. "Mu")  .or. &
-      (j2e_type .eq. "Mur") .or. &
+  if( (j2e_type .eq. "Mu")        .or. &
+      (j2e_type .eq. "Mugauss")   .or. &
+      (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Murgauss")  .or. &
+      (j2e_type .eq. "Bump")      .or. &
       (j2e_type .eq. "Boys") ) then
 
     if(env_type .eq. "None") then
@@ -791,8 +844,11 @@ subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res)
   r1(2) = final_grid_points(2,ipoint)
   r1(3) = final_grid_points(3,ipoint)
 
-  if( (j2e_type .eq. "Mu")  .or. &
-      (j2e_type .eq. "Mur") .or. &
+  if( (j2e_type .eq. "Mu")        .or. &
+      (j2e_type .eq. "Mur")       .or. &
+      (j2e_type .eq. "Mugauss")  .or. &
+      (j2e_type .eq. "Murgauss")  .or. &
+      (j2e_type .eq. "Mugauss")   .or. &
       (j2e_type .eq. "Boys") ) then
 
     if(env_type .eq. "None") then
diff --git a/plugins/local/non_h_ints_mu/plot_j_gauss.irp.f b/plugins/local/non_h_ints_mu/plot_j_gauss.irp.f
new file mode 100644
index 00000000..a4030d8c
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/plot_j_gauss.irp.f
@@ -0,0 +1,59 @@
+program plot_j_gauss
+ implicit none
+ double precision :: xmin, xmax, x, dx
+ double precision :: mu_min, mu_max, mu, d_mu
+ double precision :: pot_j_gauss,j_mu_simple,j_gauss_simple,pot_j_mu
+ double precision, allocatable :: mu_tab(:),j_mu(:),j_mu_gauss(:)
+ double precision, allocatable :: w_mu(:), w_mu_gauss(:)
+
+ character*(128) :: output
+ integer :: getUnitAndOpen
+ integer :: i_unit_output_wee_gauss,i_unit_output_wee_mu
+ integer :: i_unit_output_j_gauss,i_unit_output_j_mu
+ output=trim(ezfio_filename)//'.w_ee_mu_gauss'
+ i_unit_output_wee_gauss = getUnitAndOpen(output,'w')
+ output=trim(ezfio_filename)//'.w_ee_mu'
+ i_unit_output_wee_mu = getUnitAndOpen(output,'w')
+ output=trim(ezfio_filename)//'.j_mu_gauss'
+ i_unit_output_j_gauss = getUnitAndOpen(output,'w')
+ output=trim(ezfio_filename)//'.j_mu'
+ i_unit_output_j_mu = getUnitAndOpen(output,'w')
+
+ integer :: npt, i, j, n_mu
+ n_mu = 3
+ allocate(mu_tab(n_mu),j_mu(n_mu),j_mu_gauss(n_mu),w_mu(n_mu), w_mu_gauss(n_mu))
+ mu_min = 0.5d0
+ mu_max = 2.d0
+ d_mu = (mu_max - mu_min)/dble(n_mu)
+ mu = mu_min
+ do i = 1, n_mu
+  mu_tab(i) = mu
+  print*,'mu = ',mu
+  mu += d_mu
+ enddo
+ mu_tab(1) = 0.9d0
+ mu_tab(2) = 0.95d0
+ mu_tab(3) = 1.d0
+
+ xmin = 0.01d0
+ xmax = 10.d0
+ npt = 1000
+ dx = (xmax - xmin)/dble(npt)
+ x = xmin
+ do i = 1, npt
+  do j = 1, n_mu
+   mu = mu_tab(j)
+   w_mu_gauss(j) = pot_j_gauss(x,mu)
+   w_mu(j) = pot_j_mu(x,mu)
+   j_mu(j) = j_mu_simple(x,mu)
+   j_mu_gauss(j) = j_gauss_simple(x,mu) + j_mu(j)
+  enddo
+  write(i_unit_output_wee_gauss,'(100(F16.10,X))')x,w_mu_gauss(:)
+  write(i_unit_output_wee_mu,'(100(F16.10,X))')x,w_mu(:)
+  write(i_unit_output_j_gauss,'(100(F16.10,X))')x,j_mu_gauss(:)
+  write(i_unit_output_j_mu,'(100(F16.10,X))')x,j_mu(:)
+  x += dx
+ enddo
+
+
+end
diff --git a/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f b/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f
index 3a5984bd..4a3ec0d5 100644
--- a/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f
+++ b/plugins/local/non_h_ints_mu/plot_mu_of_r.irp.f
@@ -16,15 +16,16 @@ subroutine routine_print
  integer :: ipoint,nx,i
  double precision :: xmax,xmin,r(3),dx,sigma
  double precision :: mu_val, mu_der(3),dm_a,dm_b,grad,grad_dm_a(3), grad_dm_b(3)
- xmax =  5.D0
- xmin = -5.D0
+ xmax =  3.9D0
+ xmin = -3.9D0
  nx = 10000
  dx = (xmax - xmin)/dble(nx)
  r = 0.d0
  r(1) = xmin
  do ipoint = 1, nx
-  call mu_r_val_and_grad(r, r, mu_val, mu_der)
-  call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
+!  call mu_r_val_and_grad(r, r, mu_val, mu_der)
+  call grad_mu_of_r_mean_field(r,mu_val, dm_a, mu_der, grad_dm_a)
+!  call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
   sigma = 0.d0
   do i = 1,3
    sigma += grad_dm_a(i)**2
@@ -32,7 +33,8 @@ subroutine routine_print
   sigma=dsqrt(sigma)
   grad = mu_der(1)**2 + mu_der(2)**2 + mu_der(3)**2 
   grad = dsqrt(grad)
-  write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a+dm_b,grad,sigma/dm_a
+  print*,r(1),mu_val
+  write(i_unit_output,'(100(F16.7,X))')r(1),mu_val,dm_a,grad,sigma/dm_a
   r(1) += dx
  enddo
 end
diff --git a/plugins/local/non_h_ints_mu/pot_j_gauss.irp.f b/plugins/local/non_h_ints_mu/pot_j_gauss.irp.f
new file mode 100644
index 00000000..f9a0a7bc
--- /dev/null
+++ b/plugins/local/non_h_ints_mu/pot_j_gauss.irp.f
@@ -0,0 +1,146 @@
+double precision function j_simple(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ double precision :: j_mu_simple,j_gauss_simple
+ if(j2e_type .eq. "Mu".or.j2e_type .eq. "Mur") then
+  j_simple = j_mu_simple(x,mu)
+ else if(j2e_type .eq. "Mugauss".or.j2e_type .eq. "Murgauss") then
+  j_simple = j_gauss_simple(x,mu) + j_mu_simple(x,mu)
+ endif
+end
+
+
+double precision function j_mu_simple(x,mu)
+ implicit none
+ double precision, intent(in):: x,mu
+ include 'constants.include.F'
+ BEGIN_DOC
+! j_mu(mu,x) = 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
+ END_DOC
+ j_mu_simple = 0.5d0 * x * (1.D0 - derf(mu*x)) - 0.5d0 * inv_sq_pi/mu *  dexp(-x*mu*x*mu)
+
+end
+
+double precision function j_gauss_simple(x,mu)
+ implicit none
+ double precision, intent(in):: x,mu
+ include 'constants.include.F'
+ BEGIN_DOC
+! j_mu(mu,x) = c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
+!      with c = 27/(8 sqrt(pi)), alpha=3/2
+ END_DOC
+ double precision :: x_tmp
+ x_tmp = alpha_mu_gauss * mu * x
+ j_gauss_simple = 0.25d0 * c_mu_gauss / (alpha_mu_gauss*alpha_mu_gauss*mu) * dexp(-x_tmp*x_tmp)
+
+end
+
+double precision function j_mu_deriv(x,mu)
+ implicit none
+ BEGIN_DOC
+! d/dx j_mu(mu,x) = d/dx 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
+!              = 0.5*(1 - erf(mu x))
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ j_mu_deriv = 0.5d0 * (1.d0 - derf(mu*x))
+end
+
+double precision function j_mu_deriv_2(x,mu)
+ implicit none
+ BEGIN_DOC
+! d^2/dx^2 j_mu(mu,x) = d^2/dx^2 0.5 x (1 - erf(mu x)) - 1/[2 sqrt(pi)mu] exp(-(x*mu)^2)
+!                  = -mu/sqrt(pi) * exp(-(mu x)^2)
+ END_DOC
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ j_mu_deriv_2 = - mu * inv_sq_pi * dexp(-x*mu*x*mu)
+end
+
+double precision function j_gauss_deriv(x,mu)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
+!      with c = 27/(8 sqrt(pi)), alpha=3/2
+!                    = -0.5 * mu *  c * x * exp(-(alpha * mu * x)^2)
+ END_DOC
+ double precision :: x_tmp
+ x_tmp = alpha_mu_gauss * mu * x
+ j_gauss_deriv = -0.5d0 * mu * c_mu_gauss * x * exp(-x_tmp*x_tmp)
+end
+
+double precision function j_gauss_deriv_2(x,mu)
+ implicit none
+ include 'constants.include.F'
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d/dx j_gauss(mu,x) = d/dx c/[4 alpha^2 mu] exp(-(alpha * mu * x)^2)
+!      with c = 27/(8 sqrt(pi)), alpha=3/2
+!                    = 0.5 * mu *  c * exp(-(alpha * mu * x)^2) * (2 (alpha*mu*x)^2 - 1)
+ END_DOC
+ double precision :: x_tmp
+ x_tmp = alpha_mu_gauss * mu * x
+ x_tmp = x_tmp * x_tmp
+ j_gauss_deriv_2 = 0.5d0 * mu * c_mu_gauss * exp(-x_tmp) * (2.d0*x_tmp - 1.d0)
+end
+
+double precision function j_erf_gauss_deriv(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d/dx (j_gauss(mu,x)+j_mu(mu,x)) 
+ END_DOC
+ double precision :: j_gauss_deriv,j_mu_deriv
+ j_erf_gauss_deriv = j_gauss_deriv(x,mu)+j_mu_deriv(x,mu)
+end
+
+double precision function j_erf_gauss_deriv_2(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+! d^2/dx^2 (j_gauss(mu,x)+j_mu(mu,x)) 
+ END_DOC
+ double precision :: j_gauss_deriv_2,j_mu_deriv_2
+ j_erf_gauss_deriv_2 = j_gauss_deriv_2(x,mu)+j_mu_deriv_2(x,mu)
+end
+
+
+double precision function pot_j_gauss(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+ ! effective scalar potential associated with the erf_gauss correlation factor
+ ! 
+ ! 1/x( 1 - 2 * d/dx j_erf_gauss(x,mu)) - d^2/dx^2 j_erf_gauss(x,mu)) - d/dx d/dx (j_erf_gauss(x,mu))^2
+ END_DOC
+ double precision :: j_erf_gauss_deriv_2,j_erf_gauss_deriv
+ double precision :: deriv_1, deriv_2
+ pot_j_gauss = 0.d0
+ if(x.ne.0.d0)then
+  deriv_1 = j_erf_gauss_deriv(x,mu)
+  deriv_2 = j_erf_gauss_deriv_2(x,mu)
+  pot_j_gauss = 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2
+ endif
+
+end
+
+double precision function pot_j_mu(x,mu)
+ implicit none
+ double precision, intent(in) :: x,mu
+ BEGIN_DOC
+ ! effective scalar potential associated with the correlation factor
+ ! 
+ ! 1/x( 1 - 2 * d/dx j_erf(x,mu)) - d^2/dx^2 j_erf(x,mu)) - d/dx d/dx (j_erf(x,mu))^2
+ END_DOC
+ double precision :: j_mu_deriv_2,j_mu_deriv
+ double precision :: deriv_1, deriv_2
+ pot_j_mu = 0.d0
+ if(x.ne.0.d0)then
+  deriv_1 = j_mu_deriv(x,mu)
+  deriv_2 = j_mu_deriv_2(x,mu)
+  pot_j_mu= 1.d0/x * (1.d0 - 2.d0 * deriv_1) - deriv_1 * deriv_1 - deriv_2
+ endif
+
+end
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 b55e4f6f..295d58c2 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
@@ -57,6 +57,9 @@ end
 
 subroutine 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)
  implicit none
+ BEGIN_DOC
+ ! gradient of mu(r) mean field, together with the gradient of the one- and two-body HF density. 
+ END_DOC
  double precision, intent(in)  :: r(3)
  double precision, intent(out) :: f_mf_ab, two_bod_dens
  double precision, intent(out) :: grad_two_bod_dens(3), grad_f_mf_ab(3)
diff --git a/src/mu_of_r/test_proj_op.irp.f b/src/mu_of_r/test_proj_op.irp.f
index b03d958c..9ddc2f21 100644
--- a/src/mu_of_r/test_proj_op.irp.f
+++ b/src/mu_of_r/test_proj_op.irp.f
@@ -19,9 +19,9 @@ program projected_operators
 !  call test_f_ii_ia_aa_valence_ab
 ! call test
 !  call test_f_mean_field
-! call test_grad_f_mean_field
-! call test_grad_mu_mf
- call plot_mu_of_r_mf
+ call test_grad_f_mean_field
+ call test_grad_mu_mf
+! call plot_mu_of_r_mf
 end
 
 
diff --git a/src/utils/constants.include.F b/src/utils/constants.include.F
index 7b01f888..830b71a1 100644
--- a/src/utils/constants.include.F
+++ b/src/utils/constants.include.F
@@ -9,6 +9,9 @@ double precision, parameter :: pi_5_2 =  34.9868366552d0
 double precision, parameter :: dfour_pi =  4.d0*dacos(-1.d0)
 double precision, parameter :: dtwo_pi =  2.d0*dacos(-1.d0)
 double precision, parameter :: inv_sq_pi =  1.d0/dsqrt(dacos(-1.d0))
+double precision, parameter :: c_mu_gauss =  27.d0/(8.d0*dsqrt(dacos(-1.d0)))
+double precision, parameter :: c_mu_gauss_tot =  1.5d0*27.d0/(8.d0*dsqrt(dacos(-1.d0)))+3.d0/dsqrt(dacos(-1.d0))
+double precision, parameter :: alpha_mu_gauss =  1.5d0
 double precision, parameter :: inv_sq_pi_2 = 0.5d0/dsqrt(dacos(-1.d0))
 double precision, parameter :: thresh = 1.d-15
 double precision, parameter :: cx_lda = -0.73855876638202234d0