Commit d223655c authored by Matthias Redies's avatar Matthias Redies

apply correction to KS-kinED

parents 3a4ae40d ff6793c4
......@@ -136,7 +136,7 @@
CALL juDFT_error("mapatom",calledby="mapatom")
ENDIF
WRITE (6,FMT=8010) nat1,na,atoms%ngopr(na)
8010 FORMAT (5x,'atom',i3,' can be mapped into atom',i3,&
8010 FORMAT (5x,'atom',i5,' can be mapped into atom',i5,&
& ' through group operation',i4)
!
! end of equivalent atoms
......
......@@ -22,6 +22,7 @@ MODULE m_xmlOutput
INTEGER, SAVE :: currentElementIndex
INTEGER, SAVE :: maxNumElements
INTEGER, SAVE :: xmlOutputUnit
LOGICAL, SAVE :: xmlOpened = .FALSE.
CHARACTER(LEN= 40), ALLOCATABLE :: elementList(:)
PUBLIC startXMLOutput, endXMLOutput
......@@ -86,6 +87,7 @@ MODULE m_xmlOutput
ELSE
OPEN (xmlOutputUnit,file='out.xml',form='formatted',status='unknown')
ENDIF
xmlOpened = .TRUE.
WRITE (xmlOutputUnit,'(a)') '<?xml version="1.0" encoding="UTF-8" standalone="no"?>'
WRITE (xmlOutputUnit,'(a)') '<fleurOutput fleurOutputVersion="0.27">'
CALL openXMLElement('programVersion',(/'version'/),(/version_const/))
......@@ -150,6 +152,7 @@ MODULE m_xmlOutput
CHARACTER(LEN=10) :: dateString
CHARACTER(LEN=10) :: timeString
IF (.NOT.xmlOpened) RETURN
DO WHILE (currentElementIndex.NE.0)
CALL closeXMLElement(elementList(currentElementIndex))
END DO
......
......@@ -245,7 +245,7 @@ CONTAINS
ENDIF
IF (time<min_time*globaltimer%time) RETURN !do not print parts that take less than min_time of the totaltime
IF (timer%no_calls>1) THEN
WRITE(fid,"(a,T7,a,T46,a,T85,a,i0,a,f9.0,a,f9.0,a)") TRIM(indentstring),TRIM(timername),TRIM(timestring(time,globaltimer%time,level)),"(",timer%no_calls,"calls:",timer%mintime,"sec -",timer%maxtime,"sec)"
WRITE(fid,"(a,T7,a,T46,a,T85,a,i0,a,f9.3,a,f9.3,a)") TRIM(indentstring),TRIM(timername),TRIM(timestring(time,globaltimer%time,level)),"% (",timer%no_calls,"calls:",timer%mintime,"sec -",timer%maxtime,"sec)"
ELSE
WRITE(fid,"(a,T7,a,T46,a)") TRIM(indentstring),TRIM(timername),timestring(time,globaltimer%time,level)
END IF
......@@ -304,6 +304,7 @@ CONTAINS
OPEN(2,FILE ='juDFT_times',STATUS="replace")
ENDIF
globaltimer%time=cputime()-globaltimer%starttime
globaltimer%starttime=cputime()
WRITE(fn,"('Total execution time: ',i0,'sec')") INT(globaltimer%time)
CALL add_usage_data("Runtime",globaltimer%time)
CALL priv_writetimes_longest(globaltimer,fid=fn)
......
......@@ -259,7 +259,7 @@ CONTAINS
END SUBROUTINE xcpot_get_vxc
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen)
SUBROUTINE xcpot_get_exc(xcpot,jspins,rh,exc,grad, kinEnergyDen_KS)
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins
......@@ -267,7 +267,16 @@ CONTAINS
REAL, INTENT (OUT) :: exc(:) !points
! optional arguments for GGA
TYPE(t_gradients),OPTIONAL,INTENT(IN)::grad
REAL, INTENT(IN), OPTIONAL :: kinEnergyDen(:,:)
! 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(:,:)
! tau = 0.5 * sum[|grad phi_i(r)|²]
! see eq (3) in https://doi.org/10.1063/1.1565316
REAL, ALLOCATABLE :: kinEnergyDen_libXC(:,:)
#ifdef CPP_LIBXC
TYPE(xc_f03_func_info_t) :: xc_info
......@@ -288,12 +297,14 @@ CONTAINS
END IF
ELSEIF(xcpot%exc_is_MetaGGA()) THEN
IF(PRESENT(kinEnergyDen)) THEN
! apply correction in eq (4) in https://doi.org/10.1063/1.1565316
kinEnergyDen_libXC = transpose(0.25 * grad%laplace - 0.5 * kinEnergyDen_KS)
call xc_f03_mgga_exc(xcpot%exc_func_x, SIZE(rh,1), TRANSPOSE(rh), grad%sigma, &
transpose(grad%laplace), transpose(kinEnergyDen), exc)
transpose(grad%laplace), kinEnergyDen_libXC, exc)
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), transpose(kinEnergyDen), excc)
transpose(grad%laplace), kinEnergyDen_libXC, excc)
exc=exc+excc
END IF
......
......@@ -51,9 +51,9 @@ CONTAINS
IF (noco%l_mtnocoPot) CALL rotate_mt_den_from_local(atoms,sphhar,sym,denRot,vtot)
ENDIF
write (*,*) "Set vTot to zero in vgen_finalize()"
vTot%pw_w = 0.0
vTot%pw = 0.0
!write (*,*) "Set vTot to zero in vgen_finalize()"
!vTot%pw_w = 0.0
!vTot%pw = 0.0
! Rescale vCoul%pw_w with number of stars
DO js = 1, SIZE(vCoul%pw_w,2)
......
......@@ -89,10 +89,10 @@ CONTAINS
vTot_pw_norm(i,js)=vTot%pw(i,js) / stars%nstr(i)
END DO
END DO
write (*,*) "set vtot = 0 again in vis_xc"
vTot_pw_norm = 0.0
vTot%pw = 0.0
vTot%pw_w = 0.0
!write (*,*) "set vtot = 0 again in vis_xc"
!vTot_pw_norm = 0.0
!vTot%pw = 0.0
!vTot%pw_w = 0.0
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, vTot_pw_norm, tmp_grad, vTot_rs)
......
......@@ -5,121 +5,89 @@
!--------------------------------------------------------------------------------
MODULE m_libxc_postprocess_gga
CONTAINS
SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad)
USE m_mt_tofrom_grid
USE m_types
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sphhar),INTENT(IN) :: sphhar
INTEGER,INTENT(IN) :: n
REAL,INTENT(INOUT) :: v_xc(:,:)
TYPE(t_gradients),INTENT(IN):: grad
SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sphhar,n,v_xc,grad)
USE m_mt_tofrom_grid
USE m_types
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sphhar),INTENT(IN) :: sphhar
INTEGER,INTENT(IN) :: n
REAL,INTENT(INOUT) :: v_xc(:,:)
TYPE(t_gradients),INTENT(IN):: grad
INTEGER :: nsp,n_sigma,i
REAL,ALLOCATABLE:: vsigma(:,:),vsigma_mt(:,:,:)
TYPE(t_gradients)::grad_vsigma
n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise
nsp=SIZE(v_xc,1) !no of points
ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:sphhar%nlhd,n_sigma))
vsigma_mt=0.0
vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
CALL mt_from_grid(atoms,sphhar,nsp/atoms%jmtd,n,n_sigma,vsigma,vsigma_mt)
DO i=1,atoms%jri(n)
vsigma_mt(i,:,:)=vsigma_mt(i,:,:)*atoms%rmsh(i,n)**2
ENDDO
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL mt_to_grid(xcpot,n_sigma,atoms,sphhar,vsigma_mt,nsp/atoms%jmtd,n,grad=grad_vsigma)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
END SUBROUTINE libxc_postprocess_gga_mt
INTEGER :: nsp,n_sigma,i
REAL,ALLOCATABLE:: vsigma(:,:),vsigma_mt(:,:,:)
TYPE(t_gradients)::grad_vsigma
SUBROUTINE libxc_postprocess_gga_pw(xcpot,stars,cell,v_xc,grad)
USE m_pw_tofrom_grid
USE m_types
USE m_npy
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
REAL,INTENT(INOUT) :: v_xc(:,:)
TYPE(t_gradients),INTENT(IN):: grad
n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise
nsp=SIZE(v_xc,1) !no of points
ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:sphhar%nlhd,n_sigma))
vsigma_mt=0.0
vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
CALL mt_from_grid(atoms,sphhar,nsp/atoms%jmtd,n,n_sigma,vsigma,vsigma_mt)
DO i=1,atoms%jri(n)
vsigma_mt(i,:,:)=vsigma_mt(i,:,:)*atoms%rmsh(i,n)**2
ENDDO
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL mt_to_grid(xcpot,n_sigma,atoms,sphhar,vsigma_mt,nsp/atoms%jmtd,n,grad=grad_vsigma)
COMPLEX,ALLOCATABLE:: vsigma_g(:,:)
REAL,ALLOCATABLE:: vsigma(:,:)
TYPE(t_gradients)::grad_vsigma
INTEGER :: nsp,n_sigma
nsp=SIZE(v_xc,1) !no of points
n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine
ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma));vsigma_g=0.0
vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g)
!vsigma_g(:,1)=vsigma_g(:,1)*stars%nstr(:)
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL pw_to_grid(xcpot,n_sigma,.false.,stars,cell,vsigma_g,grad_vsigma)
call save_npy("lapl.npy", grad%laplace)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
END SUBROUTINE libxc_postprocess_gga_mt
END SUBROUTINE libxc_postprocess_gga_pw
SUBROUTINE libxc_postprocess_gga_pw(xcpot,stars,cell,v_xc,grad)
USE m_pw_tofrom_grid
USE m_types
USE m_npy
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
REAL,INTENT(INOUT) :: v_xc(:,:)
TYPE(t_gradients),INTENT(IN):: grad
SUBROUTINE libxc_postprocess_gga_pw_alt(xcpot,stars,cell,v_xc,grad)
USE m_pw_tofrom_grid
USE m_types
IMPLICIT NONE
CLASS(t_xcpot),INTENT(IN) :: xcpot
TYPE(t_stars),INTENT(IN) :: stars
TYPE(t_cell),INTENT(IN) :: cell
REAL,INTENT(INOUT) :: v_xc(:,:)
TYPE(t_gradients),INTENT(IN):: grad
COMPLEX,ALLOCATABLE:: vsigma_g(:,:)
REAL,ALLOCATABLE:: vsigma(:,:)
TYPE(t_gradients)::grad_vsigma
INTEGER :: nsp,n_sigma
COMPLEX,ALLOCATABLE:: vsigma_g(:,:)
REAL,ALLOCATABLE:: vsigma(:,:)
TYPE(t_gradients)::grad_vsigma
INTEGER :: nsp,n_sigma,i
nsp=SIZE(v_xc,1) !no of points
n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine
ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma));vsigma_g=0.0
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
DO i=1,3
vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
!Multiply with gradient
vsigma(:,1)=vsigma(:,1)*grad%gr(i,:,1)
CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g)
CALL pw_to_grid(xcpot,n_sigma,.FALSE.,stars,cell,vsigma_g,grad_vsigma)
v_xc(:,1)=v_xc(:,1)-2*grad_vsigma%gr(i,:,1)
ENDDO
END SUBROUTINE libxc_postprocess_gga_pw_alt
SUBROUTINE libxc_postprocess_gga(vsigma,grad,grad_vsigma,v_xc)
USE m_types
IMPLICIT NONE
REAL,INTENT(IN) :: vsigma(:,:)
TYPE(t_gradients),INTENT(IN):: grad,grad_vsigma
REAL,INTENT(INOUT) :: v_xc(:,:)
INTEGER:: i
IF (SIZE(v_xc,2)==1) THEN !Single spin
DO i=1,SIZE(v_xc,1) !loop over points
v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)
ENDDO
ELSE !two spins
DO i=1,SIZE(v_xc,1) !loop over points
v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)-&
dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,2))-vsigma(i,2)*grad%laplace(i,2)
v_xc(i,2)=v_xc(i,2)-2*dot_PRODUCT(grad_vsigma%gr(:,i,3),grad%gr(:,i,2))-2*vsigma(i,3)*grad%laplace(i,2)-&
dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,1))-vsigma(i,2)*grad%laplace(i,1)
ENDDO
END IF
END SUBROUTINE libxc_postprocess_gga
nsp=SIZE(v_xc,1) !no of points
n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine
ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma)); vsigma_g=0.0
vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
CALL pw_from_grid(xcpot,stars,.FALSE.,vsigma,vsigma_g)
!vsigma_g(:,1)=vsigma_g(:,1)*stars%nstr(:)
ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
CALL pw_to_grid(xcpot,n_sigma,.false.,stars,cell,vsigma_g,grad_vsigma)
call save_npy("lapl.npy", grad%laplace)
call save_npy("grad.npy", grad%gr)
CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
END SUBROUTINE libxc_postprocess_gga_pw
SUBROUTINE libxc_postprocess_gga(vsigma,grad,grad_vsigma,v_xc)
USE m_types
IMPLICIT NONE
REAL,INTENT(IN) :: vsigma(:,:)
TYPE(t_gradients),INTENT(IN):: grad,grad_vsigma
REAL,INTENT(INOUT) :: v_xc(:,:)
INTEGER:: i
IF (SIZE(v_xc,2)==1) THEN !Single spin
DO i=1,SIZE(v_xc,1) !loop over points
v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)
ENDDO
ELSE !two spins
DO i=1,SIZE(v_xc,1) !loop over points
v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)-&
dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,2))-vsigma(i,2)*grad%laplace(i,2)
v_xc(i,2)=v_xc(i,2)-2*dot_PRODUCT(grad_vsigma%gr(:,i,3),grad%gr(:,i,2))-2*vsigma(i,3)*grad%laplace(i,2)-&
dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,1))-vsigma(i,2)*grad%laplace(i,1)
ENDDO
END IF
END SUBROUTINE libxc_postprocess_gga
END MODULE m_libxc_postprocess_gga
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