From e061ab18d96f56e09540936aca8fcbc3b0e9d5c5 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Tue, 14 Apr 2020 17:45:01 +0200
Subject: [PATCH] added some functions in on_top_from_ueg.irp.f

---
 src/basis_correction/weak_corr_func.irp.f |  2 +-
 src/dft_utils_func/on_top_from_ueg.irp.f  | 27 +++++++++++++++++++++++
 src/mo_basis/mos_in_r.irp.f               | 20 ++++++++---------
 3 files changed, 38 insertions(+), 11 deletions(-)

diff --git a/src/basis_correction/weak_corr_func.irp.f b/src/basis_correction/weak_corr_func.irp.f
index d3e8bf0a..6af5e49d 100644
--- a/src/basis_correction/weak_corr_func.irp.f
+++ b/src/basis_correction/weak_corr_func.irp.f
@@ -8,8 +8,8 @@
  implicit none
  integer :: ipoint,istate
  double precision :: rho_a, rho_b, ec
- logical :: dospin
  double precision :: wall0,wall1,weight,mu
+ logical :: dospin
  dospin = .true. ! JT dospin have to be set to true for open shell
  print*,'Providing ecmd_lda_mu_of_r ...'
 
diff --git a/src/dft_utils_func/on_top_from_ueg.irp.f b/src/dft_utils_func/on_top_from_ueg.irp.f
index a9cc22c5..70560a7a 100644
--- a/src/dft_utils_func/on_top_from_ueg.irp.f
+++ b/src/dft_utils_func/on_top_from_ueg.irp.f
@@ -91,4 +91,31 @@ double precision function h_func(zeta)
 end
 
 
+!-------------------------------------------------------------------------------------------------------------------------------------------
+  subroutine g0_dg0(rho, rho_a, rho_b, g0, dg0drho)
+
+  implicit none
+  BEGIN_DOC
+  ! Give the on-top pair distribution function g0 and its derivative according to rho dg0drho
+  END_DOC
+
+  double precision, intent (in) :: rho, rho_a, rho_b
+  double precision, intent (out) :: g0, dg0drho
+  double precision :: pi
+  double precision :: g0_UEG_mu_inf, dg0drs
+  double precision :: C1, F1, D1, E1, B1, rs
+
+  pi = dacos(-1.d0)
+  C1 = 0.0819306d0
+  F1 = 0.752411d0
+  D1 = -0.0127713d0
+  E1 = 0.00185898d0
+  B1 = 0.7317d0 - F1
+  rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0) 
+   
+  g0 = g0_UEG_mu_inf(rho_a, rho_b)
+  dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*exp(-F1*rs)
+  dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
+ 
+  end subroutine g0_dg0 
 
diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f
index 7759b222..049db8aa 100644
--- a/src/mo_basis/mos_in_r.irp.f
+++ b/src/mo_basis/mos_in_r.irp.f
@@ -37,18 +37,18 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la
  integer :: i,j,k
  double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3)
  call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
- mos_array=0d0
- mos_grad_array=0d0
- mos_lapl_array=0d0
+ mos_array      = 0.d0
+ mos_grad_array = 0.d0
+ mos_lapl_array = 0.d0
  do j = 1, mo_num
   do k=1, ao_num
-   mos_array(j) += mo_coef(k,j)*aos_array(k)
-   mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1)
-   mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
-   mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
-   mos_lapl_array(j,1) += mo_coef(k,j)*aos_lapl_array(k,1)
-   mos_lapl_array(j,2) += mo_coef(k,j)*aos_lapl_array(k,2)
-   mos_lapl_array(j,3) += mo_coef(k,j)*aos_lapl_array(k,3)
+   mos_array(j)        += mo_coef(k,j) * aos_array(k)
+   mos_grad_array(j,1) += mo_coef(k,j) * aos_grad_array(k,1)
+   mos_grad_array(j,2) += mo_coef(k,j) * aos_grad_array(k,2)
+   mos_grad_array(j,3) += mo_coef(k,j) * aos_grad_array(k,3)
+   mos_lapl_array(j,1) += mo_coef(k,j) * aos_lapl_array(k,1)
+   mos_lapl_array(j,2) += mo_coef(k,j) * aos_lapl_array(k,2)
+   mos_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3)
   enddo
  enddo
 end