src/dft_utils_one_e/ec_scan.irp.f

@@ -37,8 +37,7 @@ double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2)
  gama     = 0.031091d0
  ! correlation energy lsda1
  call ec_only_lda_sr(0.d0,nup,ndo,e_c_lsda1)
- ! correlation energy per particle 
- e_c_lsda1 = e_c_lsda1/rho         
+          
  xi       = spin_d/rho
  rs       = (cst_43 * pi * rho)**(-cst_13)
  s        = drho/( 2.d0 * cst_3pi2**(cst_13) * rho**cst_43  )
@@ -62,11 +61,7 @@ double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2)
  g_at2    = 1.d0/(1.d0 + 4.d0 * a*t*t)**0.25d0
  h1       = gama * phi_3 * dlog(1.d0 + w_1 * (1.d0 - g_at2))
  ! interpolation function 
- if(cst_1alph.gt.0.d0)then
-  fc_alpha = dexp(-c_1c * alpha * inv_1alph) 
- else
-  fc_alpha = - d_c * dexp(c_2c * inv_1alph) 
- endif
+ fc_alpha = dexp(-c_1c * alpha * inv_1alph) * step_f(cst_1alph) - d_c * dexp(c_2c * inv_1alph) * step_f(-cst_1alph) 
  ! first part of the correlation energy 
  e_c_1    = e_c_lsda1 + h1
  
@@ -87,10 +82,19 @@ double precision function ec_scan(rho_a,rho_b,tau,grad_rho_2)
  ec_scan = e_c_1 + fc_alpha * (e_c_0 - e_c_1)
 end
 
+double precision function step_f(x)
+ implicit none
+ double precision, intent(in) :: x
+ if(x.lt.0.d0)then
+  step_f = 0.d0
+ else
+  step_f = 1.d0 
+ endif
+end             
+
 double precision function beta_rs(rs)
  implicit none
  double precision, intent(in) ::rs
- beta_rs = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs)
+ beta_rs(rs) = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs)
 
 end
-