Commit bb5381be authored by Matthias Redies's avatar Matthias Redies

BUGFIX: only cut exc in MT, not IS

parent 58008f04
......@@ -162,7 +162,7 @@ CONTAINS
call juDFT_error("Can't use XC-parrent class")
END SUBROUTINE xcpot_get_vxc
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS, mt_call)
USE m_types_misc
USE m_judft
USE, INTRINSIC :: IEEE_ARITHMETIC
......@@ -176,6 +176,7 @@ CONTAINS
!---> xc energy density
REAL, INTENT (OUT) :: exc (:)
TYPE(t_gradients),OPTIONAL,INTENT(IN) :: grad
LOGICAL, OPTIONAL, INTENT(IN) :: mt_call
REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:)
exc = 0.0
......
......@@ -279,7 +279,7 @@ CONTAINS
END SUBROUTINE xcpot_get_vxc
!***********************************************************************
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS, mt_call)
!***********************************************************************
USE m_xcxal, ONLY : excxal
USE m_xcwgn, ONLY : excwgn
......@@ -297,6 +297,7 @@ CONTAINS
REAL,INTENT (IN) :: rh(:,:)
REAL, INTENT (OUT) :: exc(:)
TYPE(t_gradients),OPTIONAL,INTENT(IN) ::grad
LOGICAL, OPTIONAL, INTENT(IN) :: mt_call
REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:)
!c
......
......@@ -309,34 +309,35 @@ CONTAINS
END SUBROUTINE xcpot_get_vxc
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen_KS)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen_KS, mt_call)
use m_constants
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins
REAL,INTENT (IN) :: rh(:,:) !points,spin
REAL, INTENT (OUT) :: exc(:) !points
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins
REAL,INTENT (IN) :: rh(:,:) !points,spin
REAL, INTENT (OUT) :: exc(:) !points
! optional arguments for GGA
TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
TYPE(t_gradients),OPTIONAL,INTENT(IN) :: grad
LOGICAL, OPTIONAL, INTENT(IN) :: mt_call
! kinED from Kohn-Sham equations:
! tau = sum[phi_i(r)^dag nabla phi_i(r)]
! see eq (2) in https://doi.org/10.1063/1.1565316
! (-0.5 is applied below)
REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:)
character(len=200) :: filename
#ifdef CPP_LIBXC
TYPE(xc_f03_func_info_t) :: xc_info
REAL :: excc(SIZE(exc))
REAL :: cut_ratio = 0.1
INTEGER :: cut_idx
LOGICAL :: is_mt
! tau = 0.5 * sum[|grad phi_i(r)|²]
! see eq (3) in https://doi.org/10.1063/1.1565316
REAL, ALLOCATABLE :: kinEnergyDen_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
is_mt = merge(mt_call, .False., present(mt_call))
IF (xcpot%exc_is_gga()) THEN
IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_exc for a GGA potential without providing derivatives")
CALL xc_f03_gga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh),grad%sigma,exc)
......@@ -352,10 +353,12 @@ CONTAINS
END IF
ELSEIF(xcpot%exc_is_MetaGGA()) THEN
IF(PRESENT(kinEnergyDen_KS)) THEN
cut_idx = NINT(size(rh,1) * cut_ratio)
! apply correction in eq (4) in https://doi.org/10.1063/1.1565316
kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25 * grad%laplace)
!only cut core of muffin tin
cut_idx = MERGE(NINT(size(rh,1) * cut_ratio), 0, is_mt)
exc = 0.0
excc = 0.0
call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh(cut_idx+1:,:),1),&
......
......@@ -92,9 +92,9 @@ CONTAINS
ALLOCATE ( e_xc(SIZE(rho,1),1) ); e_xc=0.0
IF(ALLOCATED(EnergyDen%pw) .AND. xcpot%exc_is_MetaGGA()) THEN
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, kinED_rs)
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, kinED_rs, mt_call=.False.)
ELSE
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad)
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, mt_call=.False.)
ENDIF
CALL pw_from_grid(xcpot,stars,.TRUE.,e_xc,exc%pw,exc%pw_w)
ENDIF
......
......@@ -172,9 +172,11 @@
!
IF(perform_MetaGGA) THEN
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),e_xc(:nsp*atoms%jri(n),1),grad, kinED_rs)
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),&
e_xc(:nsp*atoms%jri(n),1),grad, kinED_rs, mt_call=.True.)
ELSE
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),e_xc(:nsp*atoms%jri(n),1),grad)
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),&
e_xc(:nsp*atoms%jri(n),1),grad, mt_call=.True.)
ENDIF
!write (*,*) "cut first ", cut_ratio, " number of points"
......
......@@ -110,7 +110,7 @@ CONTAINS
!i calculate the exchange-correlation energy density in real space
!
IF (ALLOCATED(exc%vacz)) THEN
call xcpot%get_exc(input%jspins,af2,e_xc)
call xcpot%get_exc(input%jspins,af2,e_xc, mt_call=.False.)
!
! ----> 2-d back fft to g space
!
......@@ -159,7 +159,7 @@ CONTAINS
! calculate the ex.-corr. energy density now beyond warping region
!
IF (ALLOCATED(exc%vacz)) THEN
CALL xcpot%get_exc(input%jspins,af2(0:nmzdiff-1,:),exc%vacz(vacuum%nmzxy+1:,ivac,1))
CALL xcpot%get_exc(input%jspins,af2(0:nmzdiff-1,:),exc%vacz(vacuum%nmzxy+1:,ivac,1), mt_call=.False.)
END IF
ENDDO
IF (noco%l_noco) THEN
......
......@@ -422,7 +422,7 @@ CONTAINS
!
IF (ALLOCATED(exc%vacz)) THEN
CALL xcpot%get_exc(input%jspins,af2,e_xc,grad)
CALL xcpot%get_exc(input%jspins,af2,e_xc,grad, mt_call=.False.)
! ----> 2-d back fft to g space
!
......@@ -545,7 +545,7 @@ CONTAINS
! calculate the ex-corr. energy density now beyond warping region
!
IF (ALLOCATED(exc%vacz)) THEN
CALL xcpot%get_exc(input%jspins,rhtxc(:nmzdiff,:),exc%vacz(vacuum%nmzxy+1:,ivac,1),grad)
CALL xcpot%get_exc(input%jspins,rhtxc(:nmzdiff,:),exc%vacz(vacuum%nmzxy+1:,ivac,1),grad, mt_call=.False.)
ENDIF
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment