added test_ueg_self_contained.irp.f

src/basis_correction/pbe_ueg_self_contained.irp.f

@@ -2,13 +2,8 @@ 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
+ ! e_PBE = PBE correlation (mu=0) energy evaluated at  dens,spin_pol (and grad_rho) 
+ ! e_PBE = epsilon_PBE * dens which means that it is not the energy density but the energy density X the density
  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)
@@ -31,7 +26,7 @@ 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).                                                                                                    
+! 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
@@ -46,7 +41,7 @@ double precision function g0_UEG_func(rho_a,rho_b)
  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
+  rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) 
   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)

src/basis_correction/print_su_pbe_ot.irp.f

@@ -22,8 +22,8 @@ subroutine print_su_pbe_ot
   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)
- 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)
+! enddo
 
 end

src/basis_correction/test_ueg_self_contained.irp.f

@@ -0,0 +1,84 @@
+program test_sc
+ implicit none
+ integer :: m
+ double precision :: r(3),f_hf,on_top,mu,sqpi
+ double precision :: rho_a,rho_b,w_hf,dens,delta_rho,e_pbe
+ 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,spin_pol
+ double precision :: eps_c_md_PBE , ecmd_pbe_ueg_self_cont
+ r = 0.D0
+ r(3) = 1.D0
+ call f_HF_valence_ab(r,r,f_hf,on_top)
+ sqpi = dsqrt(dacos(-1.d0))
+ if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then
+  w_hf   = 1.d+10
+ else
+  w_hf  = f_hf /  on_top
+ endif                            
+ mu = sqpi * 0.5d0 * w_hf
+ call density_and_grad_alpha_beta(r,rho_a,rho_b, grad_rho_a, grad_rho_b)                                                                              
+ dens = rho_a + rho_b
+ delta_rho = rho_a - rho_b
+ spin_pol = delta_rho/(max(1.d-10,dens))
+ 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 
+ print*,'f_hf,on_top = ',f_hf,on_top
+ print*,'mu = ',mu
+ print*,'dens,spin_pol',dens,spin_pol
+ call ec_pbe_only(0.d0,dens,delta_rho,sigmacc,sigmaco,sigmaoo,e_PBE)
+ print*,'e_PBE = ',e_PBE
+ eps_c_md_PBE  = ecmd_pbe_ueg_self_cont(dens,spin_pol,mu,e_PBE)
+ print*,'eps_c_md_PBE = ',eps_c_md_PBE
+
+ print*,''
+ print*,''
+ print*,''
+ print*,'energy_c' ,energy_c
+
+ integer::ipoint
+ double precision :: weight , accu
+ accu = 0.d0
+ do ipoint = 1, n_points_final_grid
+  r =  final_grid_points(:,ipoint)
+  weight = final_weight_at_r_vector(ipoint)
+  call f_HF_valence_ab(r,r,f_hf,on_top)
+  sqpi = dsqrt(dacos(-1.d0))
+  if(on_top.le.1.d-12.or.f_hf.le.0.d0.or.f_hf * on_top.lt.0.d0)then
+   w_hf   = 1.d+10
+  else
+   w_hf  = f_hf /  on_top
+  endif
+  mu = sqpi * 0.5d0 * w_hf
+  call density_and_grad_alpha_beta(r,rho_a,rho_b, grad_rho_a, grad_rho_b)
+  dens = rho_a + rho_b
+  delta_rho = rho_a - rho_b
+  spin_pol = delta_rho/(max(1.d-10,dens))
+  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)
+  write(33,'(100(F16.10,X))')r(:), weight, w_hf, on_top, mu, dens, spin_pol, e_PBE, eps_c_md_PBE
+  accu += weight * eps_c_md_PBE
+ enddo
+ print*,'accu = ',accu
+  write(*, *) '  ECMD PBE-UEG       ',ecmd_pbe_ueg_mu_of_r(1)
+  write(*, *) '  ecmd_pbe_ueg_test  ',ecmd_pbe_ueg_test(1)
+ 
+end

src/becke_numerical_grid/grid_becke.irp.f

@@ -64,7 +64,8 @@ END_PROVIDER
 
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular,n_points_radial_grid,nucl_num)]
+ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_angular,n_points_radial_grid,nucl_num)]
+&BEGIN_PROVIDER [double precision, radial_points_per_atom, (n_points_radial_grid,nucl_num)]
   BEGIN_DOC
   ! x,y,z coordinates of grid points used for integration in 3d space
   END_DOC
@@ -72,6 +73,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
   integer                        :: i,j,k
   double precision               :: dr,x_ref,y_ref,z_ref
   double precision               :: knowles_function
+  radial_points_per_atom = 0.D0
   do i = 1, nucl_num
     x_ref = nucl_coord(i,1)
     y_ref = nucl_coord(i,2)
@@ -83,7 +85,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
 
       ! value of the radial coordinate for the integration
       r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x)
-
+      radial_points_per_atom(j,i) = r
       ! explicit values of the grid points centered around each atom
       do k = 1, n_points_integration_angular
         grid_points_per_atom(1,k,j,i) =                              &