Commit 8fdfcdb3 authored by Matthias Redies's avatar Matthias Redies

enable type-based kinED

parent bd5f2ab2
......@@ -18,6 +18,7 @@ MODULE m_types_xcpot
PUBLIC :: t_xcpot,t_gradients
TYPE t_kinED
logical :: set
real, allocatable :: is(:,:) ! (nsp*jmtd, jspins)
real, allocatable :: mt(:,:,:) ! (nsp*jmtd, jspins, local num of types)
contains
......@@ -74,7 +75,7 @@ CONTAINS
class(t_kinED), intent(inout) :: kED
integer, intent(in) :: nsp_x_jmtd, jspins, n_start, n_types, n_stride
integer :: cnt, n
if(.not. allocated(kED%mt)) then
cnt = 0
do n = n_start,n_types,n_stride
......@@ -150,7 +151,7 @@ CONTAINS
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN):: xcpot
xcpot_needs_grad= xcpot%vc_is_gga()
xcpot_needs_grad= xcpot%vc_is_gga() !.or. xcpot%vx_is_MetaGGA()
END FUNCTION xcpot_needs_grad
LOGICAL FUNCTION xcpot_is_hybrid(xcpot)
......@@ -167,7 +168,7 @@ CONTAINS
a_ex=-1
END FUNCTION xcpot_get_exchange_weight
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh,vxc,vx,grad)
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh,vxc,vx,grad, kinED_KS)
USE m_judft
IMPLICIT NONE
......@@ -178,13 +179,14 @@ CONTAINS
!---> xc potential
REAL, INTENT (OUT) :: vxc (:,:),vx(:,:)
TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:)
vxc = 0.0
vx = 0.0
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, mt_call)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinED_KS, mt_call)
USE m_types_misc
USE m_judft
USE, INTRINSIC :: IEEE_ARITHMETIC
......@@ -199,7 +201,7 @@ CONTAINS
REAL, INTENT (OUT) :: exc (:)
TYPE(t_gradients),OPTIONAL,INTENT(IN) :: grad
LOGICAL, OPTIONAL, INTENT(IN) :: mt_call
REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:)
REAL, INTENT(IN), OPTIONAL :: kinED_KS(:,:)
exc = 0.0
call juDFT_error("Can't use XC-parrent class")
......
......@@ -177,7 +177,7 @@ CONTAINS
IF (xcpot%is_name("vhse")) a_ex=amix_hse
END FUNCTION xcpot_get_exchange_weight
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS)
!
USE m_xcxal, ONLY : vxcxal
USE m_xcwgn, ONLY : vxcwgn
......@@ -203,6 +203,7 @@ CONTAINS
!c
REAL, INTENT (OUT) :: vx (:,:)
REAL, INTENT (OUT) :: vxc(:,:)
REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:)
! optional arguments for GGA
TYPE(t_gradients),INTENT(INOUT),OPTIONAL::grad
......@@ -279,7 +280,7 @@ CONTAINS
END SUBROUTINE xcpot_get_vxc
!***********************************************************************
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinEnergyDen_KS, mt_call)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad,kinED_KS, mt_call)
!***********************************************************************
USE m_xcxal, ONLY : excxal
USE m_xcwgn, ONLY : excwgn
......@@ -298,7 +299,7 @@ CONTAINS
REAL, INTENT (OUT) :: exc(:)
TYPE(t_gradients),OPTIONAL,INTENT(IN) ::grad
LOGICAL, OPTIONAL, INTENT(IN) :: mt_call
REAL, INTENT(IN), OPTIONAL :: kinEnergyDen_KS(:,:)
REAL, INTENT(IN), OPTIONAL :: kinED_KS(:,:)
!c
!c ---> local scalars
......
......@@ -30,9 +30,6 @@ MODULE m_types_xcpot_libxc
CONTAINS
!PROCEDURE :: vxc_is_LDA => xcpot_vxc_is_LDA
!PROCEDURE :: vxc_is_gga => xcpot_vxc_is_gga
PROCEDURE :: vx_is_LDA => xcpot_vx_is_LDA
PROCEDURE :: vx_is_GGA => xcpot_vx_is_GGA
PROCEDURE :: vx_is_MetaGGA => xcpot_vx_is_MetaGGA
......@@ -102,12 +99,8 @@ CONTAINS
END IF
!check if any potental is a MetaGGA
IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_x))) THEN
call juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
ELSEIF(xcpot%func_vxc_id_c>0) THEN
IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN
call juDFT_error("vxc_x: MetaGGA is not implemented for potentials")
ENDIF
IF(ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA] == xc_get_family(xcpot%vxc_func_c))) THEN
call juDFT_error("vxc_x: MetaGGA is not implemented for correclation potentials")
ENDIF
CALL write_xc_info(xcpot%vxc_func_x)
......@@ -130,6 +123,8 @@ CONTAINS
ELSE
write (*,*) "Using same functional for VXC and EXC"
END IF
xcpot%kinED%set = .False.
#else
CALL judft_error("You specified a libxc-exchange correlation potential but FLEUR is not linked against libxc", &
hint="Please recompile FLEUR with libxc support")
......@@ -169,8 +164,8 @@ CONTAINS
#ifdef CPP_LIBXC
TYPE(xc_f03_func_info_t) :: xc_info
xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
xcpot_exc_is_LDA = XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info)
xc_info = xc_f03_func_get_info(xcpot%exc_func_x)
xcpot_exc_is_LDA = (XC_FAMILY_LDA == xc_f03_func_info_get_family(xc_info))
#else
xcpot_exc_is_LDA = .false.
#endif
......@@ -209,7 +204,7 @@ CONTAINS
#ifdef CPP_LIBXC
TYPE(xc_f03_func_info_t) :: xc_info
xc_info = xc_f03_func_get_info(xcpot%vxc_func_c)
xc_info = xc_f03_func_get_info(xcpot%vxc_func_x)
xcpot_vx_is_MetaGGA = ANY([XC_FAMILY_MGGA, XC_FAMILY_HYB_MGGA]==xc_f03_func_info_get_family(xc_info))
#else
xcpot_vx_is_MetaGGA=.false.
......@@ -267,7 +262,7 @@ CONTAINS
END FUNCTION xcpot_get_exchange_weight
!***********************************************************************
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad)
SUBROUTINE xcpot_get_vxc(xcpot,jspins,rh, vxc,vx, grad, kinED_KS)
!***********************************************************************
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
......@@ -275,15 +270,18 @@ CONTAINS
REAL,INTENT (IN) :: rh(:,:) !Dimensions here
REAL, INTENT (OUT) :: vx (:,:) !points,spin
REAL, INTENT (OUT ) :: vxc(:,:) !
REAL, INTENT(IN),OPTIONAL:: kinED_KS(:,:)
! optional arguments for GGA
TYPE(t_gradients),OPTIONAL,INTENT(INOUT)::grad
#ifdef CPP_LIBXC
REAL,ALLOCATABLE::vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:)
REAL,ALLOCATABLE :: vxc_tmp(:,:),vx_tmp(:,:),vsigma(:,:), &
tmp_vsig(:,:), tmp_vlapl(:,:), tmp_vtau(:,:), &
kinED_libxc(:,:)
!libxc uses the spin as a first index, hence we have to transpose....
ALLOCATE(vxc_tmp(SIZE(vxc,2),SIZE(vxc,1)));vxc_tmp=0.0
ALLOCATE(vx_tmp(SIZE(vx,2),SIZE(vx,1)));vx_tmp=0.0
IF (xcpot%needs_grad()) THEN
IF (xcpot%vc_is_GGA()) THEN
IF (.NOT.PRESENT(grad)) CALL judft_error("Bug: You called get_vxc for a GGA potential without providing derivatives")
ALLOCATE(vsigma,mold=grad%vsigma)
!where(abs(grad%sigma)<1E-9) grad%sigma=1E-9
......@@ -296,7 +294,22 @@ CONTAINS
vxc_tmp=vx_tmp
ENDIF
ELSE !LDA potentials
CALL xc_f03_lda_vxc(xcpot%vxc_func_x, SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
if(xcpot%vx_is_MetaGGA()) then
if(present(kinED_KS)) then
kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)
allocate(tmp_vsig, mold=vx_tmp)
allocate(tmp_vlapl, mold=vx_tmp)
allocate(tmp_vtau, mold=vx_tmp)
call xc_f03_mgga_vxc(xcpot%vxc_func_x, size(rh,1), transpose(rh), &
grad%sigma, transpose(grad%laplace), kinED_libxc,&
vx_tmp, tmp_vsig, tmp_vlapl, tmp_vtau)
write (*,*) "did some TB09 stuff"
endif
else
CALL xc_f03_lda_vxc(inital_lda_func(jspins), SIZE(rh,1), TRANSPOSE(rh), vx_tmp)
endif
IF (xcpot%func_vxc_id_c>0) THEN
CALL xc_f03_lda_vxc(xcpot%vxc_func_c, SIZE(rh,1), TRANSPOSE(rh), vxc_tmp)
vxc_tmp=vxc_tmp+vx_tmp
......@@ -308,8 +321,7 @@ CONTAINS
#endif
END SUBROUTINE xcpot_get_vxc
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen_KS, mt_call)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinED_KS, mt_call)
use m_constants
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
......@@ -324,7 +336,7 @@ CONTAINS
! 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(:,:)
REAL, INTENT(IN), OPTIONAL :: kinED_KS(:,:)
#ifdef CPP_LIBXC
TYPE(xc_f03_func_info_t) :: xc_info
......@@ -335,7 +347,7 @@ CONTAINS
! 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(:,:)
REAL, ALLOCATABLE :: kinED_libXC(:,:), pkzb_ratio(:,:), pkzb_zaehler(:,:), pkzb_nenner(:,:)
is_mt = merge(mt_call, .False., present(mt_call))
IF (xcpot%exc_is_gga()) THEN
......@@ -352,9 +364,9 @@ CONTAINS
exc=exc+excc
END IF
ELSEIF(xcpot%exc_is_MetaGGA()) THEN
IF(PRESENT(kinEnergyDen_KS)) THEN
IF(PRESENT(kinED_KS)) THEN
! apply correction in eq (4) in https://doi.org/10.1063/1.1565316
kinEnergyDen_libXC = transpose(kinEnergyDen_KS + 0.25 * grad%laplace)
kinED_libXC = transpose(kinED_KS + 0.25 * grad%laplace)
!only cut core of muffin tin
cut_idx = MERGE(NINT(size(rh,1) * cut_ratio), 0, is_mt)
......@@ -365,7 +377,7 @@ CONTAINS
TRANSPOSE(rh(cut_idx+1:,:)), &
grad%sigma(:,cut_idx+1:), &
transpose(grad%laplace(cut_idx+1:,:)), &
kinEnergyDen_libXC(:,cut_idx+1:), &
kinED_libXC(:,cut_idx+1:), &
exc(cut_idx+1:))
call xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh(:cut_idx,:),1), &
......@@ -378,7 +390,7 @@ CONTAINS
TRANSPOSE(rh(cut_idx+1:,:)), &
grad%sigma(:,cut_idx+1:), &
transpose(grad%laplace(cut_idx+1:,:)), &
kinEnergyDen_libXC(:,cut_idx+1:), &
kinED_libXC(:,cut_idx+1:), &
excc(cut_idx+1:))
call xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh(:cut_idx,:),1), &
......@@ -478,6 +490,19 @@ CONTAINS
family = xc_f03_func_info_get_family(xc_f03_func_get_info(xc_func))
END FUNCTION xc_get_family
function inital_lda_func(jspins) result(lda)
use, intrinsic :: iso_c_binding
implicit none
INTEGER, intent(in) :: jspins
TYPE(xc_f03_func_t) :: lda
if(jspins == 1) then
CALL xc_f03_func_init(lda, 1, XC_UNPOLARIZED)
else
CALL xc_f03_func_init(lda, 1, XC_POLARIZED)
ENDIF
end function inital_lda_func
#endif
END MODULE m_types_xcpot_libxc
......@@ -37,7 +37,7 @@ CONTAINS
USE m_metagga
IMPLICIT NONE
CLASS(t_xcpot),INTENT(INOUT) :: xcpot
CLASS(t_xcpot),INTENT(INOUT) :: xcpot
TYPE(t_input),INTENT(IN) :: input
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
......@@ -52,7 +52,10 @@ CONTAINS
COMPLEX, ALLOCATABLE :: den_pw_w(:,:), EnergyDen_pw_w(:,:), vtot_pw_norm(:,:)
REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:,:)
INTEGER :: jspin, i, js
LOGICAL :: perform_MetaGGA
perform_MetaGGA = ALLOCATED(EnergyDen%mt) &
.AND. (xcpot%exc_is_MetaGGA() .or. xcpot%vx_is_MetaGGA())
CALL init_pw_grid(xcpot,stars,sym,cell)
!Put the charge on the grid, in GGA case also calculate gradients
......@@ -61,7 +64,11 @@ CONTAINS
ALLOCATE(v_xc,mold=rho)
ALLOCATE(v_x,mold=rho)
CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad)
if(perform_MetaGGA .and. xcpot%kinED%set) then
CALL xcpot%get_vxc(input%jspins,rho,v_xc, v_x,grad, kinED_KS=xcpot%kinED%is)
else
CALL xcpot%get_vxc(input%jspins,rho,v_xc,v_x,grad)
endif
IF (xcpot%needs_grad()) THEN
SELECT TYPE(xcpot)
......@@ -85,6 +92,7 @@ CONTAINS
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, vTot%pw, tmp_grad, vTot_rs)
CALL calc_kinEnergyDen_pw(ED_rs, vTot_rs, rho, xcpot%kinED%is)
xcpot%kinED%set = .True.
ENDIF
!calculate the ex.-cor energy density
......
......@@ -121,7 +121,15 @@
!
! calculate the ex.-cor. potential
CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:),v_x(:nsp*atoms%jri(n),:),grad)
write (*,*) "perform_MGGA = ", perform_MetaGGA
write (*,*) "xcpot%kinED%set = ", xcpot%kinED%set
if(perform_MetaGGA .and. xcpot%kinED%set) then
CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:)&
, v_x(:nsp*atoms%jri(n),:),grad, kinED_KS=xcpot%kinED%mt(:,:,loc_n))
else
CALL xcpot%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),v_xc(:nsp*atoms%jri(n),:)&
, v_x(:nsp*atoms%jri(n),:),grad)
endif
IF (lda_atom(n)) THEN
! Use local part of pw91 for this atom
CALL xcpot_tmp%get_vxc(input%jspins,ch(:nsp*atoms%jri(n),:),xcl(:nsp*atoms%jri(n),:),v_x(:nsp*atoms%jri(n),:),grad)
......@@ -167,6 +175,7 @@
xcpot%val_den%mt(:,0:,n,:), n, tmp_grad, val_den_rs)
CALL calc_kinEnergyDen_mt(ED_rs, vTot_rs, vTot0_rs, &
core_den_rs, val_den_rs, n, nsp, xcpot%kinED%mt(:,:,loc_n))
xcpot%kinED%set = .True.
ENDIF
IF (ALLOCATED(exc%mt)) THEN
......@@ -174,9 +183,10 @@
! calculate the ex.-cor energy density
!
IF(perform_MetaGGA) THEN
IF(perform_MetaGGA .and. xcpot%kinED%set) THEN
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),&
e_xc(:nsp*atoms%jri(n),1),grad, xcpot%kinED%mt(:,:,loc_n), mt_call=.True.)
e_xc(:nsp*atoms%jri(n),1),grad, &
kinED_KS=xcpot%kinED%mt(:,:,loc_n), mt_call=.True.)
ELSE
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),&
e_xc(:nsp*atoms%jri(n),1),grad, mt_call=.True.)
......
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