From 2125cd69aba5f560cdc4e697cc2b871f8f7fe0ad Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Thu, 7 Oct 2021 23:05:43 +0200
Subject: [PATCH] added mu(rsc)

---
 src/dft_one_e/mu_erf_dft.irp.f  | 57 +++++++++++++++++++++++++++++----
 src/dft_utils_func/mu_rsc.irp.f | 13 ++++++++
 2 files changed, 63 insertions(+), 7 deletions(-)
 create mode 100644 src/dft_utils_func/mu_rsc.irp.f

diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f
index 8161324b..a1d7a5c1 100644
--- a/src/dft_one_e/mu_erf_dft.irp.f
+++ b/src/dft_one_e/mu_erf_dft.irp.f
@@ -12,11 +12,54 @@ END_PROVIDER
 BEGIN_PROVIDER [double precision, mu_of_r_dft, (n_points_final_grid)]
  implicit none
  integer :: i
- do i = 1, n_points_final_grid
-  if(mu_dft_type == "cst")then
-   mu_of_r_dft(i) = mu_erf_dft
-  else
-   mu_of_r_dft(i) = mu_of_r_hf(i)
-  endif
- enddo
+ if(mu_dft_type == "Read")then
+   call ezfio_get_mu_of_r_mu_of_r_disk(mu_of_r_dft)
+ else
+  do i = 1, n_points_final_grid
+   if(mu_dft_type == "cst")then
+    mu_of_r_dft(i) = mu_erf_dft
+   else if(mu_dft_type == "hf")then
+    mu_of_r_dft(i) = mu_of_r_hf(i)
+   else if(mu_dft_type == "rsc")then
+    mu_of_r_dft(i) = mu_rsc_of_r(i)
+   else 
+    print*,'mu_dft_type is not of good type = ',mu_dft_type
+    print*,'it must be of type Read, cst, hf, rsc'
+    print*,'Stopping ...'
+    stop
+   endif
+  enddo
+ endif
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, mu_rsc_of_r, (n_points_final_grid)]
+ implicit none
+ integer :: i
+ double precision :: mu_rs_c,rho,r(3), dm_a, dm_b
+ 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)
+  call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
+  rho = dm_a + dm_b
+  mu_rsc_of_r(i) = mu_rs_c(rho)
+ enddo
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [double precision, mu_of_r_dft_average]
+ implicit none
+ integer :: i
+ double precision :: mu_rs_c,rho,r(3), dm_a, dm_b
+ mu_of_r_dft_average = 0.d0
+ 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)
+  call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
+  rho = dm_a + dm_b
+  mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i)
+ enddo
+ mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num)
+ print*,'mu_of_r_dft_average = ',mu_of_r_dft_average
 END_PROVIDER 
diff --git a/src/dft_utils_func/mu_rsc.irp.f b/src/dft_utils_func/mu_rsc.irp.f
new file mode 100644
index 00000000..386a1c2d
--- /dev/null
+++ b/src/dft_utils_func/mu_rsc.irp.f
@@ -0,0 +1,13 @@
+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
+