Commit fc0eff0b authored by Matthias Redies's avatar Matthias Redies

Merge branch 'TranBlaha' into 'develop'

Tran blaha

See merge request fleur/fleur!28
parents 01c6e734 7dfccdc9
......@@ -74,6 +74,54 @@ CONTAINS
END DO ! loop over spins
END SUBROUTINE integrate_cdn
SUBROUTINE integrate_realspace(xcpot, atoms, sym, sphhar, input, &
stars, cell, oneD, vacuum, noco, mt, is, hint)
use m_types
use m_mt_tofrom_grid
use m_pw_tofrom_grid
use m_constants
implicit none
CLASS(t_xcpot), INTENT(inout) :: xcpot
TYPE(t_atoms),INTENT(IN) :: atoms
TYPE(t_sym), INTENT(in) :: sym
TYPE(t_sphhar), INTENT(IN) :: sphhar
TYPE(t_input), INTENT(IN) :: input
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_oneD), INTENT(in) :: oneD
TYPE(t_vacuum), INTENT(in) :: vacuum
TYPE(t_noco), INTENT(in) :: noco
real, intent(inout) :: mt(:,:,:), is(:,:)
character(len=*), intent(in), optional :: hint
integer :: n_atm, i
TYPE(t_potden) :: tmp_potden
REAL :: q(input%jspins), qis(input%jspins), &
qmt(atoms%ntype,input%jspins), qvac(2,input%jspins),&
qtot, qistot
call tmp_potden%init(stars, atoms, sphhar, vacuum, noco, input%jspins, POTDEN_TYPE_DEN)
call init_mt_grid(input%jspins, atoms, sphhar, xcpot, sym)
do n_atm =1,atoms%ntype
call mt_from_grid(atoms, sphhar, n_atm, input%jspins, mt(:,:,n_atm), &
tmp_potden%mt(:,0:,n_atm,:))
do i=1,atoms%jri(n_atm)
tmp_potden%mt(i,:,n_atm,:) = tmp_potden%mt(i,:,n_atm,:) * atoms%rmsh(i,n_atm)**2
enddo
enddo
call finish_mt_grid()
call init_pw_grid(xcpot, stars, sym, cell)
call pw_from_grid(xcpot, stars, .False., is, tmp_potden%pw)
call finish_pw_grid()
call integrate_cdn(stars,atoms,sym,vacuum,input,cell,oneD, tmp_potden, &
q, qis, qmt, qvac, qtot, qistot)
call print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
END SUBROUTINE integrate_realspace
SUBROUTINE cdntot(stars,atoms,sym,vacuum,input,cell,oneD,&
den,l_printData,qtot,qistot)
......@@ -154,4 +202,29 @@ CONTAINS
CALL timestop("cdntot")
END SUBROUTINE cdntot
SUBROUTINE print_cdn_inte(q, qis, qmt, qvac, qtot, qistot, hint)
use ieee_arithmetic
implicit none
REAL, INTENT(in) :: q(:), qis(:), qmt(:,:), qvac(:,:), qtot, qistot
character(len=*), intent(in), optional :: hint
integer :: n_mt
if(present(hint)) write (*,*) "DEN of ", hint
write (*,*) "q = ", q
write (*,*) "qis = ", qis
write (*,*) "qmt"
do n_mt = 1,size(qmt, dim=1)
write (*,*) "mt = ", n_mt, qmt(n_mt,:)
enddo
if(.not. any(ieee_is_nan(qvac))) then
write (*, *) "qvac", qvac
endif
write (*, *) "qtot", qtot
write (*, *) "qis_tot", qistot
write (*, *) "-------------------------"
END SUBROUTINE print_cdn_inte
END MODULE m_cdntot
Subproject commit 3cb2231abf1d47fbd8b3e21c8478e9f26a73ce5f
Subproject commit ca6f7114b9fffe0964ca9e8d24e09ef15b300316
......@@ -81,6 +81,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
TYPE(t_slab) :: slab
TYPE(t_orbcomp) :: orbcomp
TYPE(t_cdnvalJob) :: cdnvalJob
TYPE(t_potden) :: val_den, core_den
!Local Scalars
......@@ -91,7 +92,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
#ifdef CPP_HDF
INTEGER(HID_T) :: banddosFile_id
#endif
LOGICAL :: l_error
LOGICAL :: l_error, perform_MetaGGA
CALL regCharges%init(input,atoms)
CALL dos%init(input,atoms,dimension,kpts,vacuum)
......@@ -117,7 +118,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
CALL cdnval(eig_id,mpi,kpts,jspin,noco,input,banddos,cell,atoms,enpara,stars,vacuum,dimension,&
sphhar,sym,vTot,oneD,cdnvalJob,outDen,regCharges,dos,results,moments,coreSpecInput,mcd,slab,orbcomp)
END DO
call xcpot%val_den%copyPotDen(outDen)
call val_den%copyPotDen(outDen)
! calculate kinetic energy density for MetaGGAs
if(xcpot%exc_is_metagga()) then
......@@ -170,7 +171,7 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
CALL cdncore(mpi,dimension,oneD,input,vacuum,noco,sym,&
stars,cell,sphhar,atoms,vTot,outDen,moments,results)
endif
call xcpot%core_den%subPotDen(outDen, xcpot%val_den)
call core_den%subPotDen(outDen, val_den)
CALL timestop("cdngen: cdncore")
CALL enpara%calcOutParams(input,atoms,vacuum,regCharges)
......@@ -198,7 +199,13 @@ SUBROUTINE cdngen(eig_id,mpi,input,banddos,sliceplot,vacuum,&
END IF
END IF ! mpi%irank == 0
perform_MetaGGA = ALLOCATED(EnergyDen%mt) &
.AND. (xcpot%exc_is_MetaGGA() .or. xcpot%vx_is_MetaGGA())
if(perform_MetaGGA) then
call set_kinED(mpi, sphhar, atoms, sym, core_den, val_den, xcpot, &
input, noco, stars, cell, outDen, EnergyDen, vTot)
endif
#ifdef CPP_MPI
CALL MPI_BCAST(noco%l_ss,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
CALL MPI_BCAST(noco%l_mperp,1,MPI_LOGICAL,0,mpi%mpi_comm,ierr)
......
......@@ -415,7 +415,7 @@ CONTAINS
! mix input and output densities
CALL mix_charge(field2,DIMENSION,mpi,(iter==input%itmax.OR.judft_was_argument("-mix_io")),&
stars,atoms,sphhar,vacuum,input,&
sym,cell,noco,oneD,archiveType,inDen,outDen,results)
sym,cell,noco,oneD,archiveType,xcpot,iter,inDen,outDen,results)
IF(mpi%irank == 0) THEN
WRITE (6,FMT=8130) iter
......
......@@ -17,7 +17,7 @@ contains
SUBROUTINE mix_charge( field, DIMENSION, mpi, l_writehistory,&
stars, atoms, sphhar, vacuum, input, sym, cell, noco, &
oneD, archiveType, inDen, outDen, results )
oneD, archiveType, xcpot, iteration, inDen, outDen, results )
use m_juDFT
use m_constants
......@@ -47,10 +47,11 @@ contains
type(t_dimension), intent(in) :: dimension
type(t_mpi), intent(in) :: mpi
TYPE(t_atoms),TARGET,INTENT(in) :: atoms
class(t_xcpot), intent(in) :: xcpot
type(t_potden), intent(inout) :: outDen
type(t_results), intent(inout) :: results
type(t_potden), intent(inout) :: inDen
integer, intent(in) :: archiveType
integer, intent(in) :: archiveType, iteration
LOGICAL, INTENT(IN) :: l_writehistory
real :: fix
......@@ -129,7 +130,7 @@ contains
CALL judft_error("Unknown Mixing schema")
END SELECT
CALL timestop("Mixing")
CALL timestart("Postprocessing")
!extracte mixed density
......@@ -146,6 +147,11 @@ contains
CALL mixvector_reset()
ENDIF
if(iteration == 1 .and. xcpot%vx_is_MetaGGA()) then
CALL mixing_history_reset(mpi)
CALL mixvector_reset()
endif
!fix charge of the new density
IF (mpi%irank==0) CALL qfix(mpi,stars,atoms,sym,vacuum, sphhar,input,cell,oneD,inDen,noco%l_noco,.FALSE.,.FALSE., fix)
......
......@@ -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
......@@ -26,7 +27,6 @@ MODULE m_types_xcpot
TYPE,ABSTRACT :: t_xcpot
REAL :: gmaxxc
TYPE(t_potden) :: core_den, val_den
TYPE(t_kinED) :: kinED
CONTAINS
PROCEDURE :: vxc_is_LDA => xcpot_vxc_is_LDA
......@@ -74,13 +74,13 @@ 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
cnt = cnt + 1
enddo
allocate(kED%mt(nsp_x_jmtd, jspins, cnt))
allocate(kED%mt(nsp_x_jmtd, jspins, cnt), source=0.0)
endif
end subroutine kED_alloc_mt
......@@ -150,7 +150,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 +167,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 +178,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 +200,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
......
......@@ -27,12 +27,7 @@ MODULE m_types_xcpot_libxc
INTEGER :: func_vxc_id_c, func_vxc_id_x !> functionals to be used for potential & density convergence
INTEGER :: func_exc_id_c, func_exc_id_x !> functionals to be used in exc- & totale-calculations
INTEGER :: jspins
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
......@@ -112,12 +107,8 @@ CONTAINS
IF(errors(4) /= 0) call juDFT_error("Correlation energy functional not in LibXC")
!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)
......@@ -140,6 +131,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")
......@@ -179,8 +172,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
......@@ -219,7 +212,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.
......@@ -277,23 +270,27 @@ 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)
USE, INTRINSIC :: IEEE_ARITHMETIC
IMPLICIT NONE
CLASS(t_xcpot_libxc),INTENT(IN) :: xcpot
INTEGER, INTENT (IN) :: jspins
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(:,:)
integer :: idx
!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
......@@ -306,7 +303,24 @@ 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() .and. 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)
idx = find_first_normal(vx_tmp)+1
vx_tmp(:,:idx) = 0.0
CALL xc_f03_lda_vxc(initial_lda_func(jspins), idx, TRANSPOSE(rh(:idx,:)), vx_tmp(:,:idx))
else
CALL xc_f03_lda_vxc(initial_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
......@@ -314,12 +328,24 @@ CONTAINS
ENDIF
vx=TRANSPOSE(vx_tmp)
vxc=TRANSPOSE(vxc_tmp)
#endif
END SUBROUTINE xcpot_get_vxc
FUNCTION find_first_normal(vec) result(idx)
USE, INTRINSIC :: IEEE_ARITHMETIC
implicit none
real, intent(in) :: vec(:,:)
integer :: idx
idx = size(vec, dim=2)
do while(all(.not. ieee_is_nan(vec(:,idx))))
idx = idx - 1
if(idx == 0) exit
enddo
END FUNCTION find_first_normal
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
......@@ -334,7 +360,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
......@@ -345,7 +371,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
......@@ -362,9 +388,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)
......@@ -375,7 +401,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), &
......@@ -388,7 +414,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), &
......@@ -424,6 +450,10 @@ CONTAINS
ALLOCATE(grad%laplace(ngrid,jspins))
ALLOCATE(grad%vsigma(MERGE(1,3,jspins==1),ngrid))
grad%gr = 0.0
grad%laplace = 0.0
grad%sigma = 0.0
grad%vsigma = 0.0
END SUBROUTINE xcpot_alloc_gradients
#ifdef CPP_LIBXC
......@@ -488,6 +518,23 @@ CONTAINS
family = xc_f03_func_info_get_family(xc_f03_func_get_info(xc_func))
END FUNCTION xc_get_family
function initial_lda_func(jspins) result(lda)
use, intrinsic :: iso_c_binding
implicit none
INTEGER, intent(in) :: jspins
TYPE(xc_f03_func_t) :: lda
integer :: ierr
if(jspins == 1) then
CALL xc_f03_func_init(lda, 1, XC_UNPOLARIZED, err=ierr)
else
CALL xc_f03_func_init(lda, 1, XC_POLARIZED, err=ierr)
ENDIF
if(ierr /= 0) call juDFT_error("can't find lda_x for init")
end function initial_lda_func
#endif
END MODULE m_types_xcpot_libxc
This diff is collapsed.
......@@ -51,10 +51,6 @@ CONTAINS
IF (noco%l_mtnocoPot) CALL rotate_mt_den_from_local(atoms,sphhar,sym,denRot,vtot)
ENDIF
!write (*,*) "Set vTot%is 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)
DO i = 1, stars%ng3
......
......@@ -32,6 +32,7 @@ CONTAINS
USE m_checkdopall
USE m_cdn_io
USE m_convol
use m_cdntot
IMPLICIT NONE
......@@ -56,6 +57,7 @@ CONTAINS
! Local type instances
TYPE(t_potden) :: workDen, exc, veff
real, allocatable :: tmp_mt(:,:,:), tmp_is(:,:)
! Local Scalars
INTEGER ifftd, ifftd2, ifftxc3d, ispin, i
#ifdef CPP_MPI
......@@ -133,10 +135,8 @@ CONTAINS
CALL timestart("Vxc in MT")
END IF
CALL vmt_xc(DIMENSION, mpi, sphhar, atoms, den, xcpot, input, sym, &
obsolete, EnergyDen, vTot, vx, exc)
!
CALL vmt_xc(mpi, sphhar, atoms, den, xcpot, input, sym, &
EnergyDen, vTot, vx, exc)
! add MT EXX potential to vr
IF (mpi%irank == 0) THEN
......
......@@ -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
......@@ -49,10 +49,12 @@ CONTAINS
TYPE(t_gradients) :: grad, tmp_grad
REAL, ALLOCATABLE :: rho(:,:), ED_rs(:,:), vTot_rs(:,:)
REAL, ALLOCATABLE :: rho_conv(:,:), ED_conv(:,:), vTot_conv(:,:)
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 +63,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)
......@@ -73,25 +79,11 @@ CONTAINS
CALL pw_from_grid(xcpot,stars,input%total,v_xc,vTot%pw,vTot%pw_w)
CALL pw_from_grid(xcpot,stars,input%total,v_x,vx%pw,vx%pw_w)
! use updated vTot for exc calculation
IF(ALLOCATED(EnergyDen%pw) .AND. xcpot%exc_is_MetaGGA()) THEN
allocate(den_pw_w, mold=vTot%pw_w)
allocate(EnergyDen_pw_w, mold=vTot%pw_w)
allocate(vtot_pw_norm, mold=vTot%pw)
CALL pw_to_grid(xcpot, input%jspins, noco%l_noco, stars, &
cell, EnergyDen%pw, tmp_grad, ED_rs)
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)
ENDIF
!calculate the ex.-cor energy density
IF (ALLOCATED(exc%pw_w)) THEN
ALLOCATE ( e_xc(SIZE(rho,1),1) ); e_xc=0.0
IF(ALLOCATED(EnergyDen%pw) .AND. xcpot%exc_is_MetaGGA()) THEN
IF(xcpot%kinED%set) THEN
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, xcpot%kinED%is, mt_call=.False.)
ELSE
CALL xcpot%get_exc(input%jspins,rho,e_xc(:,1),grad, mt_call=.False.)
......@@ -100,6 +92,5 @@ CONTAINS
ENDIF
CALL finish_pw_grid()
END SUBROUTINE vis_xc
END MODULE m_vis_xc
......@@ -23,8 +23,8 @@
! *********************************************************
CONTAINS
SUBROUTINE vmt_xc(DIMENSION,mpi,sphhar,atoms,&
den,xcpot,input,sym, obsolete,EnergyDen,vTot,vx,exc)
SUBROUTINE vmt_xc(mpi,sphhar,atoms,&
den,xcpot,input,sym,EnergyDen,vTot,vx,exc)
#include"cpp_double.h"
use m_libxc_postprocess_gga
USE m_mt_tofrom_grid
......@@ -34,10 +34,8 @@
USE m_juDFT_string
IMPLICIT NONE
CLASS(t_xcpot),INTENT(INOUT) :: xcpot
TYPE(t_dimension),INTENT(IN) :: dimension
CLASS(t_xcpot),INTENT(INOUT) :: xcpot
TYPE(t_mpi),INTENT(IN) :: mpi
TYPE(t_obsolete),INTENT(IN) :: obsolete
TYPE(t_input),INTENT(IN) :: input
TYPE(t_sym),INTENT(IN) :: sym
TYPE(t_sphhar),INTENT(IN) :: sphhar
......@@ -53,8 +51,7 @@
TYPE(t_xcpot_inbuild) :: xcpot_tmp
TYPE(t_potden) :: vTot_tmp
TYPE(t_sphhar) :: tmp_sphhar
REAL, ALLOCATABLE :: ch(:,:), core_den_rs(:,:), val_den_rs(:,:), ED_rs(:,:), &
vTot_rs(:,:), vTot0_rs(:,:)
REAL, ALLOCATABLE :: ch(:,:)
INTEGER :: n,nsp,nt,jr, loc_n
INTEGER :: i, j, idx, cnt
REAL :: divi
......@@ -90,15 +87,6 @@
ALLOCATE(ch(nsp*atoms%jmtd,input%jspins))
IF (xcpot%needs_grad()) CALL xcpot%alloc_gradients(SIZE(ch,1),input%jspins,grad)
IF (perform_MetaGGA) THEN
IF (xcpot%needs_grad()) CALL xcpot%alloc_gradients(SIZE(ch,1),input%jspins,tmp_grad)
ALLOCATE(ED_rs, mold=ch)
ALLOCATE(vTot_rs, mold=ch)
ALLOCATE(vTot0_rs, mold=vTot_rs)
ALLOCATE(core_den_rs, mold=ch)
ALLOCATE(val_den_rs, mold=ch)
ENDIF
CALL init_mt_grid(input%jspins,atoms,sphhar,xcpot,sym)
#ifdef CPP_MPI
......@@ -121,7 +109,13 @@
!
! 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)
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)
......@@ -144,39 +138,16 @@
CALL mt_from_grid(atoms,sphhar,n,input%jspins,v_xc,vTot%mt(:,0:,n,:))
CALL mt_from_grid(atoms,sphhar,n,input%jspins,v_x,vx%mt(:,0:,n,:))
IF(perform_MetaGGA) THEN
CALL mt_to_grid(xcpot, input%jspins, atoms, sphhar, EnergyDen%mt(:,0:,n,:), &
n, tmp_grad, ED_rs)