scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

metaGGA vectorized and libxc incorporated

This commit is contained in:
fabrizio22 2019-06-24 13:42:51 +02:00
parent 687b16c0bd
commit 4a8bd335b5
12 changed files with 1913 additions and 1796 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CPV/src/make.depend
View File

@ -717,6 +717,7 @@ manycp.o : ../../Modules/read_input.o
manycp.o : input.o
metaxc.o : ../../Modules/funct.o
metaxc.o : ../../Modules/kind.o
metaxc.o : ../../Modules/xc_mgga_drivers.o
modules.o : ../../Modules/kind.o
modules.o : ../../Modules/uspp.o
move_electrons.o : ../../Modules/cell_base.o

242
CPV/src/metaxc.f90
View File

@ -10,190 +10,96 @@
SUBROUTINE tpssmeta(nnr, nspin,grho,rho,kedtau,etxc)
! ===================
!--------------------------------------------------------------------
use kinds, only: dp
use funct, only: tau_xc, tau_xc_spin, tau_xc_array, tau_xc_array_spin, get_meta
use kinds, only: dp
use funct, only: get_meta
use xc_mgga, only: xc_metagcx, change_threshold_mgga !, &
!tau_xc_array, tau_xc_array_spin
!
IMPLICIT NONE
!
! input
integer nspin , nnr
real(dp) grho(3,nnr,nspin), rho(nnr,nspin),kedtau(nnr,nspin)
integer :: nspin, nnr
real(dp) :: grho(3,nnr,nspin), rho(nnr,nspin), kedtau(nnr,nspin)
! output: excrho: exc * rho ; E_xc = \int excrho(r) d_r
! output: rhor: contains the exchange-correlation potential
real(dp) etxc
real(dp) :: etxc
REAL(dp) :: zeta, rh, grh2
INTEGER :: k, ipol, is
REAL(dp), PARAMETER :: epsr = 1.0d-6, epsg = 1.0d-10
INTEGER :: imeta
!
etxc = 0.d0
! calculate the gradient of rho+rho_core in real space
imeta = get_meta()
if (imeta.eq.5.or.imeta.eq.6.or.imeta.eq.7) then
call exch_corr_meta_array_mode() !HK/MCA: currently only implmented for SCAN
else
call exch_corr_meta_scalar_mode() !HK/MCA: compatibility for the original implementation
end if
!
call exch_corr_meta() !HK/MCA
!
RETURN
contains
subroutine exch_corr_meta_array_mode()
implicit none
real(dp) :: grho_(3,nnr,nspin) !MCA/HK : store grho only in nspin=2
REAL(dp) :: arho(nnr), segno(nnr), grho2 (nnr), &
& sx(nnr), sc(nnr), &
& v1x(nnr,nspin), v2x(nnr,nspin*2-1), v3x(nnr,nspin), & !MCA/HK
& v1c(nnr,nspin), v2c(nnr,nspin*2-1), v3c(nnr,nspin) !MCA/HK
IF (nspin == 1) THEN
!
!$omp parallel do
do k = 1, nnr
!
grho2(k) = grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2
arho(k) = ABS (rho (k,1) )
segno(k) = SIGN (1.d0, rho (k,1) )
!
end do !k
!$omp end parallel do
!
CALL tau_xc_array (nnr,arho,grho2,kedtau,sx,sc,v1x,v2x,v3x,v1c,v2c,v3c)
!
! store potentials
!
rho (:, 1) = ( v1x(:,1) + v1c(:,1) )
kedtau(:,1) = ( v3x(:,1) + v3c(:,1) ) *0.5_dp
!
! v2 contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
DO ipol = 1, 3
grho(ipol,:,1) = ( v2x(:,1) + v2c(:,1) )*grho (ipol,:,1)
ENDDO
!
ELSE
!
!MCA/HK: only SCAN is available
CALL tau_xc_array_spin (nnr, rho, grho, kedtau, sx, sc, v1x, v2x, &
& v3x, v1c, v2c, v3c)
!
! MCA/HK : store grho to compute v2x cross terms
!
grho_ = grho
!
DO is = 1,nspin
!
rho(:, is) = v1x(:,is) + v1c(:,is)
!
DO ipol = 1, 3 !MCA/HK: second line is the cross term
grho(ipol,:,is) = ( v2x(:,2*is-1) + v2c(:,2*is-1) ) * grho(ipol,:,is) &
& + 0.5_dp * ( v2x(:,2) + v2c(:,2) ) * grho_(ipol,:,MOD(is,2)+1)
ENDDO
!
kedtau(:,is)= ( v3x(:,is) + v3c(:,is) ) *0.5d0
!
segno = 1.0 !MCA: not the most efficient way
!
ENDDO
!
ENDIF
!
!
contains
!
!
subroutine exch_corr_meta()
!
! compute exc energy contribution from the current process
!
etxc = 0.0_dp
!$omp parallel do reduction(+:etxc)
do k = 1, nnr
etxc = etxc + (sx(k) + sc(k)) * segno(k)
end do !k
!$omp end parallel do
return
end subroutine exch_corr_meta_array_mode
subroutine exch_corr_meta_scalar_mode()
implicit none
REAL(dp) :: grho2 (2), sx, sc, v1x, v2x, v3x,v1c, v2c, v3c, &
v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw ,v2cup(3),v2cdw(3), &
v3xup, v3xdw,grhoup(3),grhodw(3),v3cup, v3cdw, segno, arho, atau
DO k = 1, nnr
DO is = 1, nspin
grho2 (is) = grho(1,k, is)**2 + grho(2,k,is)**2 + grho(3,k,is)**2
ENDDO
IF (nspin == 1) THEN
!
! This is the spin-unpolarised case
!
arho = ABS (rho (k, 1) )
segno = SIGN (1.d0, rho (k, 1) )
atau = kedtau(k,1)
IF (arho.GT.epsr.AND.grho2 (1) .GT.epsg.AND.ABS(atau).GT.epsr) THEN
CALL tau_xc (arho, grho2(1), atau, sx, sc, &
v1x, v2x, v3x, v1c, v2c, v3c)
rho (k, 1) = (v1x + v1c )
kedtau(k,1)= (v3x + v3c) *0.5d0
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
DO ipol = 1, 3
grho(ipol,k,1) = (v2x + v2c)*grho (ipol,k,1)
ENDDO
etxc = etxc + (sx + sc) * segno
ELSE
DO ipol = 1, 3
grho (ipol,k,1) = 0.d0
ENDDO
kedtau(k,1)=0.d0
ENDIF
ELSE
!
! spin-polarised case
!
!CALL tpsscx_spin(rho (k, 1), rho (k, 2), grho2 (1), grho2 (2), &
! kedtau(k,1),kedtau(k,2),sx, &
! v1xup,v1xdw,v2xup,v2xdw,v3xup,v3xdw)
rh = rho (k, 1) + rho (k, 2)
IF (rh.GT.epsr) THEN
!zeta = (rho (k, 1) - rho (k, 2) ) / rh
DO ipol=1,3
grhoup(ipol)=grho(ipol,k,1)
grhodw(ipol)=grho(ipol,k,2)
END DO
! atau=kedtau(k,1)+kedtau(k,2)
call tau_xc_spin (rho(k,1), rho(k,2), grhoup, grhodw, &
kedtau(k,1), kedtau(k,2), sx, sc, v1xup, v1xdw, v2xup, v2xdw, &
v3xup, v3xdw, v1cup, v1cdw, v2cup, v2cdw,&
v3cup, v3cdw)
!CALL tpsscc_spin(rh,zeta,grhoup,grhodw, &
! atau,sc,v1cup,v1cdw,v2cup,v2cdw,v3c)
ELSE
sx = 0.d0
sc = 0.d0
v1xup = 0.d0
v1xdw = 0.d0
v2xup=0.d0
v2xdw=0.d0
v3xup=0.d0
v3xdw=0.d0
v1cup = 0.d0
v1cdw = 0.d0
v2cup=0.d0
v2cdw=0.d0
v3cup=0.d0
v3cdw=0.d0
!
ENDIF
!
! first term of the gradient correction : D(rho*Exc)/D(rho)
!
rho(k, 1) = (v1xup + v1cup)
rho(k, 2) = (v1xdw + v1cdw)
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
DO ipol = 1, 3
grho(ipol,k,1) = (v2xup*grho(ipol,k,1) + v2cup(ipol))
grho(ipol,k,2) = (v2xdw*grho(ipol,k,2) + v2cdw(ipol))
ENDDO
kedtau(k,1)= (v3xup + v3cup) *0.5d0
kedtau(k,2)= (v3xdw + v3cdw) *0.5d0
etxc = etxc + (sx + sc)
ENDIF
ENDDO
!
INTEGER :: np
REAL(DP), ALLOCATABLE :: sx(:), v1x(:,:), v2x(:,:), v3x(:,:), &
sc(:), v1c(:,:), v2c(:,:,:), v3c(:,:)
!
np=1
if (nspin==2) np=3
!
allocate( sx(nnr), v1x(nnr,nspin), v2x(nnr,nspin), v3x(nnr,nspin) )
allocate( sc(nnr), v1c(nnr,nspin), v2c(np,nnr,nspin), v3c(nnr,nspin) )
!
if (nspin==1) then
!
call change_threshold_mgga( epsr, epsg, epsr )
!
call xc_metagcx( nnr, 1, np, rho, grho, kedtau, sx, sc, &
v1x, v2x, v3x, v1c, v2c, v3c )
!
rho(:,1) = v1x(:,1) + v1c(:,1)
kedtau(:,1) = (v3x(:,1) + v3c(:,1)) * 0.5d0
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
do ipol = 1, 3
grho(ipol,:,1) = (v2x(:,1) + v2c(1,:,1))*grho(ipol,:,1)
enddo
etxc = SUM( (sx(:) + sc(:)) * SIGN(1.d0,rho(:,1)) )
!
else
!
call change_threshold_mgga( epsr )
!
call xc_metagcx( nnr, 2, np, rho, grho, kedtau, sx, sc, &
v1x, v2x, v3x, v1c, v2c, v3c )
!
rho(:,1) = v1x(:,1) + v1c(:,1)
rho(:,2) = v1x(:,2) + v1c(:,2)
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
do ipol = 1, 3
grho(ipol,:,1) = (v2x(:,1)*grho(ipol,:,1) + v2c(ipol,:,1))
grho(ipol,:,2) = (v2x(:,2)*grho(ipol,:,2) + v2c(ipol,:,2))
enddo
!
kedtau(:,1) = (v3x(:,1) + v3c(:,1)) * 0.5d0
kedtau(:,2) = (v3x(:,2) + v3c(:,2)) * 0.5d0
etxc = etxc + SUM( sx(:) + sc(:) )
!
endif
!
deallocate( sx, v1x, v2x, v3x )
deallocate( sc, v1c, v2c, v3c )
!
!
return
end subroutine exch_corr_meta_scalar_mode
!
end subroutine exch_corr_meta
!
END SUBROUTINE tpssmeta
!-----------------------------------------------------------------------

1
Modules/Makefile
View File

@ -111,6 +111,7 @@ xc_vdW_DF.o \
xc_rVV10.o \
xc_lda_lsda_drivers.o \
xc_gga_drivers.o \
xc_mgga_drivers.o \
io_base.o \
qes_types_module.o \
qes_libs_module.o \

377
Modules/funct.f90
View File

@ -50,7 +50,7 @@ MODULE funct
PUBLIC :: enforce_input_dft, write_dft_name
PUBLIC :: get_dft_name, get_dft_short, get_dft_long,&
get_nonlocc_name
PUBLIC :: get_iexch, get_icorr, get_igcx, get_igcc, get_meta, get_inlc
PUBLIC :: get_iexch, get_icorr, get_igcx, get_igcc, get_meta, get_metac, get_inlc
PUBLIC :: dft_is_gradient, dft_is_meta, dft_is_hybrid, dft_is_nonlocc, igcc_is_lyp
PUBLIC :: set_auxiliary_flags
!
@ -67,8 +67,6 @@ MODULE funct
!
! driver subroutines computing XC
PUBLIC :: init_xc, is_libxc
PUBLIC :: tau_xc , tau_xc_spin
PUBLIC :: tau_xc_array, tau_xc_array_spin
PUBLIC :: nlc
!
! PRIVATE variables defining the DFT functional
@ -327,10 +325,11 @@ MODULE funct
INTEGER :: igcx = notset
INTEGER :: igcc = notset
INTEGER :: imeta = notset
INTEGER :: imetac= notset
INTEGER :: inlc = notset
!
! is_libxc(i)==.TRUE. if the the i-th term of xc is from libxc
LOGICAL :: is_libxc(6)
LOGICAL :: is_libxc(7)
!
REAL(DP):: exx_fraction = 0.0_DP
REAL(DP):: screening_parameter = 0.0_DP
@ -390,7 +389,8 @@ CONTAINS
CHARACTER(len=50):: dftout
LOGICAL :: dft_defined = .FALSE.
CHARACTER(LEN=1), EXTERNAL :: capital
INTEGER :: save_iexch, save_icorr, save_igcx, save_igcc, save_meta, save_inlc
INTEGER :: save_iexch, save_icorr, save_igcx, save_igcc, save_meta, &
save_metac, save_inlc
#if defined(__LIBXC)
INTEGER :: fkind
TYPE(xc_f90_pointer_t) :: xc_func, xc_info
@ -407,6 +407,7 @@ CONTAINS
save_igcx = igcx
save_igcc = igcc
save_meta = imeta
save_metac = imetac
save_inlc = inlc
!
! convert to uppercase
@ -621,6 +622,26 @@ CONTAINS
END SELECT
!
!
! ... A temporary fix in order to keep the q-e input notation for SCAN-functionls
! valid.
IF (imeta==5 .OR. imeta==6) THEN
#if defined(__LIBXC)
imeta = 263
imetac = 267
is_libxc(5:6) = .TRUE.
#else
CALL errore( 'set_dft_from_name', 'libxc needed for this functional', 2 )
#endif
ELSEIF (imeta==3) THEN
#if defined(__LIBXC)
imeta = 208
imetac = 231
is_libxc(5:6) = .TRUE.
#else
CALL errore( 'set_dft_from_name', 'libxc needed for this functional', 2 )
#endif
ENDIF
!
!----------------------------------------------------------------
! If the DFT was not yet defined, check every part of the string
!----------------------------------------------------------------
@ -632,7 +653,12 @@ CONTAINS
igcx = matching( 3, dftout, ngcx, gradx, is_libxc(3) )
igcc = matching( 4, dftout, ngcc, gradc, is_libxc(4) )
imeta = matching( 5, dftout, nmeta, meta, is_libxc(5) )
inlc = matching( 6, dftout, ncnl, nonlocc, is_libxc(6) )
IF ( is_libxc(5) ) THEN
imetac = matching( 6, dftout, nmeta, meta, is_libxc(6) )
ELSE
imetac = 0
ENDIF
inlc = matching( 7, dftout, ncnl, nonlocc, is_libxc(7) )
!
#if defined(__LIBXC)
fkind = -100
@ -660,6 +686,10 @@ CONTAINS
IF (ANY(is_libxc(1:2)) .AND. ANY(is_libxc(3:4))) &
CALL errore( 'set_dft_from_name', 'An LDA functional has been found, but &
&libxc GGA functionals already include the LDA part)', 4 )
!
IF (imeta/=0 .AND. (.NOT. is_libxc(5)) .AND. imetac/=0) &
CALL errore( 'set_dft_from_name', 'Two conflicting metaGGA functionals &
&have been found', 5 )
#endif
!
ENDIF
@ -718,6 +748,10 @@ CONTAINS
WRITE (stdout,*) inlc, save_meta
CALL errore( 'set_dft_from_name', ' conflicting values for imeta', 1 )
ENDIF
IF (save_metac /= notset .AND. save_metac /= imetac) THEN
WRITE (stdout,*) imetac, save_metac
CALL errore( 'set_dft_from_name', ' conflicting values for imetac', 1 )
ENDIF
IF (save_inlc /= notset .AND. save_inlc /= inlc) THEN
WRITE (stdout,*) inlc, save_inlc
CALL errore( 'set_dft_from_name', ' conflicting values for inlc', 1 )
@ -835,9 +869,16 @@ CONTAINS
IF (family==XC_FAMILY_HYB_GGA .AND. fkind==XC_EXCHANGE_CORRELATION) RETURN
CASE( 4 )
IF (family==XC_FAMILY_GGA .AND. fkind==XC_CORRELATION) RETURN
CASE( 5 )
IF (family==XC_FAMILY_MGGA .AND. fkind==XC_EXCHANGE) RETURN
IF (family==XC_FAMILY_MGGA .AND. fkind==XC_EXCHANGE_CORRELATION) RETURN
IF (family==XC_FAMILY_HYB_MGGA .AND. fkind==XC_EXCHANGE) RETURN
IF (family==XC_FAMILY_HYB_MGGA .AND. fkind==XC_EXCHANGE_CORRELATION) RETURN
CASE( 6 )
IF (family==XC_FAMILY_MGGA .AND. fkind==XC_CORRELATION) RETURN
END SELECT
!
#endif
!
slot_match_libxc=.FALSE.
!
RETURN
@ -903,6 +944,7 @@ CONTAINS
igcc = i4
inlc = i5
imeta = i6
imetac= 0
!
set_dft_values = .TRUE.
!
@ -1077,6 +1119,12 @@ CONTAINS
get_meta = imeta
RETURN
END FUNCTION get_meta
!
FUNCTION get_metac()
INTEGER get_metac
get_metac = imetac
RETURN
END FUNCTION get_metac
!-----------------------------------------------------------------------
FUNCTION get_inlc()
INTEGER get_inlc
@ -1218,23 +1266,23 @@ CONTAINS
!
IF ( iexch==1 .AND. igcx==0 .AND. igcc==0) THEN
shortname = TRIM(corr(icorr))
ELSEIF (iexch==4 .AND. icorr== 0 .AND. igcx==0 .AND. igcc== 0) THEN
ELSEIF (iexch==4 .AND. icorr==0 .AND. igcx==0 .AND. igcc== 0) THEN
shortname = 'OEP'
ELSEIF (iexch==1 .AND. icorr==11 .AND. igcx==0 .AND. igcc== 0) THEN
shortname = 'VWN-RPA'
ELSEIF (iexch==1 .AND. icorr== 3 .AND. igcx==1 .AND. igcc== 3) THEN
ELSEIF (iexch==1 .AND. icorr==3 .AND. igcx==1 .AND. igcc== 3) THEN
shortname = 'BLYP'
ELSEIF (iexch==1 .AND. icorr== 1 .AND. igcx==1 .AND. igcc== 0) THEN
ELSEIF (iexch==1 .AND. icorr==1 .AND. igcx==1 .AND. igcc== 0) THEN
shortname = 'B88'
ELSEIF (iexch==1 .AND. icorr== 1 .AND. igcx==1 .AND. igcc== 1) THEN
ELSEIF (iexch==1 .AND. icorr==1 .AND. igcx==1 .AND. igcc== 1) THEN
shortname = 'BP'
ELSEIF (iexch==1 .AND. icorr== 4 .AND. igcx==2 .AND. igcc== 2) THEN
ELSEIF (iexch==1 .AND. icorr==4 .AND. igcx==2 .AND. igcc== 2) THEN
shortname = 'PW91'
ELSEIF (iexch==1 .AND. icorr== 4 .AND. igcx==3 .AND. igcc== 4) THEN
ELSEIF (iexch==1 .AND. icorr==4 .AND. igcx==3 .AND. igcc== 4) THEN
shortname = 'PBE'
ELSEIF (iexch==6 .AND. icorr== 4 .AND. igcx==8 .AND. igcc== 4) THEN
ELSEIF (iexch==6 .AND. icorr==4 .AND. igcx==8 .AND. igcc== 4) THEN
shortname = 'PBE0'
ELSEIF (iexch==6 .AND. icorr== 4 .AND. igcx==41.AND. igcc== 4) THEN
ELSEIF (iexch==6 .AND. icorr==4 .AND. igcx==41.AND. igcc== 4) THEN
shortname = 'B86BPBEX'
ELSEIF (iexch==6 .AND. icorr==4 .AND. igcx==42.AND. igcc== 3) THEN
shortname = 'BHANDHLYP'
@ -1371,9 +1419,11 @@ SUBROUTINE init_xc( family )
USE xc_lda_lsda, ONLY: libxc_switches_lda, iexch_l, icorr_l, &
exx_started_l, is_there_finite_size_corr, &
exx_fraction_l, finite_size_cell_volume_l
USE xc_gga, ONLY: libxc_switches_gga, igcx_l, igcc_l, &
exx_started_g, exx_fraction_g, &
USE xc_gga, ONLY: libxc_switches_gga, igcx_l, igcc_l, &
exx_started_g, exx_fraction_g, &
screening_parameter_l, gau_parameter_l
USE xc_mgga, ONLY: libxc_switches_mgga, imeta_l, imetac_l, &
exx_started_mg, exx_fraction_mg
!
IMPLICIT NONE
!
@ -1387,8 +1437,8 @@ SUBROUTINE init_xc( family )
iexch_l = get_iexch()
icorr_l = get_icorr()
!
IF (iexch_l==-1 .OR. icorr_l==-1) CALL errore( 'init_xc', 'LDA functional &
& indexes not well defined', 1 )
IF (iexch_l==notset .OR. icorr_l==notset) CALL errore( 'init_xc', 'LDA functional &
& indexes not defined', 1 )
!
! hybrid exchange vars
exx_started_l = exx_started !is_active()
@ -1408,8 +1458,8 @@ SUBROUTINE init_xc( family )
igcx_l = get_igcx()
igcc_l = get_igcc()
!
IF (igcx_l==-1 .OR. igcc_l==-1) CALL errore( 'init_xc', 'GGA functional &
& indexes not well defined', 2 )
IF (igcx_l==notset .OR. igcc_l==notset) CALL errore( 'init_xc', 'GGA functional &
& indexes not defined', 2 )
!
! hybrid exchange vars
exx_started_g = exx_started !is_active()
@ -1420,8 +1470,26 @@ SUBROUTINE init_xc( family )
gau_parameter_l = get_gau_parameter()
ENDIF
!
IF (family.NE.'LDA' .AND. family.NE.'GGA' .AND. family.NE.'ALL') &
CALL errore( 'init_xc', 'family not found', 3 )
IF (family.EQ.'MGGA' .OR. family.EQ.'ALL') THEN
! =1 if libxc active, =0 otherwise
IF (is_libxc(5)) libxc_switches_mgga(1) = 1
IF (is_libxc(6)) libxc_switches_mgga(2) = 1
! exchange-correlation indexes
imeta_l = get_meta()
imetac_l = get_metac()
!
IF (imeta_l==notset .OR. imetac_l==notset) CALL errore( 'init_xc', 'MGGA functional &
& indexes not defined', 3 )
!
! hybrid exchange vars
exx_started_mg = exx_started !is_active()
exx_fraction_mg = 0._DP
IF ( exx_started_mg ) exx_fraction_mg = get_exx_fraction()
!
ENDIF
!
IF ( family.NE.'LDA' .AND. family.NE.'GGA' .AND. family.NE.'MGGA' .AND. family.NE.'ALL') &
CALL errore( 'init_xc', 'family not found', 4 )
!
RETURN
!
@ -1479,269 +1547,6 @@ SUBROUTINE nlc (rho_valence, rho_core, nspin, enl, vnl, v)
!
RETURN
END SUBROUTINE nlc
!
!-----------------------------------------------------------------------
!------- META CORRECTIONS DRIVERS ----------------------------------
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE tau_xc (rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
!-----------------------------------------------------------------------
! gradient corrections for exchange and correlation - Hartree a.u.
! See comments at the beginning of module for implemented cases
!
! input: rho, grho=|\nabla rho|^2
!
! definition: E_x = \int e_x(rho,grho) dr
!
! output: sx = e_x(rho,grho) = grad corr
! v1x= D(E_x)/D(rho)
! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
! v3x= D(E_x)/D(tau)
!
! sc, v1c, v2c as above for correlation
!
IMPLICIT NONE
REAL(DP) :: rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c
!_________________________________________________________________________
if (imeta == 1) then
CALL tpsscxc (rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
elseif (imeta == 2) then
CALL m06lxc (rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
elseif (imeta == 3) then
CALL tb09cxc (rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
elseif (imeta == 4) then
! do nothing
elseif (imeta == 5) then
CALL SCANcxc (rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
else
CALL errore('tau_xc','wrong igcx and/or igcc',1)
end if
RETURN
END SUBROUTINE tau_xc
SUBROUTINE tau_xc_array (nnr, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
! HK/MCA : the xc_func_init is slow and is called too many times
! HK/MCA : we modify this subroutine so that the overhead could be minimized
!-----------------------------------------------------------------------
! gradient corrections for exchange and correlation - Hartree a.u.
! See comments at the beginning of module for implemented cases
!
! input: rho, grho=|\nabla rho|^2
!
! definition: E_x = \int e_x(rho,grho) dr
!
! output: sx = e_x(rho,grho) = grad corr
! v1x= D(E_x)/D(rho)
! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
! v3x= D(E_x)/D(tau)
!
! sc, v1c, v2c as above for correlation
!
IMPLICIT NONE
INTEGER, intent(in) :: nnr
REAL(DP) :: rho(nnr), grho(nnr), tau(nnr), ex(nnr), ec(nnr)
REAL(DP) :: v1x(nnr), v2x(nnr), v3x(nnr), v1c(nnr), v2c(nnr), v3c(nnr)
!_________________________________________________________________________
if (imeta == 5) then
CALL scancxc_array (nnr, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
elseif (imeta == 6 ) then ! HK/MCA: SCAN0
CALL scancxc_array (nnr, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
if (exx_started) then
ex = (1.0_DP - exx_fraction) * ex
v1x = (1.0_DP - exx_fraction) * v1x
v2x = (1.0_DP - exx_fraction) * v2x
v3x = (1.0_DP - exx_fraction) * v3x
end if
else
CALL errore('v_xc_meta_array','(CP only) array mode only works for SCAN',1)
end if
RETURN
END SUBROUTINE tau_xc_array
!
!
!-----------------------------------------------------------------------
SUBROUTINE tau_xc_spin (rhoup, rhodw, grhoup, grhodw, tauup, taudw, ex, ec, &
& v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw, v1cup, v1cdw,&
& v2cup, v2cdw, v3cup, v3cdw)
!-----------------------------------------------------------------------
!
!
IMPLICIT NONE
real(dp), intent(in) :: rhoup, rhodw, tauup, taudw
real(dp), dimension (3), intent(in) :: grhoup, grhodw
real(dp), intent(out) :: ex, ec, v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw, &
& v1cup, v1cdw, v3cup, v3cdw
real(dp), dimension(3), intent(out) :: v2cup, v2cdw
!
! Local variables
!
INTEGER :: ipol
real(dp) :: rh, zeta, atau, grhoup2, grhodw2
real(dp), parameter :: epsr=1.0d-08, zero=0._dp
!
!_____________________________
grhoup2 = zero
grhodw2 = zero
v2cup = zero
v2cdw = zero
! FIXME: for SCAN, this will be calculated later
if (imeta /= 4) then
do ipol=1,3
grhoup2 = grhoup2 + grhoup(ipol)**2
grhodw2 = grhodw2 + grhodw(ipol)**2
end do
end if
if (imeta == 1) then
CALL tpsscx_spin(rhoup, rhodw, grhoup2, grhodw2, tauup, &
& taudw, ex, v1xup,v1xdw,v2xup,v2xdw,v3xup,v3xdw)
rh = rhoup + rhodw
zeta = (rhoup - rhodw) / rh
atau = tauup + taudw ! KE-density in Hartree
CALL tpsscc_spin(rh,zeta,grhoup,grhodw, atau,ec, &
& v1cup,v1cdw,v2cup,v2cdw,v3cup, v3cdw)
elseif (imeta == 2) then
CALL m06lxc_spin (rhoup, rhodw, grhoup2, grhodw2, tauup, taudw, &
& ex, ec, v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw, &
& v1cup, v1cdw, v2cup(1), v2cdw(1), v3cup, v3cdw)
elseif (imeta == 5) then
! FIXME: not the most efficient use of libxc
CALL scanxc_spin(rhoup, rhodw, grhoup, grhodw, tauup, taudw, &
& ex, v1xup,v1xdw,v2xup,v2xdw,v3xup,v3xdw, &
& ec, v1cup,v1cdw,v2cup,v2cdw,v3cup,v3cdw )
else
CALL errore('tau_xc_spin','This case not implemented',imeta)
end if
END SUBROUTINE tau_xc_spin
SUBROUTINE tau_xc_array_spin (nnr, rho, grho, tau, ex, ec, v1x, v2x, v3x, &
& v1c, v2c, v3c)
! HK/MCA : the xc_func_init (LIBXC) is slow and is called too many times
! HK/MCA : we modify this SUBROUTINE so that the overhead could be minimized
!-----------------------------------------------------------------------
! gradient corrections for exchange and correlation - Hartree a.u.
! See comments at the beginning of module for implemented cases
!
! input: rho,rho, grho=\nabla rho
!
! definition: E_x = \int e_x(rho,grho) dr
!
! output: sx = e_x(rho,grho) = grad corr
! v1x= D(E_x)/D(rho)
! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
! v3x= D(E_x)/D(tau)
!
! sc, v1cup, v2cup as above for correlation
!
IMPLICIT NONE
INTEGER, intent(in) :: nnr
REAL(DP) :: rho(nnr,2), grho(3,nnr,2), tau(nnr,2), ex(nnr), ec(nnr)
REAL(DP) :: v1x(nnr,2), v2x(nnr,3), v3x(nnr,2), v1c(nnr,2), v2c(nnr,3), v3c(nnr,2)
!Local variables
INTEGER :: ipol, k, is
REAL(DP) :: grho2(3,nnr)
!MCA: Libxc format
REAL(DP) :: rho_(2,nnr), tau_(2,nnr)
REAL(DP) :: v1x_(2,nnr), v2x_(3,nnr), v3x_(2,nnr), v1c_(2,nnr), v2c_(3,nnr), v3c_(2,nnr)
!_________________________________________________________________________
grho2 = 0.0
!MCA/HK: contracted gradient of density, same format as in libxc
do k=1,nnr
do ipol=1,3
grho2(1,k) = grho2(1,k) + grho(ipol,k,1)**2
grho2(2,k) = grho2(2,k) + grho(ipol,k,1) * grho(ipol,k,2)
grho2(3,k) = grho2(3,k) + grho(ipol,k,2)**2
end do
!MCA: transforming to libxc format (DIRTY HACK)
do is=1,2
rho_(is,k) = rho(k,is)
tau_(is,k) = tau(k,is)
enddo
end do
if (imeta == 5) then
!MCA/HK: using the arrays in libxc format
CALL scancxc_array_spin (nnr, rho_, grho2, tau_, ex, ec, &
& v1x_, v2x_, v3x_, &
& v1c_, v2c_, v3c_ )
do k=1,nnr
!MCA: from libxc to QE format (DIRTY HACK)
do is=1,2
v1x(k,is) = v1x_(is,k)
v2x(k,is) = v2x_(is,k) !MCA/HK: v2x(:,2) contains the cross terms
v3x(k,is) = v3x_(is,k)
v1c(k,is) = v1c_(is,k)
v2c(k,is) = v2c_(is,k) !MCA/HK: same as v2x
v3c(k,is) = v3c_(is,k)
enddo
v2c(k,3) = v2c_(3,k)
v2x(k,3) = v2x_(3,k)
end do
elseif (imeta == 6 ) then ! HK/MCA: SCAN0
CALL scancxc_array (nnr, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c)
if (exx_started) then
ex = (1.0_DP - exx_fraction) * ex
v1x = (1.0_DP - exx_fraction) * v1x
v2x = (1.0_DP - exx_fraction) * v2x
v3x = (1.0_DP - exx_fraction) * v3x
end if
else
CALL errore('v_xc_meta_array','(CP only) array mode only works for SCAN',1)
end if
RETURN
END SUBROUTINE tau_xc_array_spin
!
!
#if defined(__LIBXC)

3
Modules/make.depend
View File

@ -113,6 +113,7 @@ funct.o : kind.o
funct.o : libxc.o
funct.o : xc_gga_drivers.o
funct.o : xc_lda_lsda_drivers.o
funct.o : xc_mgga_drivers.o
funct.o : xc_rVV10.o
funct.o : xc_vdW_DF.o
generate_function.o : ../UtilXlib/mp.o
@ -445,6 +446,8 @@ xc_gga_drivers.o : kind.o
xc_gga_drivers.o : libxc.o
xc_lda_lsda_drivers.o : kind.o
xc_lda_lsda_drivers.o : libxc.o
xc_mgga_drivers.o : kind.o
xc_mgga_drivers.o : libxc.o
xc_rVV10.o : ../FFTXlib/fft_interfaces.o
xc_rVV10.o : ../UtilXlib/mp.o
xc_rVV10.o : cell_base.o

1953
Modules/metagga.f90
View File

File diff suppressed because it is too large Load Diff

420
Modules/xc_mgga_drivers.f90 Normal file
View File

@ -0,0 +1,420 @@
MODULE xc_mgga
!
USE kinds, ONLY: DP
!
IMPLICIT NONE
!
PRIVATE
SAVE
!
! GGA exchange-correlation drivers
PUBLIC :: xc_metagcx
PUBLIC :: tau_xc, tau_xc_spin
PUBLIC :: change_threshold_mgga, select_mgga_functionals
!
PUBLIC :: libxc_switches_mgga
PUBLIC :: imeta_l, imetac_l
PUBLIC :: exx_started_mg, exx_fraction_mg
!
! libxc on/off
INTEGER :: libxc_switches_mgga(2)
!
! indexes defining xc functionals
INTEGER :: imeta_l, imetac_l
!
! input thresholds (default values)
REAL(DP) :: rho_threshold = 1.0E-8_DP
REAL(DP) :: grho2_threshold = 1.0E-12_DP
REAL(DP) :: tau_threshold = 1.0E-8_DP
!
! variables for hybrid exchange
LOGICAL :: exx_started_mg
REAL(DP) :: exx_fraction_mg
!
!
CONTAINS
!
!
!----------------------------------------------------------------------------
!----- Select functionals by the corresponding indexes ----------------------
!----------------------------------------------------------------------------
SUBROUTINE select_mgga_functionals( imeta, imetac, exx_fraction )
!-----------------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: imeta, imetac
REAL(DP), INTENT(IN), OPTIONAL :: exx_fraction
!
! exchange-correlation indexes
imeta_l = imeta
imetac_l = imetac
!
! hybrid exchange vars
exx_started_mg = .FALSE.
exx_fraction_mg = 0._DP
IF ( PRESENT(exx_fraction) ) THEN
exx_started_mg = .TRUE.
exx_fraction_mg = exx_fraction
ENDIF
!
RETURN
!
END SUBROUTINE select_mgga_functionals
!
!
!
!-------------------------------------------------------------------------------------
SUBROUTINE change_threshold_mgga( rho_thr_in, grho2_thr_in, tau_thr_in )
!------------------------------------------------------------------------------------
!! Change rho, grho and tau thresholds.
!
REAL(DP), INTENT(IN) :: rho_thr_in
REAL(DP), INTENT(IN), OPTIONAL :: grho2_thr_in
REAL(DP), INTENT(IN), OPTIONAL :: tau_thr_in
!
rho_threshold = rho_thr_in
IF ( PRESENT(grho2_thr_in) ) grho2_threshold = grho2_thr_in
IF ( PRESENT(tau_thr_in) ) tau_threshold = tau_thr_in
!
RETURN
!
END SUBROUTINE change_threshold_mgga
!
!
!
!----------------------------------------------------------------------------------------
SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!-------------------------------------------------------------------------------------
!! Wrapper routine. Calls metaGGA drivers from internal libraries
!! of q-e or from the external libxc, depending on the input choice.
!
#if defined(__LIBXC)
USE xc_f90_types_m
USE xc_f90_lib_m
#endif
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!! length of the I/O arrays
INTEGER, INTENT(IN) :: ns
!! spin components
INTEGER, INTENT(IN) :: np
!! first dimension of v2c
REAL(DP), INTENT(IN) :: rho(length,ns)
!! the charge density
REAL(DP), INTENT(IN) :: grho(3,length,ns)
!! grho = \nabla rho
REAL(DP), INTENT(IN) :: tau(length,ns)
!! kinetic energy density
REAL(DP), INTENT(OUT) :: ex(length)
!! sx = E_x(rho,grho)
REAL(DP), INTENT(OUT) :: ec(length)
!! sc = E_c(rho,grho)
REAL(DP), INTENT(OUT) :: v1x(length,ns)
!! v1x = D(E_x)/D(rho)
REAL(DP), INTENT(OUT) :: v2x(length,ns)
!! v2x = D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
REAL(DP), INTENT(OUT) :: v3x(length,ns)
!! v3x = D(E_x)/D(tau)
REAL(DP), INTENT(OUT) :: v1c(length,ns)
!! v1c = D(E_c)/D(rho)
REAL(DP), INTENT(OUT) :: v2c(np,length,ns)
!! v2c = D(E_c)/D( D rho/D r_alpha ) / |\nabla rho|
REAL(DP), INTENT(OUT) :: v3c(length,ns)
!! v3c = D(E_c)/D(tau)
!
! ... local variables
!
INTEGER :: is
REAL(DP), ALLOCATABLE :: grho2(:,:)
!
#if defined(__LIBXC)
TYPE(xc_f90_pointer_t) :: xc_func
TYPE(xc_f90_pointer_t) :: xc_info1, xc_info2
!
REAL(DP), ALLOCATABLE :: rho_lxc(:), sigma(:), tau_lxc(:)
REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
REAL(DP), ALLOCATABLE :: vx_rho(:), vx_sigma(:), vx_tau(:)
REAL(DP), ALLOCATABLE :: vc_rho(:), vc_sigma(:), vc_tau(:)
REAL(DP), ALLOCATABLE :: lapl_rho(:), vlapl_rho(:) ! not used in TPSS
!
INTEGER :: k, ipol, pol_unpol
LOGICAL :: POLARIZED
!
POLARIZED = .FALSE.
IF (ns == 2) THEN
POLARIZED = .TRUE.
ENDIF
!
pol_unpol = ns
!
ALLOCATE( rho_lxc(length*ns), sigma(length*ns), tau_lxc(length*ns) )
ALLOCATE( lapl_rho(length*ns) )
!
ALLOCATE( ex_lxc(length) , ec_lxc(length) )
ALLOCATE( vx_rho(length*ns) , vx_sigma(length*ns), vx_tau(length*ns) )
ALLOCATE( vc_rho(length*ns) , vc_sigma(length*ns), vc_tau(length*ns) )
ALLOCATE( vlapl_rho(length*ns) )
!
!
IF ( ns == 1 ) THEN
!
DO k = 1, length
rho_lxc(k) = ABS( rho(k,1) )
IF ( rho_lxc(k) > rho_threshold ) &
sigma(k) = grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2
tau_lxc(k) = tau(k,1)
ENDDO
!
ELSE
!
DO k = 1, length
rho_lxc(2*k-1) = rho(k,1)
rho_lxc(2*k) = rho(k,2)
!
sigma(3*k-2) = grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2
sigma(3*k-1) = grho(1,k,1) * grho(1,k,2) + grho(2,k,1) * grho(2,k,2) + &
grho(3,k,1) * grho(3,k,2)
sigma(3*k) = grho(1,k,2)**2 + grho(2,k,2)**2 + grho(3,k,2)**2
!
tau_lxc(2*k-1) = tau(k,1)
tau_lxc(2*k) = tau(k,2)
ENDDO
!
ENDIF
!
IF (SUM(libxc_switches_mgga(:)) /= 2) THEN
!
ALLOCATE( grho2(length,ns) )
!
DO is = 1, ns
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
IF (ns == 1) THEN
CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, v1x(:,1), &
v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1), v3c(:,1) )
ELSEIF (ns == 2) THEN
CALL tau_xc_spin( length, rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, &
v2c, v3c )
ENDIF
!
DEALLOCATE( grho2 )
!
ENDIF
!
! META EXCHANGE
!
IF (libxc_switches_mgga(1) == 1) THEN
CALL xc_f90_func_init( xc_func, xc_info1, imeta_l, pol_unpol )
CALL xc_f90_mgga_exc_vxc( xc_func, length, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
ex_lxc(1), vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
CALL xc_f90_func_end( xc_func )
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
ex(k) = ex_lxc(k) * rho_lxc(k)
v1x(k,1) = vx_rho(k)
v2x(k,1) = vx_sigma(k) * 2.0_DP
v3x(k,1) = vx_tau(k)
ENDDO
ELSE
DO k = 1, length
ex(k) = ex_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
v1x(k,1) = vx_rho(2*k-1)
v1x(k,2) = vx_rho(2*k)
v2x(k,1) = vx_sigma(3*k-2)*2.d0
v2x(k,2) = vx_sigma(3*k)*2.d0
v3x(k,1) = vx_tau(2*k-1)
v3x(k,2) = vx_tau(2*k)
ENDDO
ENDIF
!
ENDIF
!
! META CORRELATION
!
IF ( libxc_switches_mgga(2) == 1 ) THEN
!
CALL xc_f90_func_init( xc_func, xc_info1, imetac_l, pol_unpol )
CALL xc_f90_mgga_exc_vxc( xc_func, length, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
ec_lxc(1), vc_rho(1), vc_sigma(1), vlapl_rho(1), vc_tau(1) )
CALL xc_f90_func_end( xc_func )
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
ec(k) = ec_lxc(k) * rho_lxc(k) !* SIGN(1.0_DP, rho(k,1))
v1c(k,1) = vc_rho(k)
v2c(1,k,1) = vc_sigma(k) * 2.0_DP
v3c(k,1) = vc_tau(k)
ENDDO
ELSE
DO k = 1, length
ec(k) = ec_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
v1c(k,1) = vc_rho(2*k-1)
v1c(k,2) = vc_rho(2*k)
DO ipol = 1, 3
v2c(ipol,k,1) = vc_sigma(3*k-2)*grho(ipol,k,1)*2.D0 + vc_sigma(3*k-1)*grho(ipol,k,2)
v2c(ipol,k,2) = vc_sigma(3*k) *grho(ipol,k,2)*2.D0 + vc_sigma(3*k-1)*grho(ipol,k,1)
ENDDO
v3c(k,1) = vc_tau(2*k-1)
v3c(k,2) = vc_tau(2*k)
ENDDO
ENDIF
!
ENDIF
!
DEALLOCATE( rho_lxc, sigma, tau_lxc, lapl_rho )
DEALLOCATE( ex_lxc , ec_lxc )
DEALLOCATE( vx_rho , vx_sigma, vx_tau )
DEALLOCATE( vc_rho , vc_sigma, vc_tau, vlapl_rho )
!
#else
!
ALLOCATE( grho2(length,ns) )
!
DO is = 1, ns
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
IF (ns == 1) THEN
!
CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, v1x(:,1), &
v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1), v3c(:,1) )
!
ELSEIF (ns == 2) THEN
!
CALL tau_xc_spin( length, rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, &
v2c, v3c )
!
ENDIF
!
DEALLOCATE( grho2 )
!
#endif
!
RETURN
!
END SUBROUTINE xc_metagcx
!
!
!---------------------------------------------------------------------------------
SUBROUTINE tau_xc( length, rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!-------------------------------------------------------------------------------
! gradient corrections for exchange and correlation - Hartree a.u.
! See comments at the beginning of module for implemented cases
!
! input: rho, grho=|\nabla rho|^2
!
! definition: E_x = \int e_x(rho,grho) dr
!
! output: sx = e_x(rho,grho) = grad corr
! v1x= D(E_x)/D(rho)
! v2x= D(E_x)/D( D rho/D r_alpha ) / |\nabla rho|
! v3x= D(E_x)/D(tau)
!
! sc, v1c, v2c as above for correlation
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!
INTEGER :: k
REAL(DP) :: arho
REAL(DP), DIMENSION(length) :: rho, grho2, tau, &
ex, ec, v1x, v2x, v3x, v1c, v2c, v3c
!
v1x=0.d0 ; v2x=0.d0 ; v3x=0.d0 ; ex=0.d0
v1c=0.d0 ; v2c=0.d0 ; v3c=0.d0 ; ec=0.d0
!
DO k = 1, length
!
arho = ABS(rho(k))
!
IF ( (arho<=rho_threshold).OR.(grho2(k)<=grho2_threshold).OR.(ABS(tau(k))<=rho_threshold) ) CYCLE
!
SELECT CASE( imeta_l )
CASE( 1 )
CALL tpsscxc( arho, grho2(k), tau(k), ex(k), ec(k), v1x(k), v2x(k), v3x(k), v1c(k), v2c(k), v3c(k) )
CASE( 2 )
CALL m06lxc( arho, grho2(k), tau(k), ex(k), ec(k), v1x(k), v2x(k), v3x(k), v1c(k), v2c(k), v3c(k) )
CASE( 4 )
! do nothing
CASE DEFAULT
CALL errore( 'tau_xc', 'wrong igcx and/or igcc', 1 )
END SELECT
!
ENDDO
!
RETURN
!
END SUBROUTINE tau_xc
!
!
!----------------------------------------------------------------------------------------
SUBROUTINE tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!------------------------------------------------------------------------------------
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN) :: rho(length,2), tau(length,2)
REAL(DP), INTENT(IN) :: grho(3,length,2)
!
REAL(DP), INTENT(OUT) :: ex(length), ec(length), v1x(length,2), v2x(length,2), &
v3x(length,2), v1c(length,2), v3c(length,2)
REAL(DP), INTENT(OUT) :: v2c(3,length,2)
!
! ... local variables
!
INTEGER :: k, ipol
REAL(DP) :: rh, zeta, atau, grho2(2), ggrho2
!
ex=0.0_DP ; v1x=0.0_DP ; v2x=0.0_DP ; v3x=0.0_DP
ec=0.0_DP ; v1c=0.0_DP ; v2c=0.0_DP ; v3c=0.0_DP
!
! FIXME: for SCAN, this will be calculated later
!
DO k = 1, length
!
rh = rho(k,1) + rho(k,2)
atau = tau(k,1) + tau(k,2) ! KE-density in Hartree
grho2(1) = SUM( grho(:,k,1)**2 )
grho2(2) = SUM( grho(:,k,2)**2 )
ggrho2 = grho2(1) + grho2(2)
!
IF ((rh <= rho_threshold) .OR. (ggrho2 <= grho2_threshold) .OR. (ABS(atau) <= tau_threshold)) CYCLE
!
SELECT CASE( imeta_l )
CASE( 1 )
!
CALL tpsscx_spin( rho(k,1), rho(k,2), grho2(1), grho2(2), tau(k,1), &
tau(k,2), ex(k), v1x(k,1), v1x(k,2), v2x(k,1), v2x(k,2), v3x(k,1), v3x(k,2) )
!
zeta = (rho(k,1) - rho(k,2)) / rh
!
CALL tpsscc_spin( rh, zeta, grho(:,k,1), grho(:,k,2), atau, ec(k), &
v1c(k,1), v1c(k,2), v2c(:,k,1), v2c(:,k,2), v3c(k,1), v3c(k,2) )
!
CASE( 2 )
!
CALL m06lxc_spin( rho(k,1), rho(k,2), grho(:,k,1), grho(:,k,2), tau(k,1), tau(k,2), ex(k), ec(k), &
v1x(k,1), v1x(k,2), v2x(k,1), v2x(k,2), v3x(k,1), v3x(k,2), &
v1c(k,1), v1c(k,2), v2c(:,k,1), v2c(:,k,2), v3c(k,1), v3c(k,2) )
!
CASE DEFAULT
!
CALL errore( 'tau_xc_spin', 'This case not implemented', imeta_l )
!
END SELECT
!
ENDDO
!
RETURN
!
END SUBROUTINE tau_xc_spin
!
!
END MODULE xc_mgga

388
PP/src/benchmark_libxc.f90
View File

@ -5,8 +5,17 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!-----------------------------------------------------------------------------------
PROGRAM benchmark_libxc
!--------------------------------------------------------------------------------
!! This program compares the output results (energies and potentials) from the libxc
!! routines with the ones from q-e xc internal library.
!! Available options:
!! * full LDA ;
!! * derivative of LDA pot. (dmxc) ;
!! * full GGA ;
!! * derivative of GGA pot. (dgcxc, the polarized case is not yet complete) ;
!! * full metaGGA .
!
!------------------------------------------------------------------------------------!
! To be run on a single processor
@ -19,6 +28,7 @@ PROGRAM benchmark_libxc
!
USE xc_lda_lsda
USE xc_gga
USE xc_mgga
!
IMPLICIT NONE
!
@ -26,10 +36,10 @@ PROGRAM benchmark_libxc
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(14,200)
INTEGER, PARAMETER :: nnr = 6
CHARACTER(LEN=120) :: aprx, e_q, f_q
INTEGER :: ii, ns, np, quit, i_sub, family
INTEGER :: ii, ns, np, ipol, quit, i_sub, family
REAL(DP) :: exx_frctn
LOGICAL :: LDA, GGA, POLARIZED, ENERGY_ONLY, DF_OK
REAL(DP), PARAMETER :: null = 0.0_DP, pi34 = 0.6203504908994_DP
LOGICAL :: LDA, GGA, MGGA, POLARIZED, ENERGY_ONLY, DF_OK
REAL(DP), PARAMETER :: null=0.0_DP, pi34=0.6203504908994_DP
!
!----------QE vars --------------------------
INTEGER :: iexch_qe, icorr_qe
@ -37,11 +47,13 @@ PROGRAM benchmark_libxc
REAL(DP), ALLOCATABLE :: rho_qe(:,:)
REAL(DP), ALLOCATABLE :: rho_tot(:), zeta(:)
REAL(DP), ALLOCATABLE :: grho(:,:,:), grho_ud(:), grho2(:,:), grh2(:)
REAL(DP), ALLOCATABLE :: tau_qe(:,:)
REAL(DP), ALLOCATABLE :: ex_qe(:), ec_qe(:)
REAL(DP), ALLOCATABLE :: vx_qe(:,:), vc_qe(:,:)
REAL(DP), ALLOCATABLE :: dmuxc(:,:,:)
REAL(DP), ALLOCATABLE :: v1x(:,:), v2x(:,:)
REAL(DP), ALLOCATABLE :: v1x(:,:), v2x(:,:), v3x(:,:)
REAL(DP), ALLOCATABLE :: v1c(:,:), v2c(:,:), v2c_ud(:)
REAL(DP), ALLOCATABLE :: v2cm(:,:,:), v3c(:,:)
REAL(DP), ALLOCATABLE :: vrrx(:,:), vsrx(:,:), vssx(:,:)
REAL(DP), ALLOCATABLE :: vrrc(:,:), vsrc(:,:), vssc(:), vrzc(:,:)
!
@ -52,21 +64,24 @@ PROGRAM benchmark_libxc
INTEGER :: iexch_lxc, icorr_lxc
INTEGER :: pol_unpol
REAL(DP), ALLOCATABLE :: rho_lxc(:)
REAL(DP), ALLOCATABLE :: tau_lxc(:), lapl_rho(:)
REAL(DP), ALLOCATABLE :: sigma(:)
REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
REAL(DP), ALLOCATABLE :: vx_lxc(:), vc_lxc(:)
REAL(DP), ALLOCATABLE :: dex_lxc(:), dcr_lxc(:), df_lxc(:)
REAL(DP), ALLOCATABLE :: vx_rho(:), vc_rho(:)
REAL(DP), ALLOCATABLE :: vx_sigma(:), vc_sigma(:)
REAL(DP), ALLOCATABLE :: vx_sigma(:), vc_sigma(:), v2c_lxc(:,:,:)
REAL(DP), ALLOCATABLE :: vx_lxc2(:), vc_lxc2(:)
REAL(DP), ALLOCATABLE :: ex_lxc2(:), ec_lxc2(:)
REAL(DP), ALLOCATABLE :: v2rho2_x(:), v2rhosigma_x(:), v2sigma2_x(:)
REAL(DP), ALLOCATABLE :: v2rho2_c(:), v2rhosigma_c(:), v2sigma2_c(:)
REAL(DP), ALLOCATABLE :: vx_tau(:), vc_tau(:), vlapl_rho(:)
!
!
! *******************************************************************************
! *-----------------------------------------------------------------------------*
! * libxc funct. indexes: http://bigdft.org/Wiki/index.php?title=XC_codes *
! * (or use function: 'xc_functional_get_number()' ) *
! * qe " " : see comments in Modules/funct.f90 *
! *-----------------------------------------------------------------------------*
! * *
@ -93,9 +108,6 @@ PROGRAM benchmark_libxc
! libxc: ec = ec_glyp (icorr=131) / vc = vc_glyp (131)
! ... same for polarized case
!
! - pbe => differences of the order of 1-5% only in the polarized case by
! adding thw pw lda part (in qe)
!
!
PRINT *, CHAR(10)//" --- BENCHMARK TEST BETWEEN QE AND LIBXC ---"//CHAR(10)//" "
!
@ -104,7 +116,7 @@ PROGRAM benchmark_libxc
DF_OK = .FALSE.
IF ( TRIM(f_q) == 'y' ) DF_OK = .TRUE.
IF ( TRIM(f_q) /= 'y' .AND. TRIM(f_q) /= 'n' ) THEN
PRINT *, CHAR(10)//"ERROR: it is yes (y) or no (n)"//CHAR(10)
PRINT *, CHAR(10)//"Wrong answer"//CHAR(10)
GO TO 10
ENDIF
!
@ -115,15 +127,15 @@ PROGRAM benchmark_libxc
ENERGY_ONLY = .FALSE.
IF ( TRIM(e_q) == 'y' ) ENERGY_ONLY = .TRUE.
IF ( TRIM(e_q) /= 'y' .AND. TRIM(e_q) /= 'n' ) THEN
PRINT *, CHAR(10)//"ERROR: it is yes (y) or no (n)"//CHAR(10)
PRINT *, CHAR(10)//"Wrong answer"//CHAR(10)
GO TO 10
ENDIF
ENDIF
!
WRITE (*,'(/,1x,a)', ADVANCE='no') "lda or gga ? "
WRITE (*,'(/,1x,a)', ADVANCE='no') "lda or gga or mgga ? "
READ(*,*) aprx
IF ( TRIM(aprx) /= 'lda' .AND. TRIM(aprx) /= 'gga' ) THEN
PRINT *, CHAR(10)//"ERROR: you can only choose lda or gga"//CHAR(10)
IF ( TRIM(aprx) /= 'lda' .AND. TRIM(aprx) /= 'gga' .AND. TRIM(aprx) /= 'mgga' ) THEN
PRINT *, CHAR(10)//"ERROR: you can only choose among lda, gga and mgga"//CHAR(10)
GO TO 10
ENDIF
WRITE (*,'(/,1x,a)', ADVANCE='no') "Polarization switch (1 unpolarized, &
@ -146,13 +158,21 @@ PROGRAM benchmark_libxc
ENDIF
!
!
IF ( TRIM(aprx) == 'lda' ) THEN
LDA = .TRUE.
GGA = .FALSE.
ELSE
LDA = .FALSE.
GGA = .TRUE.
ENDIF
SELECT CASE( TRIM(aprx) )
CASE( 'lda' )
LDA = .TRUE.
GGA = .FALSE.
MGGA = .FALSE.
CASE( 'gga' )
LDA = .FALSE.
GGA = .TRUE.
MGGA = .FALSE.
CASE( 'mgga' )
LDA = .FALSE.
GGA = .FALSE.
MGGA = .TRUE.
END SELECT
!
!
POLARIZED = .FALSE.
IF (ns == 2) THEN
@ -194,13 +214,20 @@ PROGRAM benchmark_libxc
IF ( .NOT.POLARIZED ) ALLOCATE( dmuxc(nnr,1,1) )
IF ( POLARIZED ) ALLOCATE( dmuxc(nnr,2,2) )
ENDIF
ELSEIF ( GGA ) THEN
ALLOCATE( grho(nnr,3,ns), grho2(nnr,ns), grho_ud(nnr), grh2(nnr) )
ELSEIF ( GGA .OR. MGGA ) THEN
ALLOCATE( grho(nnr,3,ns), grho_ud(nnr), grho2(nnr,ns) )
ALLOCATE( v1x(nnr,ns), v2x(nnr,ns) )
ALLOCATE( v1c(nnr,ns), v2c(nnr,ns), v2c_ud(nnr) )
IF (DF_OK) THEN
ALLOCATE( vrrx(nnr,ns), vsrx(nnr,ns), vssx(nnr,ns) )
ALLOCATE( vrrc(nnr,ns), vsrc(nnr,ns), vssc(nnr), vrzc(nnr,ns) )
ALLOCATE( v1c(nnr,ns) )
IF ( GGA ) THEN
ALLOCATE( grh2(nnr) )
ALLOCATE( v2c(nnr,ns), v2c_ud(nnr) )
IF ( DF_OK ) THEN
ALLOCATE( vrrx(nnr,ns), vsrx(nnr,ns), vssx(nnr,ns) )
ALLOCATE( vrrc(nnr,ns), vsrc(nnr,ns), vssc(nnr), vrzc(nnr,ns) )
ENDIF
ELSEIF ( MGGA ) THEN
ALLOCATE( v2cm(np,nnr,ns), tau_qe(nnr,ns) )
ALLOCATE( v3x(nnr,ns), v3c(nnr,ns) )
ENDIF
ENDIF
!
@ -214,22 +241,30 @@ PROGRAM benchmark_libxc
IF ( .NOT.POLARIZED ) ALLOCATE( dex_lxc(nnr), dcr_lxc(nnr), df_lxc(nnr) )
IF ( POLARIZED ) ALLOCATE( dex_lxc(nnr*3), dcr_lxc(nnr*3), df_lxc(nnr*3) )
ENDIF
ELSEIF ( GGA ) THEN
ELSEIF ( GGA .OR. MGGA ) THEN
ALLOCATE( sigma(nnr*np) )
ALLOCATE( vx_rho(nnr*ns), vx_sigma(nnr*np) )
ALLOCATE( vc_rho(nnr*ns), vc_sigma(nnr*np) )
ALLOCATE( vx_lxc2(nnr*ns), vc_lxc2(nnr*ns) )
ALLOCATE( ex_lxc2(nnr), ec_lxc2(nnr) )
IF ( DF_OK ) THEN
IF ( .NOT.POLARIZED ) THEN
ALLOCATE( dex_lxc(nnr), dcr_lxc(nnr), df_lxc(nnr) )
ALLOCATE( v2rho2_x(nnr), v2rhosigma_x(nnr), v2sigma2_x(nnr) )
ALLOCATE( v2rho2_c(nnr), v2rhosigma_c(nnr), v2sigma2_c(nnr) )
ELSEIF ( POLARIZED ) THEN
ALLOCATE( dex_lxc(nnr*3), dcr_lxc(nnr*3), df_lxc(nnr*3) )
ALLOCATE( v2rho2_x(3*nnr), v2rhosigma_x(6*nnr), v2sigma2_x(6*nnr) )
ALLOCATE( v2rho2_c(3*nnr), v2rhosigma_c(6*nnr), v2sigma2_c(6*nnr) )
!
IF ( GGA ) THEN
ALLOCATE( vx_lxc2(nnr*ns), vc_lxc2(nnr*ns) )
ALLOCATE( ex_lxc2(nnr), ec_lxc2(nnr) )
IF ( DF_OK ) THEN
IF ( .NOT.POLARIZED ) THEN
ALLOCATE( dex_lxc(nnr), dcr_lxc(nnr), df_lxc(nnr) )
ALLOCATE( v2rho2_x(nnr), v2rhosigma_x(nnr), v2sigma2_x(nnr) )
ALLOCATE( v2rho2_c(nnr), v2rhosigma_c(nnr), v2sigma2_c(nnr) )
ELSEIF ( POLARIZED ) THEN
ALLOCATE( dex_lxc(nnr*3), dcr_lxc(nnr*3), df_lxc(nnr*3) )
ALLOCATE( v2rho2_x(3*nnr), v2rhosigma_x(6*nnr), v2sigma2_x(6*nnr) )
ALLOCATE( v2rho2_c(3*nnr), v2rhosigma_c(6*nnr), v2sigma2_c(6*nnr) )
ENDIF
ENDIF
ELSEIF ( MGGA ) THEN
ALLOCATE( tau_lxc(nnr*ns), lapl_rho(nnr*ns*np) )
ALLOCATE( vx_tau(nnr*ns), vc_tau(nnr*ns) )
ALLOCATE( vlapl_rho(nnr*ns*np) )
ALLOCATE( v2c_lxc(np,nnr,ns) )
ENDIF
ENDIF
!
@ -243,16 +278,17 @@ PROGRAM benchmark_libxc
rho_tot = 0.0_DP
zeta = 0.0_DP
ENDIF
IF ( GGA ) THEN
IF ( GGA .OR. MGGA ) THEN
grho = 0.0_DP
grho2 = 0.0_DP
grho_ud = 0.0_DP
IF (GGA) grho_ud = 0.0_DP
ENDIF
!
! ... libcx
!
rho_lxc = 0.0_DP
IF ( GGA ) sigma = 0.0_DP
IF ( GGA .OR. MGGA ) sigma = 0.0_DP
lapl_rho = 0.0_DP
!
! -------- Setting up an arbitrary input for both qe and libxc -----
!
@ -262,36 +298,41 @@ PROGRAM benchmark_libxc
!
rho_qe(ii,1) = DBLE(ii)/DBLE(nnr+2)
!
IF ( GGA ) THEN
IF ( GGA .OR. MGGA ) THEN
grho(ii,1,1) = ABS( 0.05_DP + 0.8_DP*SIN(DBLE(ii)) )
grho(ii,2,1) = ABS( 0.05_DP + 0.7_DP*SIN(DBLE(ii)) )
grho(ii,3,1) = ABS( 0.05_DP + 0.6_DP*SIN(DBLE(ii)) )
ENDIF
grho2(ii,1) = grho(ii,1,1)**2 + grho(ii,2,1)**2 + grho(ii,3,1)**2
!
IF ( MGGA ) tau_qe(ii,1) = ABS( 0.05_DP + 0.8_DP*SIN(DBLE(ii)) )*0.5d0
!
IF ( POLARIZED ) THEN
!
rho_qe(ii,2) = (1.0_DP - rho_qe(ii,1))*0.7_DP
rho_tot(ii) = rho_qe(ii,1) + rho_qe(ii,2)
zeta(ii) = (rho_qe(ii,1) - rho_qe(ii,2)) / rho_tot(ii)
!
IF ( GGA ) THEN
IF ( GGA .OR. MGGA ) THEN
grho(ii,1,2) = ABS( (1.0_DP - grho(ii,1,1))*0.7_DP )
grho(ii,2,2) = ABS( (1.0_DP - grho(ii,2,1))*0.6_DP )
grho(ii,3,2) = ABS( (1.0_DP - grho(ii,3,1))*0.5_DP )
!
grh2(ii)= ( grho(ii,1,1) + grho(ii,1,2) )**2 + &
( grho(ii,2,1) + grho(ii,2,2) )**2 + &
( grho(ii,3,1) + grho(ii,3,2) )**2
!
grho2(ii,2) = ( grho(ii,1,2)**2 + grho(ii,2,2)**2 + grho(ii,3,2)**2 )
!
grho_ud(ii) = grho(ii,1,1) * grho(ii,1,2) + &
grho(ii,2,1) * grho(ii,2,2) + &
grho(ii,3,1) * grho(ii,3,2)
IF (GGA) THEN
grh2(ii) = ( grho(ii,1,1) + grho(ii,1,2) )**2 + &
( grho(ii,2,1) + grho(ii,2,2) )**2 + &
( grho(ii,3,1) + grho(ii,3,2) )**2
ENDIF
!
ENDIF
!
IF ( MGGA ) tau_qe(ii,2) = ABS( 0.05_DP + 0.8_DP*SIN(DBLE(ii)) )*0.2d0
!
ENDIF
!
ENDDO
@ -304,17 +345,23 @@ PROGRAM benchmark_libxc
!
rho_lxc(ii) = rho_qe(ii,1)
!
IF ( GGA ) sigma(ii) = grho2(ii,1)
IF ( GGA .OR. MGGA ) sigma(ii) = grho2(ii,1)
IF ( MGGA ) tau_lxc(ii) = tau_qe(ii,1)
!
ELSE
!
rho_lxc(2*ii-1) = rho_qe(ii,1)
rho_lxc(2*ii) = rho_qe(ii,2)
!
IF ( GGA ) THEN
sigma(3*ii-2) = grho2(ii,1)
sigma(3*ii-1) = grho_ud(ii)
sigma(3*ii) = grho2(ii,2)
IF ( GGA .OR. MGGA ) THEN
sigma(3*ii-2) = grho2(ii,1)
sigma(3*ii-1) = grho_ud(ii)
sigma(3*ii) = grho2(ii,2)
ENDIF
!
IF ( MGGA ) THEN
tau_lxc(2*ii-1) = tau_qe(ii,1)
tau_lxc(2*ii) = tau_qe(ii,2)
ENDIF
!
ENDIF
@ -411,7 +458,6 @@ PROGRAM benchmark_libxc
IF (icorr_lxc /= 131 ) then ! .AND. .NOT.(icorr_lxc == 130 .AND. POLARIZED)) THEN
! remove LDA correlation for compatibility with QE
i_sub=12 !(pw)
!IF (icorr_lxc == 132) i_sub=9 !(pz)
CALL xc_f90_func_init( xc_func, xc_info4, i_sub, pol_unpol )
CALL xc_f90_lda_exc_vxc( xc_func, nnr, rho_lxc(1), ec_lxc2(1), vc_lxc2(1) )
!
@ -448,8 +494,6 @@ PROGRAM benchmark_libxc
CALL select_gga_functionals( iexch_qe, icorr_qe, exx_fraction=exx_frctn )
!
IF ( DF_OK ) THEN
!
!CALL select_lda_functionals( iexch_qe, icorr_qe )
!
IF (.NOT. POLARIZED) THEN
CALL dgcxc( nnr, rho_qe(:,1), grho2(:,1), vrrx(:,1), vsrx(:,1), vssx(:,1), &
@ -482,15 +526,6 @@ PROGRAM benchmark_libxc
v2c(:,2) = v2c(:,1)
v2c_ud(:) = v2c(:,1)
!
!
!IF (icorr_qe == 4) THEN
! rs(:) = pi34 / rho_tot(:)**(1.d0/3.d0) ! .... trovare le combinazioni e mettere ordine
! CALL pw_spin( nnr, rs, zeta, ec_qe2, vc_qe2 )
! ec_qe(:) = ec_qe(:) + ec_qe2(:)*rho_tot(:)
! v1c(:,1) = v1c(:,1) + vc_qe2(:,1)
! v1c(:,2) = v1c(:,2) + vc_qe2(:,2)
!ENDIF
!
ELSE
CALL gcc_spin_more( nnr, rho_qe, grho2, grho_ud, ec_qe, v1c, v2c, v2c_ud )
CALL lsd_lyp( nnr, rho_tot, zeta, ec_qe2, vc_qe2 )
@ -504,12 +539,61 @@ PROGRAM benchmark_libxc
!
ENDIF
!
!ec_qe=ec_qe*2.d0
!v1c = v1c * 2.d0
ELSEIF ( MGGA ) THEN
!
ENDIF
!------ LIBXC ------
!
! exch
CALL xc_f90_func_init( xc_func, xc_info1, iexch_lxc, pol_unpol )
CALL xc_f90_mgga_exc_vxc( xc_func, nnr, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
ex_lxc(1), vx_rho(1), vx_sigma(1), vlapl_rho(1), vx_tau(1) )
CALL xc_f90_func_end( xc_func )
!
IF (.NOT. POLARIZED) THEN
ex_lxc = ex_lxc * rho_qe(:,1)
ELSE
ex_lxc = ex_lxc * rho_tot
ENDIF
vx_sigma = vx_sigma * 2.0_DP
!
! corr
CALL xc_f90_func_init( xc_func, xc_info2, icorr_lxc, pol_unpol )
CALL xc_f90_mgga_exc_vxc( xc_func, nnr, rho_lxc(1), sigma(1), lapl_rho(1), tau_lxc(1), &
ec_lxc(1), vc_rho(1), vc_sigma(1), vlapl_rho(1), vc_tau(1) )
CALL xc_f90_func_end( xc_func )
!
IF (.NOT. POLARIZED) THEN
ec_lxc = ec_lxc * rho_qe(:,1)
ELSE
ec_lxc = ec_lxc * rho_tot
ENDIF
!
IF (.NOT. POLARIZED) THEN
vc_sigma = vc_sigma * 2.0_DP
ELSE
DO ii = 1, nnr
DO ipol = 1, 3
v2c_lxc(ipol,ii,1) = vc_sigma(3*ii-2)*grho(ii,ipol,1) * 2.D0 + vc_sigma(3*ii-1)*grho(ii,ipol,2)
v2c_lxc(ipol,ii,2) = vc_sigma(3*ii) *grho(ii,ipol,2) * 2.D0 + vc_sigma(3*ii-1)*grho(ii,ipol,1)
ENDDO
ENDDO
ENDIF
!
!----- QE ----------
!
CALL select_mgga_functionals( iexch_qe, icorr_qe ) ! ... icorr_qe not used
!
IF ( .NOT. POLARIZED ) THEN
CALL tau_xc( nnr, rho_qe(:,1), grho2(:,1), tau_qe(:,1), ex_qe, ec_qe, v1x(:,1), &
v2x(:,1), v3x(:,1), v1c(:,1), v2cm(1,:,1), v3c(:,1) )
ELSE
CALL tau_xc_spin( nnr, rho_qe, grho2, tau_qe, ex_qe, ec_qe, v1x, v2x, v3x, v1c, &
v2cm, v3c )
ENDIF
!
ENDIF
!
!------------------
!--
!
CALL xc_f90_info_name( xc_info1, name1 )
CALL xc_f90_info_name( xc_info2, name2 )
@ -520,6 +604,7 @@ PROGRAM benchmark_libxc
PRINT *, "Correlation: ", TRIM(name2)
PRINT *, " "
!
!
IF ( LDA ) THEN
!
DO ii = 1, nnr, nnr-1
@ -574,14 +659,14 @@ PROGRAM benchmark_libxc
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) dmuxc(ii,1,1)
WRITE (*,201) df_lxc(ii)
PRINT *, " --- "
WRITE (*,301) dmuxc(ii,1,1)-df_lxc(ii)
ELSE
WRITE (*,104) dmuxc(ii,1,1), dmuxc(ii,2,1), dmuxc(ii,2,2), dmuxc(ii,1,2)
WRITE (*,203) df_lxc(3*ii-2), df_lxc(3*ii-1), df_lxc(3*ii)
PRINT *, " --- "
WRITE (*,303) dmuxc(ii,1,1)-df_lxc(3*ii-2), dmuxc(ii,2,1)-df_lxc(3*ii-1), &
dmuxc(ii,2,2)-df_lxc(3*ii)
PRINT *, " --- "
WRITE (*,301) dmuxc(ii,1,1)-df_lxc(ii)
ELSE
WRITE (*,104) dmuxc(ii,1,1), dmuxc(ii,2,1), dmuxc(ii,2,2), dmuxc(ii,1,2)
WRITE (*,203) df_lxc(3*ii-2), df_lxc(3*ii-1), df_lxc(3*ii)
PRINT *, " --- "
WRITE (*,303) dmuxc(ii,1,1)-df_lxc(3*ii-2), dmuxc(ii,2,1)-df_lxc(3*ii-1), &
dmuxc(ii,2,2)-df_lxc(3*ii)
ENDIF
ENDIF
!
@ -623,7 +708,7 @@ PROGRAM benchmark_libxc
WRITE (*,202) vx_rho(2*ii-1), vx_rho(2*ii)
PRINT *, " --- "
WRITE (*,302) v1x(ii,1)-vx_rho(2*ii-1), v1x(ii,2)-vx_rho(2*ii)
ENDIF
ENDIF
!
PRINT *, " "
PRINT *, "=== Exchange potential vsigma ==="
@ -738,6 +823,121 @@ PROGRAM benchmark_libxc
!
ENDDO
!
ELSEIF ( MGGA ) THEN
!
DO ii = 1, nnr !, nnr-1
WRITE(*,*) ' '
WRITE(*,*) ' '
WRITE(*,909) ii, nnr
IF (.NOT. POLARIZED ) THEN
WRITE (*,401) rho_qe(ii,1)
WRITE (*,501) grho2(ii,1)
WRITE (*,601) tau_qe(ii,1)
ELSE
WRITE (*,402) rho_qe(ii,1), rho_qe(ii,2)
WRITE (*,502) grho2(ii,1), grho2(ii,2)
WRITE (*,602) tau_qe(ii,1), tau_qe(ii,2)
ENDIF
!
PRINT *, " "
PRINT *, "=== Exchange and correlation energies: ==="
WRITE (*,102) ex_qe(ii), ec_qe(ii)
WRITE (*,202) ex_lxc(ii), ec_lxc(ii)
PRINT *, " --- "
WRITE (*,302) ex_qe(ii)-ex_lxc(ii), ec_qe(ii)-ec_lxc(ii)
!WRITE (*,302) ex_qe(ii)/ex_lxc(ii), ec_qe(ii)/ec_lxc(ii)
!
IF (.NOT. ENERGY_ONLY) THEN
!
PRINT *, " "
PRINT *, "=== Exchange potential vrho ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) v1x(ii,1)
WRITE (*,201) vx_rho(ii)
PRINT *, " --- "
WRITE (*,301) v1x(ii,1)-vx_rho(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) v1x(ii,1), v1x(ii,2)
WRITE (*,202) vx_rho(2*ii-1), vx_rho(2*ii)
PRINT *, " --- "
WRITE (*,302) v1x(ii,1)-vx_rho(2*ii-1), v1x(ii,2)-vx_rho(2*ii)
ENDIF
!
PRINT *, " "
PRINT *, "=== Exchange potential vsigma ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) v2x(ii,1)
WRITE (*,201) vx_sigma(ii)
PRINT *, " --- "
WRITE (*,301) v2x(ii,1)-vx_sigma(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,103) v2x(ii,1), null, v2x(ii,2)
WRITE (*,203) vx_sigma(3*ii-2), vx_sigma(3*ii-1), vx_sigma(3*ii)
PRINT *, " --- "
WRITE (*,303) v2x(ii,1)-vx_sigma(3*ii-2), null-vx_sigma(3*ii-1), &
v2x(ii,2)-vx_sigma(3*ii)
ENDIF
!
PRINT *, " "
PRINT *, "=== Exchange potential vtau ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) v3x(ii,1)
WRITE (*,201) vx_tau(ii)
PRINT *, " --- "
WRITE (*,301) v3x(ii,1)-vx_tau(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) v3x(ii,1), v3x(ii,2)
WRITE (*,202) vx_tau(2*ii-1), vx_tau(2*ii)
PRINT *, " --- "
WRITE (*,303) v3x(ii,1)-vx_tau(2*ii-1), v3x(ii,2)-vx_tau(2*ii)
ENDIF
!
PRINT *, " "
PRINT *, "=== Correlation potential vrho ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) v1c(ii,1)
WRITE (*,201) vc_rho(ii)
PRINT *, " --- "
WRITE (*,301) v1c(ii,1)-vc_rho(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) v1c(ii,1), v1c(ii,2)
WRITE (*,202) vc_rho(2*ii-1), vc_rho(2*ii)
PRINT *, " --- "
WRITE (*,302) v1c(ii,1)-vc_rho(2*ii-1), v1c(ii,2)-vc_rho(2*ii)
ENDIF
!
PRINT *, " "
PRINT *, "=== Correlation potential vsigma ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) v2cm(1,ii,1)
WRITE (*,201) vc_sigma(ii)
PRINT *, " --- "
WRITE (*,301) v2cm(1,ii,1)-vc_sigma(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) v2cm(1,ii,1), v2cm(1,ii,2)
WRITE (*,203) v2c_lxc(1,ii,1), v2c_lxc(1,ii,2)
PRINT *, " --- "
WRITE (*,303) v2cm(1,ii,1)-v2c_lxc(1,ii,1), v2cm(1,ii,2)-v2c_lxc(1,ii,2)
ENDIF
!
PRINT *, " "
PRINT *, "=== Correlation potential vtau ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) v3c(ii,1)
WRITE (*,201) vc_tau(ii)
PRINT *, " --- "
WRITE (*,301) v3c(ii,1)-vc_tau(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) v3c(ii,1), v3c(ii,2)
WRITE (*,202) vc_tau(2*ii-1), vc_tau(2*ii)
PRINT *, " --- "
WRITE (*,302) v3c(ii,1)-vc_tau(2*ii-1), v3c(ii,2)-vc_tau(2*ii)
ENDIF
!
ENDIF
!
ENDDO
!
ENDIF
!
101 FORMAT('qe: ',3x,F17.14)
@ -757,10 +957,15 @@ PROGRAM benchmark_libxc
402 FORMAT('rho(up,down): ',F17.14,4x,F17.14)
!
501 FORMAT('grho2: ',F17.14)
502 FORMAT('grho2(uu,dd): ',F17.14,4x,F17.14)
503 FORMAT('grho2(uu,ud,dd): ',F17.14,4x,F17.14,4x,F17.14)
!
601 FORMAT('tau: ',F17.14)
602 FORMAT('tau(up,down): ',F17.14,4x,F17.14)
!
909 FORMAT('grid-point: ',I4,' of',I4)
!
! -- qe
DEALLOCATE( rho_qe )
IF ( POLARIZED ) DEALLOCATE( rho_tot, zeta )
DEALLOCATE( ex_qe, ec_qe )
@ -769,16 +974,23 @@ PROGRAM benchmark_libxc
IF ( DF_OK ) THEN
DEALLOCATE( dmuxc )
ENDIF
ELSEIF ( GGA ) THEN
DEALLOCATE( grho, grho2, grho_ud, grh2 )
DEALLOCATE( v1x, v2x )
DEALLOCATE( v1c, v2c, v2c_ud )
IF (DF_OK) THEN
DEALLOCATE( vrrx, vsrx, vssx )
DEALLOCATE( vrrc, vsrc, vssc, vrzc )
ELSEIF ( GGA .OR. MGGA ) THEN
DEALLOCATE( grho, grho_ud, grho2 )
DEALLOCATE( v1x, v2x, v1c )
IF ( GGA ) THEN
DEALLOCATE( grh2 )
DEALLOCATE( v2c, v2c_ud )
IF (DF_OK) THEN
DEALLOCATE( vrrx, vsrx, vssx )
DEALLOCATE( vrrc, vsrc, vssc, vrzc )
ENDIF
ELSEIF ( MGGA ) THEN
DEALLOCATE( tau_qe )
DEALLOCATE( v2cm, v3x, v3c )
ENDIF
ENDIF
!
! -- libxc
DEALLOCATE( rho_lxc )
DEALLOCATE( ex_lxc, ec_lxc )
IF ( LDA ) THEN
@ -796,6 +1008,12 @@ PROGRAM benchmark_libxc
DEALLOCATE( v2rho2_x, v2rhosigma_x, v2sigma2_x )
DEALLOCATE( v2rho2_c, v2rhosigma_c, v2sigma2_c )
ENDIF
ELSEIF ( MGGA ) THEN
DEALLOCATE( sigma, tau_lxc, lapl_rho )
DEALLOCATE( vx_rho, vx_sigma, vx_tau )
DEALLOCATE( vc_rho, vc_sigma, vc_tau )
DEALLOCATE( vlapl_rho )
DEALLOCATE( v2c_lxc )
ENDIF
!
PRINT *, " "

1
PP/src/make.depend
View File

@ -91,6 +91,7 @@ bands.o : ../../UtilXlib/mp.o
benchmark_libxc.o : ../../Modules/libxc.o
benchmark_libxc.o : ../../Modules/xc_gga_drivers.o
benchmark_libxc.o : ../../Modules/xc_lda_lsda_drivers.o
benchmark_libxc.o : ../../Modules/xc_mgga_drivers.o
chdens_bspline.o : ../../Modules/bspline.o
chdens_bspline.o : ../../Modules/cell_base.o
chdens_bspline.o : ../../Modules/fft_base.o

2
PW/src/make.depend
View File

@ -1880,6 +1880,7 @@ stres_gradcorr.o : ../../Modules/funct.o
stres_gradcorr.o : ../../Modules/kind.o
stres_gradcorr.o : ../../Modules/mp_bands.o
stres_gradcorr.o : ../../Modules/xc_gga_drivers.o
stres_gradcorr.o : ../../Modules/xc_mgga_drivers.o
stres_gradcorr.o : ../../UtilXlib/mp.o
stres_har.o : ../../FFTXlib/fft_interfaces.o
stres_har.o : ../../Modules/cell_base.o
@ -2167,6 +2168,7 @@ v_of_rho.o : ../../Modules/noncol.o
v_of_rho.o : ../../Modules/recvec.o
v_of_rho.o : ../../Modules/tsvdw.o
v_of_rho.o : ../../Modules/xc_lda_lsda_drivers.o
v_of_rho.o : ../../Modules/xc_mgga_drivers.o
v_of_rho.o : ../../UtilXlib/mp.o
v_of_rho.o : Coul_cut_2D.o
v_of_rho.o : esm.o

99
PW/src/stres_gradcorr.f90
View File

@ -12,9 +12,10 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
!----------------------------------------------------------------------------
!
USE kinds, ONLY: DP
USE funct, ONLY: dft_is_gradient, dft_is_meta, get_igcc, &
tau_xc, tau_xc_spin, get_meta
USE xc_gga, ONLY: xc_gcx !gcxc, gcx_spin, gcc_spin, gcc_spin_more
USE funct, ONLY: dft_is_gradient, dft_is_meta, get_igcc, &
get_meta
USE xc_gga, ONLY: xc_gcx
USE xc_mgga, ONLY: xc_metagcx
USE mp_bands, ONLY: intra_bgrp_comm
USE mp, ONLY: mp_sum
USE fft_types, ONLY: fft_type_descriptor
@ -32,30 +33,29 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
!
! ... local variables
!
INTEGER :: k, l, m, ipol, is, nspin0
INTEGER :: k, l, m, ipol, is, nspin0, np
INTEGER :: nr1, nr2, nr3, nrxx, ngm
REAL(DP), ALLOCATABLE :: grho(:,:,:), rhoaux(:,:)
REAL(DP), ALLOCATABLE :: grho2(:,:), grho_ud(:),grhor2(:)
COMPLEX(DP), ALLOCATABLE :: rhogaux(:,:)
!^^^
!
REAL(DP), ALLOCATABLE :: zeta(:), rh(:)
REAL(DP) :: sx(dfft%nnr), sc(dfft%nnr)
REAL(DP), ALLOCATABLE :: v1x(:,:), v2x(:,:), rhos(:)
REAL(DP), ALLOCATABLE :: v1c(:,:), v2c(:,:), v2c_ud(:)
INTEGER, ALLOCATABLE :: null_v(:)
REAL(DP) :: vnull
REAL(DP), ALLOCATABLE :: v1x(:,:), v2x(:,:), v3x(:,:), rhos(:)
REAL(DP), ALLOCATABLE :: v1c(:,:), v2c(:,:,:), v3c(:,:), v2c_ud(:)
!
REAL(DP), PARAMETER :: epsr = 1.0d-6, epsg = 1.0d-10, e2 = 2.d0
REAL(DP) :: sigma_gradcorr(3, 3)
LOGICAL :: igcc_is_lyp
!
! ... dummy variables for meta-gga
REAL(DP) :: v3x, v3c, v3xup, v3xdw, v3cup, v3cdw
!
!
IF ( .NOT. dft_is_gradient() .AND. .NOT. dft_is_meta() ) RETURN
!
nspin0 = nspin
!
np = 1
IF (nspin0==2 .AND. dft_is_meta()) np=3
!
!if (nspin==4) nspin0 = 1
!if (nspin==4.and.domag) nspin0 = 2
IF ( nspin==4 ) CALL errore('stres_gradcorr', &
@ -77,8 +77,9 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
ALLOCATE( grho2(nrxx,nspin0) )
ALLOCATE( rhogaux(ngm,nspin0) )
!
ALLOCATE( v1x(nrxx,nspin0), v2x(nrxx,nspin0) )
ALLOCATE( v1c(nrxx,nspin0), v2c(nrxx,nspin0) )
ALLOCATE( v1x(nrxx,nspin0), v2x(nrxx,nspin0), v3x(nrxx,nspin0) )
ALLOCATE( v1c(nrxx,nspin0), v2c(np,nrxx,nspin0), v3c(nrxx,nspin0) )
!
! calculate the gradient of rho+rhocore in real space
! in LSDA case rho is temporarily converted in (up,down) format
@ -112,41 +113,49 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
! sigma_gradcor_{alpha,beta} ==
! omega^-1 \int (grad_alpha rho) ( D(rho*Exc)/D(grad_alpha rho) ) d3
!
ALLOCATE( null_v(nrxx) )
!
! routine computing v1x_v and v2x_v is different for GGA and meta-GGA
! FIXME : inefficient implementation
!
null_v = 1
grho2(:,1) = grho(1,:,1)**2 + grho(2,:,1)**2 + grho(3,:,1)**2
!
WHERE( .NOT. (ABS(rhoaux(:,1))>epsr .AND. grho2(:,1)>epsg) ) null_v(:) = 0
!
IF ( .NOT. (dft_is_meta() .AND. get_meta() /= 4) ) &
CALL xc_gcx( nrxx, nspin, rhoaux, grho, sx, sc, v1x, v2x, v1c, v2c )
CALL xc_gcx( nrxx, nspin, rhoaux, grho, sx, sc, v1x, v2x, v1c, v2c(1,:,:) )
!
DO k = 1, nrxx
IF ( dft_is_meta() .AND. get_meta() /= 4 .AND. null_v(k) /= 0 ) THEN
IF ( ABS(rhoaux(k,1))>epsr .AND. grho2(k,1)>epsg ) THEN
!
kedtau(k,1) = kedtau(k,1) / e2
CALL tau_xc( rhoaux(k,1), grho2(k,1), kedtau(k,1), sx(k), sc(k), &
v1x(k,1), v2x(k,1), v3x, v1c(k,1), v2c(k,1), v3c )
kedtau(k,1) = kedtau(k,1) * e2
!
ENDIF
ENDIF
!
DO l = 1, 3
DO m = 1, l
sigma_gradcorr(l,m) = sigma_gradcorr(l,m) + grho(l,k,1)*grho(m,k,1)* &
e2 * (v2x(k,1) + v2c(k,1)) !* null_v(k)
ENDDO
ENDDO
!
IF ( dft_is_meta() .AND. get_meta() /= 4 ) THEN
kedtau(:,1) = kedtau(:,1) / e2
CALL xc_metagcx( nrxx, 1, np, rhoaux, grho, kedtau, sx, sc, &
v1x, v2x, v3x, v1c, v2c, v3c )
kedtau(:,1) = kedtau(:,1) * e2
ENDIF
!
DO l = 1, 3
DO m = 1, l
sigma_gradcorr(l,m) = sigma_gradcorr(l,m) + SUM( grho(l,:,1)*grho(m,:,1)* &
e2 * (v2x(:,1) + v2c(1,:,1)) )
ENDDO
ENDDO
!
DEALLOCATE( null_v )
!
! DO k = 1, nrxx
! IF ( dft_is_meta() .AND. get_meta() /= 4 .AND. null_v(k) /= 0 ) THEN
! IF ( ABS(rhoaux(k,1))>epsr .AND. grho2(k,1)>epsg ) THEN
! !
! kedtau(k,1) = kedtau(k,1) / e2
! CALL tau_xc( rhoaux(k,1), grho2(k,1), kedtau(k,1), sx(k), sc(k), &
! v1x(k,1), v2x(k,1), v3x, v1c(k,1), v2c(k,1), v3c )
! kedtau(k,1) = kedtau(k,1) * e2
! !
! ENDIF
! ENDIF
! !
! DO l = 1, 3
! DO m = 1, l
! sigma_gradcorr(l,m) = sigma_gradcorr(l,m) + grho(l,k,1)*grho(m,k,1)* &
! e2 * (v2x(k,1) + v2c(k,1))
! ENDDO
! ENDDO
! !
! ENDDO
!
!
ELSEIF (nspin == 2) THEN
@ -167,7 +176,7 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
!
ALLOCATE( v2c_ud(dfft%nnr) )
!
CALL xc_gcx( nrxx, nspin, rhoaux, grho, sx, sc, v1x, v2x, v1c, v2c, v2c_ud )
CALL xc_gcx( nrxx, nspin, rhoaux, grho, sx, sc, v1x, v2x, v1c, v2c(1,:,:), v2c_ud )
!
DO l = 1, 3
DO m = 1, l
@ -179,8 +188,8 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
!
! ... correlation
sigma_gradcorr(l,m) = sigma_gradcorr(l,m) + &
SUM( grho(l,:,1) * grho(m,:,1) * v2c(:,1) + &
grho(l,:,2) * grho(m,:,2) * v2c(:,2) + &
SUM( grho(l,:,1) * grho(m,:,1) * v2c(1,:,1) + &
grho(l,:,2) * grho(m,:,2) * v2c(1,:,2) + &
(grho(l,:,1) * grho(m,:,2) + &
grho(l,:,2) * grho(m,:,1)) * v2c_ud(:) ) * e2
ENDDO
@ -196,8 +205,8 @@ SUBROUTINE stres_gradcorr( rho, rhog, rho_core, rhog_core, kedtau, nspin, &
DEALLOCATE( grho )
DEALLOCATE( grho2 )
!
DEALLOCATE( v1x, v2x )
DEALLOCATE( v1c, v2c )
DEALLOCATE( v1x, v2x, v3x )
DEALLOCATE( v1c, v2c, v3c )
!
!
DO l = 1, 3

222
PW/src/v_of_rho.f90
View File

@ -104,6 +104,8 @@ SUBROUTINE v_of_rho( rho, rho_core, rhog_core, &
RETURN
!
END SUBROUTINE v_of_rho
!
!
!----------------------------------------------------------------------------
SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
!----------------------------------------------------------------------------
@ -116,14 +118,15 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
USE gvect, ONLY : g, ngm
USE lsda_mod, ONLY : nspin
USE cell_base, ONLY : omega
USE funct, ONLY : tau_xc, tau_xc_spin, get_meta, dft_is_nonlocc, nlc
USE scf, ONLY : scf_type
USE funct, ONLY : get_meta, dft_is_nonlocc, nlc
USE xc_mgga, ONLY : xc_metagcx
USE scf, ONLY : scf_type, rhoz_or_updw
USE mp, ONLY : mp_sum
USE mp_bands, ONLY : intra_bgrp_comm
!
IMPLICIT NONE
!
TYPE (scf_type), INTENT(IN) :: rho
TYPE (scf_type), INTENT(INOUT) :: rho
!! the valence charge
REAL(DP), INTENT(IN) :: rho_core(dfftp%nnr)
!! the core charge in real space
@ -141,151 +144,121 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
! ... local variables
!
REAL(DP) :: zeta, rh, sgn(2)
INTEGER :: k, ipol, is
REAL(DP) :: ex, ec, v1x, v2x, v3x,v1c, v2c, v3c, &
& v1xup, v1xdw, v2xup, v2xdw, v1cup, v1cdw, &
& v3xup, v3xdw,v3cup, v3cdw, &
& arho, atau, fac, rhoup, rhodw, ggrho2, tauup,taudw
INTEGER :: k, ipol, is, np
!
REAL(DP), ALLOCATABLE :: ex(:), ec(:)
REAL(DP), ALLOCATABLE :: v1x(:,:), v2x(:,:), v3x(:,:)
REAL(DP), ALLOCATABLE :: v1c(:,:), v2c(:,:,:), v3c(:,:)
!
REAL(DP) :: fac
REAL(DP), DIMENSION(2) :: grho2, rhoneg
REAL(DP), DIMENSION(3) :: grhoup, grhodw, v2cup, v2cdw
REAL(DP), DIMENSION(2) :: grho2, rhoneg
REAL(DP), DIMENSION(3) :: grhoup, grhodw
!
REAL(DP), ALLOCATABLE :: grho(:,:,:), h(:,:,:), dh(:)
REAL(DP), ALLOCATABLE :: rhoout(:)
REAL(DP), ALLOCATABLE :: grho(:,:,:), h(:,:,:), dh(:)
REAL(DP), ALLOCATABLE :: rhoout(:)
COMPLEX(DP), ALLOCATABLE :: rhogsum(:)
REAL(DP), PARAMETER :: eps12 = 1.0d-12, zero=0._dp
!
!----------------------------------------------------------------------------
!
REAL(DP), PARAMETER :: eps12 = 1.0d-12, zero=0._dp
!
CALL start_clock( 'v_xc_meta' )
!
!
etxc = zero
vtxc = zero
v(:,:) = zero
etxc = zero
vtxc = zero
v(:,:) = zero
rhoneg(:) = zero
sgn(1) = 1._dp ; sgn(2) = -1._dp
fac = 1.D0 / DBLE( nspin )
fac = 1.D0 / DBLE( nspin )
np = 1
IF (nspin==2) np=3
!
ALLOCATE (grho(3,dfftp%nnr,nspin))
ALLOCATE (h(3,dfftp%nnr,nspin))
ALLOCATE (rhogsum(ngm))
ALLOCATE( grho(3,dfftp%nnr,nspin) )
ALLOCATE( h(3,dfftp%nnr,nspin) )
ALLOCATE( rhogsum(ngm) )
!
ALLOCATE( ex(dfftp%nnr), ec(dfftp%nnr) )
ALLOCATE( v1x(dfftp%nnr,nspin), v2x(dfftp%nnr,nspin) , v3x(dfftp%nnr,nspin) )
ALLOCATE( v1c(dfftp%nnr,nspin), v2c(np,dfftp%nnr,nspin), v3c(dfftp%nnr,nspin) )
!
! ... calculate the gradient of rho + rho_core in real space
! ... in LSDA case rhogsum is in (up,down) format
!
DO is = 1, nspin
!
rhogsum(:)=fac*rhog_core(:) + ( rho%of_g(:,1) + sgn(is)*rho%of_g(:,nspin) )*0.5D0
rhogsum(:) = fac*rhog_core(:) + ( rho%of_g(:,1) + sgn(is)*rho%of_g(:,nspin) )*0.5D0
!
CALL fft_gradient_g2r( dfftp, rhogsum, g, grho(1,1,is) )
!
END DO
ENDDO
DEALLOCATE(rhogsum)
!
DO k = 1, dfftp%nnr
!
DO is = 1, nspin
grho2 (is) = grho(1,k, is)**2 + grho(2,k,is)**2 + grho(3,k, is)**2
ENDDO
!
IF (nspin == 1) THEN
!
! This is the spin-unpolarised case
!
arho = ABS(rho%of_r(k, 1) )
!
atau = rho%kin_r(k,1) / e2 ! kinetic energy density in Hartree
!
IF ( (arho>eps8).AND.(grho2(1)>eps12).AND.(ABS(atau)>eps8) ) THEN
!
CALL tau_xc( arho, grho2(1),atau, ex, ec, v1x, v2x, &
v3x,v1c, v2c,v3c )
!
v(k, 1) = (v1x + v1c) * e2
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
h(:,k,1) = (v2x + v2c) * grho(:,k,1) * e2
!
kedtaur(k,1) = (v3x + v3c) * 0.5d0 * e2
!
etxc = etxc + (ex + ec) * e2 !* segno
vtxc = vtxc + (v1x+v1c) * e2 * arho
!
ELSE
h (:, k, 1) = zero
kedtaur(k,1)= zero
ENDIF
!
IF ( rho%of_r(k, 1) < zero ) rhoneg(1) = rhoneg(1) - rho%of_r (k, 1)
!
ELSE
!
! spin-polarised case
!
rhoup = ( rho%of_r(k, 1) + rho%of_r(k, 2) )*0.5d0
rhodw = ( rho%of_r(k, 1) - rho%of_r(k, 2) )*0.5d0
!
rh = rhoup + rhodw
!
do ipol=1,3
grhoup(ipol)=grho(ipol,k,1)
grhodw(ipol)=grho(ipol,k,2)
end do
!
ggrho2 = ( grho2 (1) + grho2 (2) ) * 4._dp
!
tauup = rho%kin_r(k,1) / e2
taudw = rho%kin_r(k,2) / e2
atau = tauup + taudw
!
IF ( (rh > eps8).AND.(ggrho2 > eps12).AND.(ABS(atau) > eps8) ) THEN
CALL tau_xc_spin( rhoup, rhodw, grhoup, grhodw, tauup, taudw, ex, ec, &
v1xup, v1xdw, v2xup, v2xdw, v3xup, v3xdw, v1cup, &
v1cdw, v2cup, v2cdw, v3cup, v3cdw )
!
IF (nspin == 1) THEN
!
CALL xc_metagcx( dfftp%nnr, 1, np, rho%of_r, grho, rho%kin_r/e2, ex, ec, &
v1x, v2x, v3x, v1c, v2c, v3c )
!
DO k = 1, dfftp%nnr
!
v(k,1) = (v1x(k,1)+v1c(k,1)) * e2
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
h(:,k,1) = (v2x(k,1)+v2c(1,k,1)) * grho(:,k,1) * e2
!
kedtaur(k,1) = (v3x(k,1)+v3c(k,1)) * 0.5d0 * e2
!
etxc = etxc + (ex(k)+ec(k)) * e2
vtxc = vtxc + (v1x(k,1)+v1c(k,1)) * e2 * ABS(rho%of_r(k,1))
!
IF (rho%of_r(k,1) < zero) rhoneg(1) = rhoneg(1)-rho%of_r(k,1)
!
ENDDO
!
ELSE
!
CALL rhoz_or_updw( rho, 'only_r', '->updw' )
!
CALL xc_metagcx( dfftp%nnr, 2, np, rho%of_r, grho, rho%kin_r/e2, ex, ec, &
v1x, v2x, v3x, v1c, v2c, v3c )
!
! first term of the gradient correction : D(rho*Exc)/D(rho)
!
DO k = 1, dfftp%nnr
!
v(k,1) = (v1x(k,1) + v1c(k,1)) * e2
v(k,2) = (v1x(k,2) + v1c(k,2)) * e2
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
!
IF ( get_meta()==1 .OR. get_meta()==5 ) THEN ! tpss, scan
!
! first term of the gradient correction : D(rho*Exc)/D(rho)
h(:,k,1) = (v2x(k,1) * grho(:,k,1) + v2c(:,k,1)) * e2
h(:,k,2) = (v2x(k,2) * grho(:,k,2) + v2c(:,k,2)) * e2
!
v(k, 1) = (v1xup + v1cup) * e2
v(k, 2) = (v1xdw + v1cdw) * e2
ELSE
!
! h contains D(rho*Exc)/D(|grad rho|) * (grad rho) / |grad rho|
h(:,k,1) = (v2x(k,1) + v2c(1,k,1)) * grho(:,k,1) * e2
h(:,k,2) = (v2x(k,2) + v2c(1,k,2)) * grho(:,k,2) * e2
!
IF (get_meta()==1 .OR. get_meta()==5 ) THEN ! tpss, scan
!
h(:,k,1) = (v2xup * grhoup(:) + v2cup(:)) * e2
h(:,k,2) = (v2xdw * grhodw(:) + v2cdw(:)) * e2
!
ELSE
!
h(:,k,1) = (v2xup + v2cup(1)) * grhoup(:) * e2
h(:,k,2) = (v2xdw + v2cdw(1)) * grhodw(:) * e2
!
ENDIF
!
kedtaur(k,1)= (v3xup + v3cup) * 0.5d0 * e2
kedtaur(k,2)= (v3xdw + v3cdw) * 0.5d0 * e2
!
etxc = etxc + (ex + ec) * e2
vtxc = vtxc + (v1xup+v1cup+v1xdw+v1cdw) * e2 * rh
!
ELSE
!
h(:,k,1) = zero
h(:,k,2) = zero
!
kedtaur(k,1) = zero
kedtaur(k,2) = zero
!
ENDIF
!
IF ( rhoup < zero ) rhoneg(1) = rhoneg(1) - rhoup
IF ( rhodw < zero ) rhoneg(2) = rhoneg(2) - rhodw
!
ENDIF
ENDDO
ENDIF
!
kedtaur(k,1) = (v3x(k,1) + v3c(k,1)) * 0.5d0 * e2
kedtaur(k,2) = (v3x(k,2) + v3c(k,2)) * 0.5d0 * e2
!
etxc = etxc + (ex(k)+ec(k)) * e2
vtxc = vtxc + (v1x(k,1)+v1c(k,1)+v1x(k,2)+v1c(k,2)) * e2 * (rho%of_r(k,1)+rho%of_r(k,2))
!
IF ( rho%of_r(k,1) < 0.d0 ) rhoneg(1) = rhoneg(1) - rho%of_r(k,1)
IF ( rho%of_r(k,2) < 0.d0 ) rhoneg(2) = rhoneg(2) - rho%of_r(k,2)
!
ENDDO
!
CALL rhoz_or_updw( rho, 'only_r', '->rhoz' )
!
ENDIF
!
DEALLOCATE( ex, ec )
DEALLOCATE( v1x, v2x, v3x )
DEALLOCATE( v1c, v2c, v3c )
!
!
ALLOCATE( dh( dfftp%nnr ) )
@ -334,6 +307,7 @@ SUBROUTINE v_xc_meta( rho, rho_core, rhog_core, etxc, vtxc, v, kedtaur )
!
END SUBROUTINE v_xc_meta
!
!
SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
!----------------------------------------------------------------------------
!! Exchange-Correlation potential Vxc(r) from n(r)