scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Added PW86 (unrevised, PRB 33, 8800 (1986)) and B86B (JCP 85 (1986) 

7184) exchange functionals. 'pw86pbe' and 'b86bpbe' are shorthand for
pw86 exchange + pbe correlation and b86b exchange + pbe
correlation. 



git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@10454 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
aoterodelaroza 2013-09-17 01:34:16 +00:00
parent db3a37a92a
commit 2a6b38c368
2 changed files with 180 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

114
Modules/funct.f90
View File

@ -89,12 +89,14 @@ module funct
! "blyp" = "sla+b88+lyp+blyp" = BLYP
! "pbe" = "sla+pw+pbx+pbc" = PBE
! "revpbe"= "sla+pw+rpb+pbc" = revPBE (Zhang-Yang)
! "pw86pbe" = "sla+pw+pw86+pbc" = PW86 exchange + PBE correlation
! "b86bpbe" = "sla+pw+b86b+pbc" = B86b exchange + PBE correlation
! "pbesol"= "sla+pw+psx+psc" = PBEsol
! "q2d" = "sla+pw+q2dx+q2dc" = PBEQ2D
! "hcth" = "nox+noc+hcth+hcth" = HCTH/120
! "olyp" = "nox+lyp+optx+blyp" = OLYP
! "wc" = "sla+pw+wcx+pbc" = Wu-Cohen
! "sogga = "sla+pw+sox+pbec" = SOGGA
! "sogga" = "sla+pw+sox+pbec" = SOGGA
! "tpss" = "sla+pw+tpss+tpss" = TPSS Meta-GGA
! "m06l" = "nox+noc+m6lx+m6lc" = M06L Meta-GGA
! "pbe0" = "pb0x+pw+pb0x+pbc" = PBE0
@ -155,6 +157,8 @@ module funct
! gau-pbe in
! "gaup" Gau-PBE hybrid exchange igcx =20
! gau-pbe out
! "pw86" Perdew-Wang (1986) exchange igcx =21
! "b86b" Becke (1986) exchange igcx =22
!
! Gradient Correction on Correlation:
! "nogc" none igcc =0 (default)
@ -265,8 +269,8 @@ module funct
integer :: nxc, ncc, ngcx, ngcc, ncnl
! gau-pbe in
! parameter (nxc = 8, ncc =11, ngcx =19, ngcc = 12, ncnl=2)
parameter (nxc = 8, ncc =11, ngcx =20, ngcc = 12, ncnl=2)
! parameter (nxc = 8, ncc =11, ngcx =21, ngcc = 12, ncnl=2)
parameter (nxc = 8, ncc =11, ngcx =22, ngcc = 12, ncnl=2)
! gau-pbe out
character (len=4) :: exc, corr
@ -279,7 +283,7 @@ module funct
data gradx / 'NOGX', 'B88', 'GGX', 'PBX', 'RPB', 'HCTH', 'OPTX',&
'TPSS', 'PB0X', 'B3LP','PSX', 'WCX', 'HSE', 'RW86', 'PBE', &
'META', 'C09X', 'SOX', 'M6LX', 'Q2DX', 'GAUP' /
'META', 'C09X', 'SOX', 'M6LX', 'Q2DX', 'GAUP', 'PW86', 'B86B' /
data gradc / 'NOGC', 'P86', 'GGC', 'BLYP', 'PBC', 'HCTH', 'TPSS',&
'B3LP', 'PSC', 'PBE', 'META', 'M6LC', 'Q2DC' /
@ -344,6 +348,24 @@ CONTAINS
call set_dft_value (inlc, 0)
dft_defined = .true.
elseif ('PW86PBE' .EQ. TRIM(dftout) ) then
! special case : PW86PBE
call set_dft_value (iexch,1) !Default
call set_dft_value (icorr,4)
call set_dft_value (igcx, 21)
call set_dft_value (igcc, 4)
call set_dft_value (inlc, 0)
dft_defined = .true.
elseif ('B86BPBE' .EQ. TRIM(dftout) ) then
! special case : B86BPBE
call set_dft_value (iexch,1) !Default
call set_dft_value (icorr,4)
call set_dft_value (igcx, 22)
call set_dft_value (igcc, 4)
call set_dft_value (inlc, 0)
dft_defined = .true.
else if ('RPBE' .EQ. TRIM(dftout)) then
! special case : RPBE
call errore('set_dft_from_name', &
@ -1421,6 +1443,10 @@ subroutine gcxc (rho, grho, sx, sc, v1x, v2x, v1c, v2c)
v1x = v1x - exx_fraction * v1xsr
v2x = v2x - exx_fraction * v2xsr
endif
elseif (igcx == 21) then ! 'pw86'
call pw86 (rho, grho, sx, v1x, v2x)
elseif (igcx == 22) then ! 'b86b'
call becke86b (rho, grho, sx, v1x, v2x)
! gau-pbe out
else
sx = 0.0_DP
@ -1641,6 +1667,44 @@ subroutine gcx_spin (rhoup, rhodw, grhoup2, grhodw2, &
v2xup = 2.0_DP * v2xup
v2xdw = 2.0_DP * v2xdw
elseif (igcx == 21) then ! 'PW86'
if (rhoup > small .and. sqrt (abs (grhoup2) ) > small) then
call pw86 (2.0_DP * rhoup, 4.0_DP * grhoup2, sxup, v1xup, v2xup)
else
sxup = 0.0_DP
v1xup = 0.0_DP
v2xup = 0.0_DP
endif
if (rhodw > small .and. sqrt (abs (grhodw2) ) > small) then
call pw86 (2.0_DP * rhodw, 4.0_DP * grhodw2, sxdw, v1xdw, v2xdw)
else
sxdw = 0.0_DP
v1xdw = 0.0_DP
v2xdw = 0.0_DP
endif
sx = 0.5_DP * (sxup + sxdw)
v2xup = 2.0_DP * v2xup
v2xdw = 2.0_DP * v2xdw
elseif (igcx == 22) then ! 'B86B'
if (rhoup > small .and. sqrt (abs (grhoup2) ) > small) then
call becke86b (2.0_DP * rhoup, 4.0_DP * grhoup2, sxup, v1xup, v2xup)
else
sxup = 0.0_DP
v1xup = 0.0_DP
v2xup = 0.0_DP
endif
if (rhodw > small .and. sqrt (abs (grhodw2) ) > small) then
call becke86b (2.0_DP * rhodw, 4.0_DP * grhodw2, sxdw, v1xdw, v2xdw)
else
sxdw = 0.0_DP
v1xdw = 0.0_DP
v2xdw = 0.0_DP
endif
sx = 0.5_DP * (sxup + sxdw)
v2xup = 2.0_DP * v2xup
v2xdw = 2.0_DP * v2xdw
! case igcx == 5 (HCTH) and 6 (OPTX) not implemented
! case igcx == 7 (meta-GGA) must be treated in a separate call to another
! routine: needs kinetic energy density in addition to rho and grad rho
@ -1797,6 +1861,48 @@ subroutine gcx_spin_vec(rhoup, rhodw, grhoup2, grhodw2, &
v2xup = 2.0_DP * v2xup
v2xdw = 2.0_DP * v2xdw
case(21) ! 'pw86'
do i=1,length
if (rhoup(i) > small .and. sqrt(abs(grhoup2(i))) > small) then
call pw86 (2.0_DP * rhoup(i), 4.0_DP * grhoup2(i), sxup(i), v1xup(i), v2xup(i))
else
sxup(i) = 0.0_DP
v1xup(i) = 0.0_DP
v2xup(i) = 0.0_DP
endif
if (rhodw(i) > small .and. sqrt(abs(grhodw2(i))) > small) then
call pw86 (2.0_DP * rhodw(i), 4.0_DP * grhodw2(i), sxdw(i), v1xdw(i), v2xdw(i))
else
sxdw(i) = 0.0_DP
v1xdw(i) = 0.0_DP
v2xdw(i) = 0.0_DP
endif
end do
sx = 0.5_DP * (sxup + sxdw)
v2xup = 2.0_DP * v2xup
v2xdw = 2.0_DP * v2xdw
case(22) ! 'b86b'
do i=1,length
if (rhoup(i) > small .and. sqrt(abs(grhoup2(i))) > small) then
call becke86b (2.0_DP * rhoup(i), 4.0_DP * grhoup2(i), sxup(i), v1xup(i), v2xup(i))
else
sxup(i) = 0.0_DP
v1xup(i) = 0.0_DP
v2xup(i) = 0.0_DP
endif
if (rhodw(i) > small .and. sqrt(abs(grhodw2(i))) > small) then
call becke86b (2.0_DP * rhodw(i), 4.0_DP * grhodw2(i), sxdw(i), v1xdw(i), v2xdw(i))
else
sxdw(i) = 0.0_DP
v1xdw(i) = 0.0_DP
v2xdw(i) = 0.0_DP
endif
end do
sx = 0.5_DP * (sxup + sxdw)
v2xup = 2.0_DP * v2xup
v2xdw = 2.0_DP * v2xdw
case default
call errore ('gcx_spin_vec', 'not implemented', igcx)
end select

70
flib/functionals.f90
View File

@ -416,6 +416,39 @@ subroutine gl (rs, ec, vc)
end subroutine gl
!
!-----------------------------------------------------------------------
subroutine becke86b(rho, grho, sx, v1x, v2x)
!-----------------------------------------------------------------------
! Becke 1986 gradient correction to exchange
! xA.D. Becke, J. Chem. Phys. 85 (1986) 7184
!
USE kinds, ONLY : DP
implicit none
real(DP) :: rho, grho, sx, v1x, v2x
real(DP) :: arho, agrho, beta, gamma
parameter (beta = 0.00375_dp, gamma=0.007_dp)
real(dp) :: sgp1, sgp1_45, sgp1_95
real(dp) :: rdg2_43, rdg2_73, rdg2_83, rdg2_4, rdg4_5
arho = 0.5d0 * rho
agrho = 0.25d0 * grho
rdg2_43 = agrho / arho**(4d0/3d0)
rdg2_73 = rdg2_43 / arho
rdg2_83 = rdg2_43 * rdg2_43 / agrho
rdg2_4 = rdg2_43 * rdg2_83 / agrho
rdg4_5 = rdg2_73 * rdg2_83
sgp1 = 1d0 + gamma * rdg2_83
sgp1_45 = sgp1**(-4d0/5d0)
sgp1_95 = sgp1_45 / sgp1
sx = -2d0 * beta * agrho / arho**(4d0/3d0) * sgp1_45
v1x = -beta * (-4d0/3d0 * rdg2_73 * sgp1_45 + 32d0/15d0 * gamma * rdg4_5 * sgp1_95)
v2x = -beta * (sgp1_45 * rdg2_43 / agrho - 4d0/5d0 * gamma * rdg2_4 * sgp1_95)
end subroutine becke86b
!
!-----------------------------------------------------------------------
subroutine becke88 (rho, grho, sx, v1x, v2x)
!-----------------------------------------------------------------------
! Becke exchange: A.D. Becke, PRA 38, 3098 (1988)
@ -527,6 +560,43 @@ subroutine rPW86 (rho, grho, sx, v1x, v2x)
end subroutine rPW86
subroutine PW86 (rho, grho, sx, v1x, v2x)
!-----------------------------------------------------------------------
! Perdew-Wang 1986 exchange gradient correction: PRB 33, 8800 (1986)
!
USE kinds
implicit none
real(DP), intent(in) :: rho, grho
real(DP), intent(out) :: sx, v1x, v2x
real(DP) :: s, s_2, s_3, s_4, s_5, s_6, fs, grad_rho, df_ds
real(DP) :: a, b, c, s_prefactor, Ax, four_thirds
parameter( a = 1.296_dp, b = 14_dp, c = 0.2_dp, s_prefactor = 6.18733545256027_dp, &
Ax = -0.738558766382022_dp, four_thirds = 4._dp/3._dp)
grad_rho = sqrt(grho)
s = grad_rho/(s_prefactor*rho**(four_thirds))
s_2 = s**2
s_3 = s_2 * s
s_4 = s_2**2
s_5 = s_3 * s_2
s_6 = s_2 * s_4
!! Calculation of energy
fs = (1 + a*s_2 + b*s_4 + c*s_6)**(1.d0/15.d0)
sx = Ax * rho**(four_thirds) * (fs-1d0)
!! Calculation of the potential
df_ds = (1.d0/(15.d0*fs**(14d0)))*(2*a*s + 4*b*s_3 + 6*c*s_5)
v1x = Ax*(four_thirds)*(rho**(1.d0/3.d0)*(fs-1d0) &
-grad_rho/(s_prefactor * rho)*df_ds)
v2x = Ax * df_ds/(s_prefactor*grad_rho)
end subroutine PW86
!
!---------------------------------------------------------------
subroutine c09x (rho, grho, sx, v1x, v2x)