added mu of r grad_rho

src/dft_one_e/mu_erf_dft.irp.f

@@ -44,9 +44,21 @@ BEGIN_PROVIDER [double precision, mu_rsc_of_r, (n_points_final_grid)]
   rho = dm_a + dm_b
   mu_rsc_of_r(i) = mu_rs_c(rho)
  enddo
+END_PROVIDER
 
+BEGIN_PROVIDER [double precision, mu_grad_rho, (n_points_final_grid)]
+ implicit none
+ integer :: i
+ double precision :: mu_grad_rho_func, r(3)
+ do i = 1, n_points_final_grid
+  r(1) = final_grid_points(1,i)
+  r(2) = final_grid_points(2,i)
+  r(3) = final_grid_points(3,i)
+  mu_grad_rho(i) = mu_grad_rho_func(r)
+ enddo
 END_PROVIDER 
 
+
 BEGIN_PROVIDER [double precision, mu_of_r_dft_average]
  implicit none
  integer :: i

src/dft_utils_func/mu_of_r_dft.irp.f

@@ -0,0 +1,37 @@
+double precision function mu_rs_c(rho)
+ implicit none
+ double precision, intent(in) :: rho
+ include 'constants.include.F'
+ double precision :: cst_rs,alpha_rs,rs
+ cst_rs   = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0)
+ alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi
+  
+ rs = cst_rs * rho**(-1.d0/3.d0)
+ mu_rs_c =  alpha_rs/dsqrt(rs)
+
+end
+
+double precision function mu_grad_rho_func(r)
+ implicit none
+ double precision , intent(in) :: r(3)
+ integer :: m
+ double precision :: rho, dm_a, dm_b, grad_dm_a(3), grad_dm_b(3)
+ double precision :: eta, grad_rho(3), grad_sqr
+ eta = 0.135d0
+  call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)  
+  rho = dm_a + dm_b
+  do m = 1,3
+   grad_rho(m) = grad_dm_a(m) + grad_dm_b(m)
+  enddo
+  grad_sqr=0.d0
+  do m = 1,3
+   grad_sqr=grad_sqr+grad_rho(m)*grad_rho(m)
+  enddo
+  grad_sqr = dsqrt(grad_sqr)
+  if (rho<1.d-12) then
+   mu_grad_rho_func = 1.d-10
+  else
+   mu_grad_rho_func = eta * grad_sqr / rho
+  endif
+  
+end

src/dft_utils_func/mu_rsc.irp.f

@@ -1,13 +0,0 @@
-double precision function mu_rs_c(rho)
- implicit none
- double precision, intent(in) :: rho
- include 'constants.include.F'
- double precision :: cst_rs,alpha_rs,rs
- cst_rs   = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0)
- alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi
-  
- rs = cst_rs * rho**(-1.d0/3.d0)
- mu_rs_c =  alpha_rs/dsqrt(rs)
-
-end
-