diff --git a/src/basis_correction/pbe_ueg_self_contained.irp.f b/src/basis_correction/pbe_ueg_self_contained.irp.f index afedfc9c..1dbc5e17 100644 --- a/src/basis_correction/pbe_ueg_self_contained.irp.f +++ b/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) diff --git a/src/basis_correction/print_su_pbe_ot.irp.f b/src/basis_correction/print_su_pbe_ot.irp.f index dcd4b6f1..1d9416e5 100644 --- a/src/basis_correction/print_su_pbe_ot.irp.f +++ b/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 diff --git a/src/basis_correction/test_ueg_self_contained.irp.f b/src/basis_correction/test_ueg_self_contained.irp.f new file mode 100644 index 00000000..7f08b875 --- /dev/null +++ b/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 diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f index 79f15c9a..48fd041a 100644 --- a/src/becke_numerical_grid/grid_becke.irp.f +++ b/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) = &