diff --git a/src/determinants/generate_cas_space.irp.f b/src/determinants/generate_cas_space.irp.f
index 47a2ca30..05201c74 100644
--- a/src/determinants/generate_cas_space.irp.f
+++ b/src/determinants/generate_cas_space.irp.f
@@ -33,7 +33,7 @@ subroutine generate_cas_space
   print *, 'CAS(', n_alpha_act+n_beta_act, ', ', n_act_orb, ')'
   print *, ''
 
-  n_det_alpha_unique = binom_int(n_act_orb, n_alpha_act)
+  n_det_alpha_unique = int(binom_int(n_act_orb, n_alpha_act),4)
   TOUCH n_det_alpha_unique
 
   n = n_alpha_act
@@ -56,7 +56,7 @@ subroutine generate_cas_space
     u = ior(t1,t2)
   enddo
 
-  n_det_beta_unique = binom_int(n_act_orb, n_beta_act)
+  n_det_beta_unique = int(binom_int(n_act_orb, n_beta_act),4)
   TOUCH n_det_beta_unique
 
   n = n_beta_act
diff --git a/src/mu_of_r/mu_of_r.irp.f b/src/mu_of_r/mu_of_r.irp.f
new file mode 100644
index 00000000..e9cf6c4a
--- /dev/null
+++ b/src/mu_of_r/mu_of_r.irp.f
@@ -0,0 +1,43 @@
+
+subroutine grad_mu_of_r_mean_field(r,mu_mf, dm, grad_mu_mf, grad_dm)
+ implicit none
+ BEGIN_DOC
+ ! returns the value and gradients of the mu(r) mean field, together with the HF density and its gradients. 
+ END_DOC
+ 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
+ 
+ double precision :: dist
+ 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)
+ 
+ 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
+  if(mu_of_r_tc=="Erfmu")then
+   mu_mf           = 0.3333333333d0 * sqpi * (f_mf_ab/two_bod_dens + 0.25d0)
+   grad_mu_mf(1:3) = 0.3333333333d0 * 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)
+  else if(mu_of_r_tc=="Standard")then
+   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)
+  else if(mu_of_r_tc=="Erfmugauss")then
+   mu_mf = (f_mf_ab/two_bod_dens + 0.25d0)/c_mu_gauss_tot 
+   grad_mu_mf(1:3) = 1.d0/c_mu_gauss_tot* (grad_f_mf_ab(1:3) * two_bod_dens - f_mf_ab * grad_two_bod_dens(1:3))& 
+                                    /(two_bod_dens*two_bod_dens)
+  else 
+   print*,'Wrong value for mu_of_r_tc !'
+   stop
+  endif
+ endif 
+
+end
+