73 lines
3.0 KiB
Fortran
73 lines
3.0 KiB
Fortran
|
|
|
|
subroutine ec_md_on_top_PBE_mu_corrected(mu,r,two_dm,eps_c_md_on_top_PBE)
|
|
implicit none
|
|
BEGIN_DOC
|
|
! enter with "r(3)", and "two_dm(N_states)" which is the on-top pair density at "r" for each states
|
|
!
|
|
! you get out with the energy density defined in J. Chem. Phys.150, 084103 (2019); doi: 10.1063/1.508263
|
|
!
|
|
! by Eq. (26), which includes the correction of the on-top pair density of Eq. (29).
|
|
END_DOC
|
|
double precision, intent(in) :: mu , r(3), two_dm
|
|
double precision, intent(out) :: eps_c_md_on_top_PBE(N_states)
|
|
double precision :: two_dm_in_r, pi, e_pbe(N_states),beta(N_states),mu_correction_of_on_top
|
|
double precision :: aos_array(ao_num), grad_aos_array(3,ao_num)
|
|
double precision :: rho_a(N_states),rho_b(N_states)
|
|
double precision :: grad_rho_a(3,N_states),grad_rho_b(3,N_states)
|
|
double precision :: grad_rho_a_2(N_states),grad_rho_b_2(N_states),grad_rho_a_b(N_states)
|
|
double precision :: rhoc,rhoo,sigmacc,sigmaco,sigmaoo,vrhoc,vrhoo,vsigmacc,vsigmaco,vsigmaoo,on_top_corrected
|
|
integer :: m, istate
|
|
|
|
pi = 4.d0 * datan(1.d0)
|
|
|
|
eps_c_md_on_top_PBE = 0.d0
|
|
call density_and_grad_alpha_beta_and_all_aos_and_grad_aos_at_r(r,rho_a,rho_b, grad_rho_a, grad_rho_b, aos_array, grad_aos_array)
|
|
grad_rho_a_2 = 0.d0
|
|
grad_rho_b_2 = 0.d0
|
|
grad_rho_a_b = 0.d0
|
|
do istate = 1, N_states
|
|
do m = 1, 3
|
|
grad_rho_a_2(istate) += grad_rho_a(m,istate)*grad_rho_a(m,istate)
|
|
grad_rho_b_2(istate) += grad_rho_b(m,istate)*grad_rho_b(m,istate)
|
|
grad_rho_a_b(istate) += grad_rho_a(m,istate)*grad_rho_b(m,istate)
|
|
enddo
|
|
enddo
|
|
do istate = 1, N_states
|
|
! convertion from (alpha,beta) formalism to (closed, open) formalism
|
|
call rho_ab_to_rho_oc(rho_a(istate),rho_b(istate),rhoo,rhoc)
|
|
call grad_rho_ab_to_grad_rho_oc(grad_rho_a_2(istate),grad_rho_b_2(istate),grad_rho_a_b(istate),sigmaoo,sigmacc,sigmaco)
|
|
|
|
! usual PBE correlation energy using the density, spin polarization and density gradients for alpha/beta electrons
|
|
call ec_pbe_only(0.d0,rhoc,rhoo,sigmacc,sigmaco,sigmaoo,e_PBE(istate))
|
|
|
|
! correction of the on-top pair density according to Eq. (29)
|
|
on_top_corrected = mu_correction_of_on_top(mu,two_dm)
|
|
|
|
! quantity of Eq. (27) with a factor two according to the difference of normalization
|
|
! between the on-top of the JCP paper and that of QP2
|
|
beta(istate) = (3.d0*e_PBE(istate))/( (-2.d0+sqrt(2d0))*sqrt(2.d0*pi)*2.d0* on_top_corrected)
|
|
|
|
! quantity of Eq. (26)
|
|
eps_c_md_on_top_PBE(istate)=e_PBE(istate)/(1.d0+beta(istate)*mu**3)
|
|
enddo
|
|
end
|
|
|
|
|
|
double precision function mu_correction_of_on_top(mu,on_top)
|
|
implicit none
|
|
BEGIN_DOC
|
|
! mu-based correction to the on-top pair density provided by the assymptotic expansion of
|
|
!
|
|
! P. Gori-Giorgi and A. Savin, Phys. Rev. A73, 032506 (2006)
|
|
!
|
|
! This is used in J. Chem. Phys.150, 084103 (2019); Eq. (29).
|
|
END_DOC
|
|
double precision, intent(in) :: mu,on_top
|
|
double precision :: pi
|
|
pi = 4.d0 * datan(1.d0)
|
|
mu_correction_of_on_top = on_top / ( 1.d0 + 2.d0/(dsqrt(pi)*mu) )
|
|
mu_correction_of_on_top = max(mu_correction_of_on_top ,1.d-15)
|
|
end
|
|
|