fixed the ao effective potential in DFT

src/dft_one_e/effective_pot.irp.f

@@ -54,31 +54,22 @@ END_PROVIDER
 &BEGIN_PROVIDER [double precision, ao_effective_one_e_potential_without_kin, (ao_num, ao_num,N_states)]
  implicit none
  integer :: i,j,istate
- effective_one_e_potential = 0.d0
+ ao_effective_one_e_potential = 0.d0
+ ao_effective_one_e_potential_without_kin = 0.d0
  BEGIN_DOC
-! Effective_one_e_potential(i,j) = $\rangle i_{MO}| v_{H}^{sr} |j_{MO}\rangle  + \rangle i_{MO}| h_{core} |j_{MO}\rangle  + \rangle i_{MO}|v_{xc} |j_{MO}\rangle$
+! Effective_one_e_potential(i,j) = $\rangle i_{AO}| v_{H}^{sr} |j_{AO}\rangle  + \rangle i_{AO}| h_{core} |j_{AO}\rangle  + \rangle i_{AO}|v_{xc} |j_{AO}\rangle$
 !
 ! on the |MO| basis
 ! 
 ! Taking the expectation value does not provide any energy, but
 !
-! effective_one_e_potential(i,j) is the potential coupling DFT and WFT parts 
+! ao_effective_one_e_potential(i,j) is the potential coupling DFT and WFT parts 
 !
 ! and it is used in any RS-DFT based calculations  
  END_DOC
  do istate = 1, N_states
-  do j = 1, mo_num
-   do i = 1, mo_num
-
-    effective_one_e_potential(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_integrals_n_e(i,j) + mo_kinetic_integrals(i,j)   &
-                                   + 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate)                          &
-                                   +          potential_x_beta_mo(i,j,istate)  + potential_c_beta_mo(i,j,istate)   )
-
-    effective_one_e_potential_without_kin(i,j,istate) = short_range_Hartree_operator(i,j,istate) + mo_integrals_n_e(i,j)                   &
-                                   + 0.5d0 * (potential_x_alpha_mo(i,j,istate) + potential_c_alpha_mo(i,j,istate)                          &
-                                   +          potential_x_beta_mo(i,j,istate)  + potential_c_beta_mo(i,j,istate)   )
-   enddo
-  enddo
+    call mo_to_ao(effective_one_e_potential(1,1,istate),mo_num,ao_effective_one_e_potential(1,1,istate),ao_num)
+    call mo_to_ao(effective_one_e_potential_without_kin(1,1,istate),mo_num,ao_effective_one_e_potential_without_kin(1,1,istate),ao_num)
  enddo
 END_PROVIDER
 

src/dft_utils_func/on_top_from_ueg.irp.f

@@ -32,7 +32,8 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
  C = 0.08193d0
  D = -0.01277d0
  E = 0.001859d0                     
- if (dabs(rho) > 1.d-12) then
+ x = -d2*rs
+ if (dabs(rho) > 1.d-12.and.dabs(x).lt.20.d0) then
   rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
   x = -d2*rs
   if(dabs(x).lt.50.d0)then