diff --git a/src/mo_basis/EZFIO.cfg b/src/mo_basis/EZFIO.cfg
index 81ffba5c..4c4f1eca 100644
--- a/src/mo_basis/EZFIO.cfg
+++ b/src/mo_basis/EZFIO.cfg
@@ -9,6 +9,12 @@ doc: Coefficient of the i-th |AO| on the j-th |MO|
 interface: ezfio
 size: (ao_basis.ao_num,mo_basis.mo_num)
 
+[mo_coef_aux]
+type: double precision
+doc: AUX Coefficient of the i-th |AO| on the j-th |MO|
+interface: ezfio
+size: (ao_basis.ao_num,mo_basis.mo_num)
+
 [mo_coef_imag]
 type: double precision
 doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO|
diff --git a/src/mo_basis/mos_aux.irp.f b/src/mo_basis/mos_aux.irp.f
new file mode 100644
index 00000000..27a874b1
--- /dev/null
+++ b/src/mo_basis/mos_aux.irp.f
@@ -0,0 +1,53 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, mo_coef_aux, (ao_num,mo_num)]
+
+  implicit none
+  integer                        :: i, j
+  logical                        :: exists
+  double precision, allocatable  :: buffer(:,:)
+
+  PROVIDE ezfio_filename
+
+  if (mpi_master) then
+    ! Coefs
+    call ezfio_has_mo_basis_mo_coef_aux(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read mo_coef_aux with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      call ezfio_get_mo_basis_mo_coef_aux(mo_coef_aux)
+      write(*,*) 'Read  mo_coef_aux'
+    endif
+    IRP_IF MPI
+      call MPI_BCAST(mo_coef_aux, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read mo_coef_aux with MPI'
+      endif
+    IRP_ENDIF
+  else
+    ! Orthonormalized AO basis
+    do i = 1, mo_num
+      do j = 1, ao_num
+        mo_coef_aux(j,i) = ao_ortho_canonical_coef(j,i)
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+
diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f
index 0adc49fb..d1a044cb 100644
--- a/src/non_h_ints_mu/jast_deriv.irp.f
+++ b/src/non_h_ints_mu/jast_deriv.irp.f
@@ -408,24 +408,71 @@ subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
   implicit none
   double precision, intent(in)  :: r1(3), r2(3)
   double precision, intent(out) :: mu_val, mu_der(3)
+  double precision              :: r(3)
+  double precision              :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
+  double precision              :: dm_tot, tmp1, tmp2, tmp3
 
   if(j1b_type .eq. 200) then
 
-    double precision :: r(3), dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
-    double precision :: dm_tot, tmp
+    !
+    ! r = 0.5 (r1 + r2)
+    !
+    ! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho)
+    !
+    ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
+    ! d mu[rho(r)] / dx  = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
+    !
+
+    PROVIDE mu_r_ct mu_erf
 
-    PROVIDE mu_r_ct
     r(1) = 0.5d0 * (r1(1) + r2(1))
     r(2) = 0.5d0 * (r1(2) + r2(2))
     r(3) = 0.5d0 * (r1(3) + r2(3))
 
     call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)
-    dm_tot    = dm_a(1) + dm_b(1)
-    mu_val    = mu_r_ct * dsqrt(dm_tot)
-    tmp       = 0.25d0 * mu_r_ct / dm_tot
-    mu_der(1) = tmp * (grad_dm_a(1,1) + grad_dm_b(1,1))
-    mu_der(2) = tmp * (grad_dm_a(2,1) + grad_dm_b(2,1))
-    mu_der(3) = tmp * (grad_dm_a(3,1) + grad_dm_b(3,1))
+
+    dm_tot = dm_a(1) + dm_b(1)
+    tmp1   = dsqrt(dm_tot)
+    tmp2   = mu_erf * dexp(-dm_tot)
+
+    mu_val = mu_r_ct * tmp1 + tmp2
+
+    mu_der = 0.d0
+    if(dm_tot .lt. 1d-7) return
+
+    tmp3      = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2
+    mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
+    mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
+    mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
+
+  elseif(j1b_type .eq. 201) then
+
+    !
+    ! r = 0.5 (r1 + r2)
+    !
+    ! mu[rho(r)] = alpha rho + mu0 exp(-rho)
+    !
+    ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
+    ! d mu[rho(r)] / dx  = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
+    !
+
+    PROVIDE mu_r_ct mu_erf
+
+    r(1) = 0.5d0 * (r1(1) + r2(1))
+    r(2) = 0.5d0 * (r1(2) + r2(2))
+    r(3) = 0.5d0 * (r1(3) + r2(3))
+
+    call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)
+
+    dm_tot = dm_a(1) + dm_b(1)
+    tmp2   = mu_erf * dexp(-dm_tot)
+
+    mu_val = mu_r_ct * dm_tot + tmp2
+
+    tmp3      = 0.5d0 * (mu_r_ct - tmp2)
+    mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
+    mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
+    mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
 
   else
 
diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f
index e5f123a7..980d12de 100644
--- a/src/tc_bi_ortho/print_tc_energy.irp.f
+++ b/src/tc_bi_ortho/print_tc_energy.irp.f
@@ -9,6 +9,10 @@ program print_tc_energy
   my_n_pt_a_grid = 50
   read_wf = .True.
   touch read_wf
+
+  PROVIDE j1b_type
+  print*, 'j1b_type = ', j1b_type
+
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
  call write_tc_energy
 end
diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg
index 2e14488e..fad6f1c2 100644
--- a/src/tc_keywords/EZFIO.cfg
+++ b/src/tc_keywords/EZFIO.cfg
@@ -176,7 +176,7 @@ default: 1.e-7
 type: logical
 doc: If |true|, the integrals of the three-body jastrow are computed with cycles
 interface: ezfio,provider,ocaml
-default: True
+default: Flase
 
 [thresh_biorthog_diag]
 type: Threshold