diff --git a/src/basis_correction/pbe_ueg_self_contained.irp.f b/src/basis_correction/pbe_ueg_self_contained.irp.f
new file mode 100644
index 00000000..afedfc9c
--- /dev/null
+++ b/src/basis_correction/pbe_ueg_self_contained.irp.f
@@ -0,0 +1,62 @@
+double precision function ecmd_pbe_ueg_self_cont(dens,spin_pol,mu,e_PBE)
+ implicit none
+ ! dens = total density 
+ ! spin_pol = spin_polarization (n_a - n_b)/dens
+ ! e_PBE = PBE correlation (mu=0) energy evaluated at (dens,spin_pol,grad_rho) 
+ ! e_PBE = epsilon_PBE * dens
+ ! dens = a + b
+ ! spin_pol = (a - b)/(a+b)
+ ! spin_pol * dens = a - b
+ ! a - b + a+b = 2 a
+ ! a - b - a - b = - 2b
+ double precision, intent(in) :: dens,spin_pol,mu,e_PBE
+ double precision :: rho_a,rho_b,pi,g0_UEG_func,denom,beta
+ pi = dacos(-1.d0)
+ rho_a = (dens * spin_pol + dens)*0.5d0
+ rho_b = (dens - dens * spin_pol)*0.5d0 
+ if(mu == 0.d0) then
+  ecmd_pbe_ueg_self_cont = e_PBE
+ else
+! note: the on-top pair density is (1-zeta^2) rhoc^2 g0 = 4 rhoa * rhob * g0
+  denom = (-2.d0+sqrt(2d0))*sqrt(2.d0*pi) * 4.d0*rho_a*rho_b*g0_UEG_func(rho_a,rho_b)
+  if (dabs(denom) > 1.d-12) then
+   beta = (3.d0*e_PBE)/denom
+   ecmd_pbe_ueg_self_cont=e_PBE/(1.d0+beta*mu**3)
+  else
+   ecmd_pbe_ueg_self_cont=0.d0
+  endif
+ endif
+end
+
+double precision function g0_UEG_func(rho_a,rho_b)
+! Pair distribution function g0(n_alpha,n_beta) of the Colombic UEG 
+!
+! Taken from Eq. (46)  P. Gori-Giorgi and A. Savin, Phys. Rev. A 73, 032506 (2006).                                                                                                    
+ implicit none
+ double precision, intent(in) :: rho_a,rho_b
+ double precision :: rho,pi,x
+ double precision :: B, C, D, E, d2, rs, ahd
+ rho = rho_a+rho_b
+ pi = 4d0 * datan(1d0)
+ ahd = -0.36583d0
+ d2 = 0.7524d0
+ B = -2d0 * ahd - d2
+ C = 0.08193d0
+ D = -0.01277d0
+ E = 0.001859d0
+ x = -d2*rs
+ if (dabs(rho) > 1.d-20) 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
+   g0_UEG_func= 0.5d0 * (1d0+ rs* (-B + rs*(C + rs*(D + rs*E))))*dexp(x)
+  else
+   g0_UEG_func= 0.d0
+  endif
+ else
+  g0_UEG_func= 0.d0
+ endif
+ g0_UEG_func = max(g0_UEG_func,1.d-14)
+
+end
+
diff --git a/src/basis_correction/print_su_pbe_ot.irp.f b/src/basis_correction/print_su_pbe_ot.irp.f
index 49f90ade..dcd4b6f1 100644
--- a/src/basis_correction/print_su_pbe_ot.irp.f
+++ b/src/basis_correction/print_su_pbe_ot.irp.f
@@ -20,6 +20,7 @@ subroutine print_su_pbe_ot
  integer :: istate 
  do istate = 1, N_states
   write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD PBE-UEG       , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
+  write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ecmd_pbe_ueg_test  , state ',istate,' = ',ecmd_pbe_ueg_test(istate)
  enddo
  do istate = 1, N_states
   write(*, '(A29,X,I3,X,A3,X,F16.10)') '  ECMD SU-PBE-OT     , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)
diff --git a/src/basis_correction/weak_corr_func.irp.f b/src/basis_correction/weak_corr_func.irp.f
index 6af5e49d..4d4f7075 100644
--- a/src/basis_correction/weak_corr_func.irp.f
+++ b/src/basis_correction/weak_corr_func.irp.f
@@ -81,3 +81,54 @@ BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_mu_of_r, (N_states)]
  print*,'Time for the ecmd_pbe_ueg_mu_of_r:',wall1-wall0
 
 END_PROVIDER
+
+BEGIN_PROVIDER [double precision, ecmd_pbe_ueg_test, (N_states)]
+BEGIN_DOC
+! test of the routines contained in pbe_ueg_self_contained.irp.f
+END_DOC
+implicit none
+double precision :: weight
+integer :: ipoint,istate,m
+double precision :: mu,rho_a,rho_b
+double precision :: dens,spin_pol,grad_rho,e_PBE,delta_rho
+double precision :: ecmd_pbe_ueg_self_cont,eps_c_md_PBE
+ecmd_pbe_ueg_test = 0.d0
+ 
+do istate = 1, N_states
+ do ipoint = 1, n_points_final_grid
+  weight=final_weight_at_r_vector(ipoint)
+
+  ! mu(r) defined by Eq. (37) of J. Chem. Phys. 149, 194301 (2018)
+  mu    = mu_of_r_prov(ipoint,istate)
+
+  ! conversion from rho_a,rho_b --> dens,spin_pol
+  rho_a = one_e_dm_and_grad_alpha_in_r(4,ipoint,istate)
+  rho_b = one_e_dm_and_grad_beta_in_r(4,ipoint,istate)
+  dens = rho_a + rho_b
+  spin_pol = (rho_a - rho_b)/(max(dens,1.d-12))
+  delta_rho = rho_a - rho_b
+
+  ! conversion from grad_rho_a ... to sigma 
+  double precision :: grad_rho_a(3),grad_rho_b(3),grad_rho_a_2(3),grad_rho_b_2(3),grad_rho_a_b(3)
+  double precision :: sigmacc,sigmaco,sigmaoo
+  grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
+  grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
+  grad_rho_a_2 = 0.d0
+  grad_rho_b_2 = 0.d0
+  grad_rho_a_b = 0.d0
+  do m = 1, 3
+   grad_rho_a_2 += grad_rho_a(m)*grad_rho_a(m)
+   grad_rho_b_2 += grad_rho_b(m)*grad_rho_b(m)
+   grad_rho_a_b += grad_rho_a(m)*grad_rho_b(m)
+  enddo
+  call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2,grad_rho_b_2,grad_rho_a_b,sigmaoo,sigmacc,sigmaco)
+
+  ! call the PBE energy 
+  call ec_pbe_only(0.d0,dens,delta_rho,sigmacc,sigmaco,sigmaoo,e_PBE)
+  eps_c_md_PBE  = ecmd_pbe_ueg_self_cont(dens,spin_pol,mu,e_PBE)
+
+  ecmd_pbe_ueg_test(istate) += eps_c_md_PBE * weight
+ enddo
+enddo
+!
+END_PROVIDER
diff --git a/src/dft_utils_func/ecmd_pbe_general.irp.f b/src/dft_utils_func/ecmd_pbe_general.irp.f
index cf58092c..48f1608d 100644
--- a/src/dft_utils_func/ecmd_pbe_general.irp.f
+++ b/src/dft_utils_func/ecmd_pbe_general.irp.f
@@ -57,4 +57,4 @@ subroutine ec_md_pbe_on_top_general(mu,rho_a,rho_b,grad_rho_a,grad_rho_b,on_top,
   endif
  end
 
- 
+