subroutine xcssbD (calcv, nspin, nx, n, rho, drho, d2rho, rhotrd, nilrho, e, v, nsplit) use Vartypes implicit NONE integer(KINT), intent(IN) :: nspin, nx, n, nsplit logical, intent(IN) :: calcv, nilrho(n,nspin) real(KREAL), intent(IN) :: rho(n,nspin), rhotrd(n,nspin) real(KREAL), intent(IN) :: drho(n,3,nspin), d2rho(n,6,nspin) real(KREAL), intent(OUT) :: e(n,nsplit), v(nx,nspin) ! ====================================================================== ! purpose: evaluate the energy density and (optionally) the potential ! according to SSB-D's gradient correction formula. ! m.swart ! ! input : calcv - (logical): calculate the potential (if false: only ! the energy density) ! nspin - nr. of independent spins ! nx - size of array (max nr. of points) ! n - actual nr. of points for which the functions have ! to be evaluated. ! rho - charge density (for all spins) ! drho - partial (cartesian) derivatives of rho. ! d2rho - partial second derivatives ! rhotrd - rho**1/3 (for all spins) ! nilrho - (logicals) function result must be set to zero ! (the density is in fact zero, irrespective of the ! value contained in array rho) ! nsplit - 1: normal, 2: energy per spin ! ! in-out : e,v - energy and potential (the latter spin-dependent) ! the values computed in this routine are added to ! the input values. ! ! scratch: grad,rho43,rhom4,x2,tsuml,tlapl,dp1,dp2,dp3,do1,do2,do3, ! dt1,dt2,grd2,vkt,v1,v2,v2,v4 - (n) ! ! remark * if nspin=1, then rho should be half the total charge density, ! i.e. it always expects the spin-density ! ! called from: xc ! ! author: M. Swart, June 2009 ! ====================================================================== real(KREAL), parameter :: ZERO = 0.0_KREAL real(KREAL), parameter :: ONE = 1.0_KREAL real(KREAL), parameter :: TWO = 2.0_KREAL real(KREAL), parameter :: THREE = 3.0_KREAL real(KREAL), parameter :: FOUR = 4.0_KREAL real(KREAL), parameter :: EIGHT = 8.0_KREAL real(KREAL), parameter :: R1O3 = ONE / THREE real(KREAL), parameter :: R4O3 = FOUR / THREE ! Funr = 2**(1/3) real(KREAL), parameter :: Funr = 1.259921049894873_KREAL ! Ax = -(3/4)*((3/pi)**(1/3)) real(KREAL), parameter :: Ax = -0.738558766382022405884230032680836_KREAL real(KREAL), parameter :: Cx = Funr * Ax real(KREAL), parameter :: UM = 0.2195149727645171_KREAL real(KREAL), parameter :: UK = 0.8040_KREAL real(KREAL), parameter :: UL = UM/UK real(KREAL), parameter :: DELTA = 0.1_KREAL real(KREAL), parameter :: XTINY = 1.0e-20_KREAL real(KREAL), allocatable :: grad(:), rho43(:), rhom4(:), x2(:), tsuml(:), tlapl(:), & dp1(:), dp2(:), dp3(:), do1(:), do2(:), do3(:), & dt1(:), dt2(:), grd2(:), v1(:), v2(:), v3(:), v4(:), vkt(:) integer(KINT) :: is, k, istat, isplit real(KREAL) :: factor, rA, rB, rC, rD, rE, rG, pi, ckf, rk, tmpX real(KREAL) :: fA, fB, fC, fD, fE, fF, fU, gamma real(KREAL) :: fV0, fV1, fV2, fV3, fV4 pi = acos(-ONE) ckf = (THREE*pi*pi)**R1O3 rk = TWO*ckf ! ----------------- ! SSB-Dx functional ! ----------------- fA = 1.079966_KREAL fB = 0.197465_KREAL fC = 0.272729_KREAL fD = 0.197465_KREAL fE = 5.873645_KREAL ! -------------------------------------------------- ! parameters for inclusion of KT gradient-functional ! -------------------------------------------------- fU = -0.749940_KREAL fF = 0.949488_KREAL fD = fD * (ONE - fU) ! -------------------------- ! gamma = -fU * fB * fF * Cx ! -------------------------- gamma = -0.130838_KREAL ! -------------------------------------------------------------------------- ! the factor rk comes from the difference between s and x ! the factor Funr comes from how the exchange is calculated (spin-polarized) ! -------------------------------------------------------------------------- rA = fA rB = fB / ((Funr*rk)**TWO) rC = fC / ((Funr*rk)**TWO) rD = fD / ((Funr*rk)**TWO) rE = fE / ((Funr*rk)**FOUR) rG = gamma / ((Funr*rk)**TWO) fV0 = R4O3*Cx fV1 = R4O3*Cx fV2 = R4O3*Cx fV3 = -TWO*Cx fV4 = EIGHT*Cx allocate (grad(n), rho43(n), rhom4(n), x2(n), dp1(n), dp2(n), dp3(n), do1(n), do2(n), do3(n), & tsuml(n), tlapl(n), dt1(n), dt2(n), grd2(n), vkt(n), v1(n), v2(n), v3(n), v4(n), stat = istat) call chckmem (0, 'xcssb-d', istat) grad = ZERO rho43 = ZERO rhom4 = ZERO x2 = ZERO dp1 = ZERO dp2 = ZERO dp3 = ZERO do1 = ZERO do2 = ZERO do3 = ZERO dt1 = ZERO dt2 = ZERO grd2 = ZERO vkt = ZERO v1 = ZERO v2 = ZERO v3 = ZERO v4 = ZERO tsuml = ZERO tlapl = ZERO factor = ONE if (nspin==1) factor = TWO is_: do is = 1, nspin isplit = min(is,nsplit) do k = 1, n grad(k) = sqrt(drho(k,1,is)**2+drho(k,2,is)**2+drho(k,3,is)**2) grd2(k) = grad(k)*grad(k) rho43(k) = rho(k,is)*rhotrd(k,is) tmpX = max(grad(k)/rho43(k),XTINY) x2(k) = tmpX*tmpX dp1(k) = ONE / (ONE + rC*x2(k)) do1(k) = ONE / (ONE + rE*x2(k)*x2(k)) dt1(k) = ONE / (rho43(k) + DELTA) end do do k = 1, n ! --------------------------- ! only non-LDA part needed !! ! therefore rA - ONE ! --------------------------- if (.not.nilrho(k,is)) e(k,isplit) = e(k,isplit) + factor*Cx*rho43(k)*( & rA - ONE + rB*x2(k)*dp1(k) - rD*x2(k)*do1(k) ) + factor*rG*grd2(k)*dt1(k) end do if (.not.calcv) cycle is_ do k = 1, n dp2(k) = dp1(k)*dp1(k) dp3(k) = dp1(k)*dp2(k) do2(k) = do1(k)*do1(k) do3(k) = do1(k)*do2(k) dt2(k) = dt1(k)*dt1(k) * R4O3 if (.not.nilrho(k,is)) rhom4(k) = ONE / (rho(k,is)**FOUR) end do do k = 1, n ! ----------------------------------------------------------- ! sum over (dRho/dL x d2Rho/dLdK x dRho/dK) where L,K=(x,y,z) ! ----------------------------------------------------------- tsuml(k) = drho(k,1,is)*(drho(k,1,is)*d2rho(k,1,is) + & two*( drho(k,2,is)*d2rho(k,2,is)+drho(k,3,is)*d2rho(k,3,is) ) ) + & drho(k,2,is)*(drho(k,2,is)*d2rho(k,4,is) + two*drho(k,3,is)*d2rho(k,5,is)) + & drho(k,3,is)**2*d2rho(k,6,is) end do do k = 1, n tlapl(k) = d2rho(k,1,is) + d2rho(k,4,is) + d2rho(k,6,is) end do do k = 1, n v1(k) = fV1 * rhotrd(k,is)*( rA + rB*x2(k)*dp1(k)*(ONE - TWO*dp1(k)) + & FOUR*rD*x2(k)*do2(k) - THREE*rD*x2(k)*do1(k)) end do do k = 1, n v2(k) = fV2 * rhotrd(k,is)*( dp3(k)*(TWO*rB*x2(k)*(ONE - THREE*rC*x2(k))) + & 36.0*rD*x2(k)*do2(k) - 32.0*rD*x2(k)*do3(k) - 6.0*rD*x2(k)*do1(k) ) end do do k = 1, n if (.not.nilrho(k,is)) v3(k) = fV3*( rB*dp2(k) -TWO*rD*do2(k) + rD*do1(k) )*tlapl(k)/rho43(k) end do do k = 1, n if (.not.nilrho(k,is)) v4(k) = fV4*( rB*rC*dp3(k) - FOUR*rD*rE*x2(k)*do3(k) + & rD*rE*x2(k)*do2(k) )*tsuml(k)*rhom4(k) end do do k = 1, n vkt(k) = rG*grd2(k)*dt2(k)*rhotrd(k,is) - TWO*rG*tlapl(k)*dt1(k) end do do k = 1, n ! --------------------------- ! only non-LDA part needed !! ! therefore the fV0 term ! --------------------------- if (.not.nilrho(k,is)) v(k,is) = v(k,is) + v1(k) + v2(k) + v3(k) + v4(k) + vkt(k) - fV0*rhotrd(k,is) end do end do is_ deallocate (grad, rho43, rhom4, x2, dp1, dp2, dp3, tsuml, tlapl, v1, v2, v3, v4, & do1, do2, do3, dt1, dt2, grd2, vkt, stat = istat) call chckmem (1, 'xcssb-d', istat) end subroutine xcssbD