Commit 3a7328ee authored by Matthias Redies's avatar Matthias Redies

cut cores und use pbe

parent 5f1d422a
......@@ -45,7 +45,7 @@ CONTAINS
CHARACTER(LEN=20) :: attributes(6)
REAL stateEnergies(29)
! ..
c = c_light(1.0)
seig = 0.
!
......
......@@ -281,11 +281,14 @@ CONTAINS
#ifdef CPP_LIBXC
TYPE(xc_f03_func_info_t) :: xc_info
REAL :: excc(SIZE(exc))
REAL :: cut_ratio = 0.0
INTEGER :: cut_idx
! 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(:,:)
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)
......@@ -300,16 +303,12 @@ CONTAINS
exc=exc+excc
END IF
ELSEIF(xcpot%exc_is_MetaGGA()) THEN
write (*,*) "cut_ration Mgga = ", cut_ratio
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)
!where(kinEnergyDen_libXC < 1d-5) kinEnergyDen_libXC = 1d-5
!write (*,*) "apply tf approx. shapes: "
!write (*,*) "shape(rh) = ", shape(rh)
!write (*,*) "shape(grad%sigma) = ", shape(transpose(grad%sigma))
!write (*,*) "shape(grad%lapl) = ", shape(grad%laplace)
write (filename, '("kED_libxc_", I0.6, ".npy")') size(kinEnergyDen_libxc, dim=2)
call save_npy(filename, transpose(kinEnergyDen_libxc))
......@@ -318,12 +317,33 @@ CONTAINS
exc = 0.0
excc = 0.0
call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, &
transpose(grad%laplace), kinEnergyDen_libXC, exc)
call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh(cut_idx+1:,:),1),&
TRANSPOSE(rh(cut_idx+1:,:)), &
grad%sigma(:,cut_idx+1:), &
transpose(grad%laplace(cut_idx+1:,:)), &
kinEnergyDen_libXC(:,cut_idx+1:), &
exc(cut_idx+1:))
call xc_f03_gga_exc(xcpot%vxc_func_x, SIZE(rh(:cut_idx,:),1), &
TRANSPOSE(rh(:cut_idx,:)), &
grad%sigma(:,:cut_idx), &
exc(:cut_idx))
write (*,*) "max(exc)", maxval(abs(exc(:cut_idx))), maxval(abs(exc(cut_idx+1:)))
IF (xcpot%func_exc_id_c>0) THEN
CALL xc_f03_mgga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, &
transpose(grad%laplace), kinEnergyDen_libXC, excc)
!CALL xc_f03_mgga_exc(xcpot%exc_func_c, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, &
!transpose(grad%laplace), kinEnergyDen_libXC, excc)
call xc_f03_mgga_exc(xcpot%exc_func_c, SIZE(rh(cut_idx+1:,:),1), &
TRANSPOSE(rh(cut_idx+1:,:)), &
grad%sigma(:,cut_idx+1:), &
transpose(grad%laplace(cut_idx+1:,:)), &
kinEnergyDen_libXC(:,cut_idx+1:), &
excc(cut_idx+1:))
call xc_f03_gga_exc(xcpot%vxc_func_c, SIZE(rh(:cut_idx,:),1), &
TRANSPOSE(rh(:cut_idx,:)), &
grad%sigma(:,:cut_idx), &
excc(:cut_idx))
exc=exc+excc
END IF
......
......@@ -58,7 +58,7 @@
INTEGER :: n,nsp,nt,jr
INTEGER :: i, j, idx
REAL :: divi
REAL, PARAMETER :: cut_ratio = 0.1
REAL, PARAMETER :: cut_ratio = 0.9
LOGICAL, allocatable :: cut_mask(:)
! ..
......@@ -147,15 +147,17 @@
CALL mt_from_grid(atoms,sphhar,n,input%jspins,v_x,vx%mt(:,0:,n,:))
! use updated vTot for exc calculation
IF(perform_MetaGGA) THEN
cut_mask = [(.False., i=1,size(atoms%rmsh,dim=1)*nsp)]
idx = 1
do i = 1, floor(size(atoms%rmsh, dim=1)*cut_ratio)
do j = 1,nsp
cut_mask(idx) = .True.
idx = idx +1
enddo
cut_mask = [(.False., i=1,size(atoms%rmsh,dim=1)*nsp)]
idx = 1
do i = 1, floor(size(atoms%rmsh, dim=1)*cut_ratio)
do j = 1,nsp
cut_mask(idx) = .True.
idx = idx +1
enddo
enddo
write (*,*) "set cut_mask", allocated(cut_mask)
IF(perform_MetaGGA) THEN
CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, EnergyDen%mt(:,0:,n,:), &
n, tmp_grad, ED_rs)
......@@ -182,7 +184,9 @@
ELSE
CALL xcpot%get_exc(input%jspins,ch(:nsp*atoms%jri(n),:),e_xc(:nsp*atoms%jri(n),1),grad)
ENDIF
!write (*,*) "cut first ", cut_ratio, " number of points"
!where(cut_mask) e_xc(:,1) = 0.0
IF (lda_atom(n)) THEN
! Use local part of pw91 for this atom
......@@ -195,6 +199,8 @@
nt=nt+nsp
END DO
ENDIF
call save_npy("mt=" // int2str(n) // "_exc.npy", &
get_radial_line(e_xc,1,nsp) )
CALL mt_from_grid(atoms,sphhar,n,1,e_xc,exc%mt(:,0:,n,:))
ENDIF
ENDDO
......
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