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

master updated

This commit is contained in:
fabrizio22 2019-05-08 12:08:43 +02:00
parent 28609e4f8b
commit 5e5884bf6a
19 changed files with 1798 additions and 1130 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
CPV/src/exch_corr.f90
View File

@ -372,23 +372,22 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
use kinds, only: DP
use funct, only: dft_is_gradient, get_igcc, &
gcxc, gcx_spin, gcc_spin, gcc_spin_more, init_lda_xc
use xc_lda_lsda, only : xc, xc_spin
use xc_lda_lsda, only : xc
implicit none
integer, intent(in) :: nnr
integer, intent(in) :: nspin
real(DP), intent(in) :: grhor( 3, nnr, nspin )
real(DP) :: h( nnr, nspin, nspin )
real(DP), intent(in) :: rhor( nnr, nspin )
real(DP) :: v( nnr, nspin )
real(DP), intent(in) :: grhor(3,nnr,nspin)
real(DP) :: h(nnr,nspin,nspin)
real(DP), intent(in) :: rhor(nnr,nspin)
real(DP) :: v(nnr,nspin)
real(DP) :: etxc
integer :: ir, is, k
real(DP) :: rup, rdw, ex, ec
real(DP), dimension(nnr) :: rhox_v,arhox_v,zeta_v
real(DP), dimension(nnr) :: ex_v, ec_v
real(DP), dimension(nnr,nspin) :: vx_v, vc_v
real(DP) :: rh, grh2, zeta
real(DP) :: rup, rdw
real(DP), dimension(nnr) :: ex, ec
real(DP), dimension(nnr,nspin) :: rhox, vx, vc
real(DP) :: rh, grh2, zeta, zetas
real(DP) :: sx, sc, v1x, v2x, v1c, v2c
real(DP) :: rhox, arhox, e2
real(DP) :: e2
real(DP) :: grho2(2), arho, segno
real(DP) :: v1xup, v1xdw, v2xup, v2xdw
real(DP) :: v1cup, v1cdw
@ -404,20 +403,18 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
e2 = 1.0d0
etxc = 0.0d0
!
call init_lda_xc()
CALL init_lda_xc()
!
if( nspin == 1 ) then
IF ( nspin == 1 ) THEN
!
! spin-unpolarized case
!
arhox_v(:) = abs(rhor(:,nspin))
CALL xc( nnr, 1, 1, rhor, ex, ec, vx, vc )
!
CALL xc( nnr, arhox_v, ex_v, ec_v, vx_v, vc_v )
v(:,1) = e2 * (vx(:,1) + vc(:,1) )
etxc = e2 * SUM( (ex + ec)*rhor(:,1) )
!
v(:,nspin) = e2 * (vx_v(:,1) + vc_v(:,1) )
etxc = e2 * SUM( (ex_v + ec_v)*rhor(:,nspin) )
!
else
ELSE
!
! spin-polarized case
!
@ -425,30 +422,26 @@ subroutine exch_corr_wrapper(nnr, nspin, grhor, rhor, etxc, v, h)
neg(2) = 0
neg(3) = 0
!
do ir = 1, nnr
rhox_v(ir) = rhor(ir,1) + rhor(ir,2)
arhox_v(ir) = abs(rhox_v(ir))
if ( arhox_v(ir) > 1.D-30 ) then
zeta_v(ir) = ( rhor(ir,1) - rhor(ir,2) ) / arhox_v(ir)
if (abs(zeta_v(ir)) > 1.d0) then
neg(3) = neg(3) + 1
zeta_v(ir) = sign(1.d0,zeta_v(ir))
endif
if (rhor(ir,1) < 0.d0) neg(1) = neg(1) + 1
if (rhor(ir,2) < 0.d0) neg(2) = neg(2) + 1
endif
enddo
rhox(:,1) = rhor(:,1) + rhor(:,2)
rhox(:,2) = rhor(:,1) - rhor(:,2)
!
call xc_spin( nnr, arhox_v, zeta_v, ex_v, ec_v, vx_v, vc_v )
CALL xc( nnr, 2, 2, rhox, ex, ec, vx, vc )
!
do ir = 1, nnr
do is = 1, nspin
v(ir,is) = e2 * (vx_v(ir,is) + vc_v(ir,is) )
enddo
etxc = etxc + e2 * (ex_v(ir) + ec_v(ir)) * rhox_v(ir)
enddo
DO ir = 1, nnr
!
DO is = 1, nspin
v(ir,is) = e2 * (vx(ir,is) + vc(ir,is))
ENDDO
etxc = etxc + e2 * (ex(ir) + ec(ir)) * rhox(ir,1)
!
zetas = rhox(ir,2) / rhox(ir,1)
IF (rhor(ir,1) < 0.d0) neg(1) = neg(1) + 1
IF (rhor(ir,2) < 0.d0) neg(2) = neg(2) + 1
IF (ABS(zetas) > 1.d0) neg(3) = neg(3) + 1
!
ENDDO
!
endif
ENDIF
if( debug_xc ) then
open(unit=17,form='unformatted')

79
LR_Modules/setup_dmuxc.f90
View File

@ -8,7 +8,7 @@
!-----------------------------------------------------------------------
SUBROUTINE setup_dmuxc
!-----------------------------------------------------------------------
!! This subroutine computes dmuxc (derivative of the XC potential)
!! This subroutine computes dmuxc (derivative of the XC potential).
!
USE kinds, ONLY : DP
USE eqv, ONLY : dmuxc
@ -17,6 +17,7 @@ SUBROUTINE setup_dmuxc
USE scf, ONLY : rho, rho_core
USE noncollin_module, ONLY : noncolin, nspin_mag
USE spin_orb, ONLY : domag
USE funct, ONLY : init_lda_xc
!
IMPLICIT NONE
!
@ -30,48 +31,54 @@ SUBROUTINE setup_dmuxc
!
ns = 1
IF ( lsda ) ns = 2
ALLOCATE(rho_aux(dfftp%nnr,ns))
IF ( (.NOT. lsda) .AND. noncolin .AND. domag ) ns = 4
!
ALLOCATE( rho_aux(dfftp%nnr,ns) )
!
dmuxc(:,:,:) = 0.d0
!
CALL init_lda_xc()
!
IF ( lsda ) THEN
!
rho_aux(:,1) = ( rho%of_r(:,1) + rho%of_r(:,2) + rho_core(:) ) * 0.5d0
rho_aux(:,2) = ( rho%of_r(:,1) - rho%of_r(:,2) + rho_core(:) ) * 0.5d0
CALL dmxc_spin(dfftp%nnr, rho_aux, dmuxc )
!
!
rho_aux(:,1) = ( rho%of_r(:,1) + rho%of_r(:,2) + rho_core(:) )*0.5_DP
rho_aux(:,2) = ( rho%of_r(:,1) - rho%of_r(:,2) + rho_core(:) )*0.5_DP
!
CALL dmxc( dfftp%nnr, 2, rho_aux, dmuxc )
!
ELSE
!
IF ( noncolin .AND. domag ) THEN
!
rho_aux(:,1) = rho%of_r (:,1) + rho_core(:)
CALL dmxc_nc( dfftp%nnr, rho_aux(:,1), rho%of_r(:,2:4), dmuxc )
!
ELSE
!
ALLOCATE(sign_r(dfftp%nnr))
rho_aux(:,1) = rho%of_r(:,1) + rho_core(:)
sign_r = 1.0_DP
DO ir = 1, dfftp%nnr
IF ( rho_aux(ir,1) < -1.d-30 ) THEN
sign_r(ir) = -1.0_DP
rho_aux(ir,1) = -rho_aux(ir,1)
ELSEIF ( rho_aux(ir,1) < 1.d-30 .AND.rho_aux(ir,1) > -1.d-30 ) THEN
sign_r(ir) = 0.0_DP
rho_aux(ir,1) = 0.5_DP
ENDIF
ENDDO
!
CALL dmxc( dfftp%nnr, rho_aux(:,1), dmuxc(:,1,1) )
dmuxc(:,1,1) = dmuxc(:,1,1)*sign_r(:)
!
DEALLOCATE(sign_r)
!
ENDIF
!
!
IF ( noncolin .AND. domag ) THEN
!
rho_aux(:,1) = rho%of_r(:,1) + rho_core(:)
rho_aux(:,2:4) = rho%of_r(:,2:4)
CALL dmxc( dfftp%nnr, 4, rho_aux, dmuxc )
!
ELSE
!
ALLOCATE(sign_r(dfftp%nnr))
rho_aux(:,1) = rho%of_r(:,1) + rho_core(:)
sign_r = 1.0_DP
DO ir = 1, dfftp%nnr
IF ( rho_aux(ir,1) < -1.d-30 ) THEN
sign_r(ir) = -1.0_DP
rho_aux(ir,1) = -rho_aux(ir,1)
ELSEIF ( rho_aux(ir,1) < 1.d-30 .AND.rho_aux(ir,1) > -1.d-30 ) THEN
sign_r(ir) = 0.0_DP
rho_aux(ir,1) = 0.5_DP
ENDIF
ENDDO
!
CALL dmxc( dfftp%nnr, 1, rho_aux, dmuxc )
dmuxc(:,1,1) = dmuxc(:,1,1)*sign_r(:)
!
DEALLOCATE(sign_r)
!
ENDIF
!
ENDIF
!
DEALLOCATE(rho_aux)
DEALLOCATE( rho_aux )
!
CALL stop_clock ('setup_dmuxc')
!

756
Modules/dmxc_drivers.f90
View File

@ -5,13 +5,170 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!---------------------------------------------------------------------
SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
!---------------------------------------------------------------------
!! Wrapper routine. Calls dmxc-driver routines from internal libraries
!! or from the external one 'libxc', depending on the input choice.
!
! Only two possibilities in the present version (LDA only):
! 1) iexch libxc + icorr libxc
! 2) iexch qe + icorr qe
!
USE kinds, ONLY: DP
#if defined(__LIBXC)
USE xc_f90_types_m
USE xc_f90_lib_m
USE xc_lda_lsda
#endif
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!! length of the I/O arrays
INTEGER, INTENT(IN) :: sr_d
!! number of spin components
REAL(DP), INTENT(IN) :: rho_in(length,sr_d)
!! charge density
REAL(DP), INTENT(OUT) :: dmuxc(length,sr_d,sr_d)
!! the derivative of the xc potential
!
! ... local variables
!
#if defined(__LIBXC)
TYPE(xc_f90_pointer_t) :: xc_func
TYPE(xc_f90_pointer_t) :: xc_info1, xc_info2
INTEGER :: pol_unpol
REAL(DP), ALLOCATABLE :: rho_lxc(:)
REAL(DP), ALLOCATABLE :: dmxc_lxc(:), dmex_lxc(:), dmcr_lxc(:)
LOGICAL :: exch_lxc_avail, corr_lxc_avail
#endif
!
INTEGER :: ir, length_lxc, length_dlxc
REAL(DP), PARAMETER :: small = 1.E-10_DP, rho_trash = 0.5_DP
!
!
#if defined(__LIBXC)
!
IF (libxc_switches(1)==1 .AND. libxc_switches(2)==1) THEN
!
length_lxc = length*sr_d
!
! ... set libxc input
SELECT CASE( sr_d )
CASE( 1 )
!
ALLOCATE( rho_lxc(length_lxc) )
pol_unpol = 1
rho_lxc = rho_in(:,1)
!
CASE( 2 )
!
ALLOCATE( rho_lxc(length_lxc) )
pol_unpol = 2
DO ir = 1, length
rho_lxc(2*ir-1) = rho_in(ir,1)
rho_lxc(2*ir) = rho_in(ir,2)
ENDDO
!
CASE( 4 )
!
! ... libxc version not available for non-collinear case with domag...
!
!
CALL dmxc_nc( length, rho_in(:,1), rho_in(:,2:4), dmuxc ) !da sistemare riconversione indici in qe da libxc
!
RETURN
!
CASE DEFAULT
!
CALL errore( 'dmxc', 'Wrong number of spin dimensions', 1 )
!
END SELECT
!
length_dlxc = length
IF (pol_unpol == 2) length_dlxc = length*3
!
!
ALLOCATE( dmex_lxc(length_dlxc), dmcr_lxc(length_dlxc), &
dmxc_lxc(length_dlxc) )
!
! ... DERIVATIVE FOR EXCHANGE
CALL xc_f90_func_init( xc_func, xc_info1, iexch_l, pol_unpol )
CALL xc_f90_lda_fxc( xc_func, length, rho_lxc(1), dmex_lxc(1) )
CALL xc_f90_func_end( xc_func )
!
! ... DERIVATIVE FOR CORRELATION
CALL xc_f90_func_init( xc_func, xc_info2, icorr_l, pol_unpol )
CALL xc_f90_lda_fxc( xc_func, length, rho_lxc(1), dmcr_lxc(1) )
CALL xc_f90_func_end( xc_func )
!
dmxc_lxc = (dmex_lxc + dmcr_lxc)*2.0_DP
!
IF (sr_d == 1) THEN
dmuxc(:,1,1) = dmxc_lxc(:)
ELSEIF (sr_d == 2) THEN
DO ir = 1, length
dmuxc(ir,1,1) = dmxc_lxc(3*ir-2)
dmuxc(ir,1,2) = dmxc_lxc(3*ir-1)
dmuxc(ir,2,1) = dmxc_lxc(3*ir-1)
dmuxc(ir,2,2) = dmxc_lxc(3*ir)
ENDDO
ENDIF
!
DEALLOCATE( dmex_lxc, dmcr_lxc, dmxc_lxc )
DEALLOCATE( rho_lxc )
!
ELSEIF (libxc_switches(1)==0 .AND. libxc_switches(2)==0 ) THEN
!
IF ( sr_d == 1 ) CALL dmxc_lda( length, rho_in(:,1), dmuxc(:,1,1) )
IF ( sr_d == 2 ) CALL dmxc_lsda( length, rho_in, dmuxc )
!
ELSE
!
CALL errore( 'dmxc', 'Derivative of exchange and correlation terms &
& must be both qe or both libxc.', 3 )
!
ENDIF
!
#else
!
SELECT CASE( sr_d )
CASE( 1 )
!
CALL dmxc_lda( length, rho_in(:,1), dmuxc(:,1,1) )
!
CASE( 2 )
!
CALL dmxc_lsda( length, rho_in, dmuxc )
!
CASE( 4 )
!
CALL dmxc_nc( length, rho_in(:,1), rho_in(:,2:4), dmuxc )
!
CASE DEFAULT
!
CALL errore( 'xc_LDA', 'Wrong ns input', 2 )
!
END SELECT
!
#endif
!
!
RETURN
!
END SUBROUTINE
!
!
!
!-----------------------------------------------------------------------
SUBROUTINE dmxc( length, rho_in, dmuxc )
SUBROUTINE dmxc_lda( length, rho_in, dmuxc )
!---------------------------------------------------------------------
!! Computes the derivative of the xc potential with respect to the
!! local density.
!
USE xc_lda_lsda, ONLY: xc, iexch_l, icorr_l
USE xc_lda_lsda, ONLY: xc_lda, iexch_l, icorr_l
USE funct, ONLY: init_lda_xc
USE kinds, ONLY: DP
!
@ -26,70 +183,77 @@ SUBROUTINE dmxc( length, rho_in, dmuxc )
!
! ... local variables
!
REAL(DP), DIMENSION(length) :: ex, vx
REAL(DP), ALLOCATABLE, DIMENSION(:) :: rs, rho, dr
REAL(DP), ALLOCATABLE, DIMENSION(:) :: ex, vx
REAL(DP), ALLOCATABLE, DIMENSION(:) :: rhoaux, dr
REAL(DP), ALLOCATABLE, DIMENSION(:) :: ec, vc
!
REAL(DP), EXTERNAL :: dpz
INTEGER :: iflg, ir
REAL(DP) :: rs, ex_s, vx_s
REAL(DP) :: dpz
INTEGER :: iflg, ir, i1, i2, f1, f2
!
REAL(DP), PARAMETER :: small = 1.E-30_DP, e2 = 2.0_DP, &
pi34 = 0.75_DP / 3.141592653589793_DP, third = 1.0_DP /3.0_DP, &
rho_trash = 0.5_DP, rs_trash = 1.0_DP
!
CALL init_lda_xc()
REAL(DP), PARAMETER :: small = 1.E-30_DP, e2 = 2.0_DP, &
pi34 = 0.75_DP/3.141592653589793_DP, &
third=1.0_DP/3.0_DP, rho_trash=0.5_DP, &
rs_trash = 1.0_DP
!
dmuxc = 0.0_DP
!
! ... first case: analytical derivatives available
!
IF (iexch_l == 1 .AND. icorr_l == 1) THEN
!
ALLOCATE( rs(length) )
!
rs = rs_trash
WHERE ( rho_in > small ) rs = (pi34 / rho_in)**third
!
! ... exchange
CALL slater( length, rs, ex, vx )
dmuxc = vx / (3.0_DP * rho_in)
!
! ... correlation
!
!$omp parallel do private( rs, ex_s, vx_s)
DO ir = 1, length
!
IF ( rho_in(ir) > small ) THEN
rs = (pi34 / rho_in(ir))**third
ELSE
dmuxc(ir) = 0.0_DP
CYCLE
ENDIF
!
CALL slater( rs, ex_s, vx_s )
dmuxc(ir) = vx_s / (3.0_DP * rho_in(ir))
!
iflg = 2
if (rs(ir) < 1.0_DP) iflg = 1
dmuxc(ir) = dmuxc(ir) + dpz(rs(ir), iflg)
if (rs < 1.0_DP) iflg = 1
dmuxc(ir) = dmuxc(ir) + dpz( rs, iflg )
!
ENDDO
!
WHERE ( rho_in < small ) dmuxc = 0.0_DP
!
DEALLOCATE( rs )
!$omp end parallel do
!
ELSE
!
! ... second case: numerical derivatives
!
ALLOCATE( ec(length), vc(length) )
ALLOCATE( dr(length), rho(length) )
ALLOCATE( ex(2*length), vx(2*length) )
ALLOCATE( ec(2*length), vc(2*length) )
ALLOCATE( dr(length), rhoaux(2*length) )
!
rho = rho_in
WHERE ( rho_in < small ) rho = rho_trash
!
i1 = 1 ; f1 = length !two blocks: [ rho+dr ]
i2 = length+1 ; f2 = 2*length ! [ rho-dr ]
!
dr(:) = MIN( 1.E-6_DP, 1.E-4_DP * rho(:) )
!WHERE ( rho_in < small ) rho = rho_trash
!
CALL xc( length, rho+dr, ex, ec, vx, vc )
dr = 0.0_DP
WHERE ( rho_in > small ) dr = MIN( 1.E-6_DP, 1.E-4_DP * rho_in )
!
dmuxc = vx + vc
rhoaux(i1:f1) = rho_in+dr
rhoaux(i2:f2) = rho_in-dr
!
CALL xc( length, rho-dr, ex, ec, vx, vc )
CALL xc_lda( length*2, rhoaux, ex, ec, vx, vc )
!
dmuxc = (dmuxc - vx - vc) / (2.0_DP * dr)
dmuxc(:) = (vx(i1:f1) + vc(i1:f1) - vx(i2:f2) - vc(i2:f2)) / &
(2.0_DP * dr(:))
!
!
DEALLOCATE( ex, vx )
DEALLOCATE( ec, vc )
DEALLOCATE( dr, rho )
DEALLOCATE( dr, rhoaux )
!
WHERE ( rho_in < small ) dmuxc = 0.0_DP
! however a higher threshold is already present in xc()
! however a higher threshold is already present in xc_lda()
!
ENDIF
!
@ -99,181 +263,174 @@ SUBROUTINE dmxc( length, rho_in, dmuxc )
!
RETURN
!
END SUBROUTINE dmxc
END SUBROUTINE dmxc_lda
!
!
!-----------------------------------------------------------------------
SUBROUTINE dmxc_spin( length, rho, dmuxc )
SUBROUTINE dmxc_lsda( length, rho_in, dmuxc )
!-----------------------------------------------------------------------
!! Computes the derivative of the xc potential with respect to the
!! local density in the spin-polarized case.
!
USE funct, ONLY: init_lda_xc
USE kinds, ONLY: DP
USE xc_lda_lsda, ONLY: xc_spin, iexch_l, icorr_l
USE xc_lda_lsda, ONLY: xc_lsda, iexch_l, icorr_l
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!! length of the input/output arrays
REAL(DP), INTENT(IN), DIMENSION(length,2) :: rho
REAL(DP), INTENT(IN), DIMENSION(length,2) :: rho_in
!! spin-up and spin-down charge density
REAL(DP), INTENT(OUT), DIMENSION(length,2,2) :: dmuxc
!! u-u, u-d, d-u, d-d derivatives of the XC functional
!
! ... local variables
!
REAL(DP), DIMENSION(length) :: zeta, rhotot
REAL(DP), DIMENSION(length) :: aux1, null_v
REAL(DP), DIMENSION(length) :: rhotot, zeta, zeta_eff
REAL(DP), DIMENSION(length) :: null_v
!
REAL(DP), ALLOCATABLE, DIMENSION(:) :: rs, zeta_eff, ds
REAL(DP), ALLOCATABLE, DIMENSION(:) :: ecu, ecp, aux2
REAL(DP), ALLOCATABLE, DIMENSION(:) :: vcu, vcp
REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: vx, vc
REAL(DP), ALLOCATABLE, DIMENSION(:) :: aux1, aux2, dr, dz
REAL(DP), ALLOCATABLE, DIMENSION(:) :: rhoaux, zetaux
REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: vx, vc, vxc
REAL(DP) :: ecu, ecp, ex_s
REAL(DP) :: vcu, vcp, vx_s
!
REAL(DP) :: fz, fz1, fz2, dmcu, dmcp, aa, bb, cc
REAL(DP) :: rs, zeta_s
!
REAL(DP), EXTERNAL :: dpz, dpz_polarized
REAL(DP) :: dpz, dpz_polarized
!
INTEGER :: ir, is, iflg
INTEGER :: i1, i2, i3, i4
INTEGER :: f1, f2, f3, f4
!
REAL(DP), PARAMETER :: small = 1.E-30_DP, e2 = 2.0_DP, &
pi34 = 0.75_DP / 3.141592653589793_DP, third = 1.0_DP/3.0_DP, &
p43 = 4.0_DP / 3.0_DP, p49 = 4.0_DP / 9.0_DP, m23 = -2.0_DP / 3.0_DP, &
rho_trash = 0.5_DP, zeta_trash = 0.5_DP
REAL(DP), PARAMETER :: small = 1.E-30_DP, e2 = 2.0_DP, &
pi34 = 0.75_DP/3.141592653589793_DP, &
third = 1.0_DP/3.0_DP, p43 = 4.0_DP/3.0_DP, &
p49 = 4.0_DP/9.0_DP, m23 = -2.0_DP/3.0_DP, &
rho_trash = 0.5_DP, zeta_trash = 0.5_DP
!
dmuxc = 0.0_DP
null_v = 1.0_DP
zeta = zeta_trash
rhotot(:) = rho(:,1) + rho(:,2)
!
CALL init_lda_xc()
!
WHERE (rhotot <= small)
rhotot = rho_trash
null_v = 0.0_DP
ELSEWHERE
zeta(:) = (rho(:,1) - rho(:,2)) / rhotot(:)
END WHERE
!
WHERE (ABS(zeta) > 1.0_DP)
zeta = zeta_trash
null_v = 0.0_DP
END WHERE
rhotot(:) = rho_in(:,1) + rho_in(:,2)
!
IF (iexch_l == 1 .AND. icorr_l == 1) THEN
!
! first case: analytical derivative available
!
ALLOCATE( rs(length) )
ALLOCATE( vcu(length), vcp(length), vx(length,1) )
ALLOCATE( ecu(length), ecp(length) )
!
! ... exchange
!
rs(:) = ( pi34 / (2.0_DP * rho(:,1)) )**third
CALL slater( length, rs, aux1, vx(:,1) )
!
dmuxc(:,1,1) = vx(:,1) / (3.0_DP * rho(:,1))
!
rs(:) = ( pi34 / (2.0_DP * rho(:,2)) )**third
CALL slater( length, rs, aux1, vx(:,1) )
!
dmuxc(:,2,2) = vx(:,1) / (3.0_DP * rho(:,2))
!
! ... correlation
!
rs = (pi34 / rhotot)**third
!
CALL pz( length, rs, 1, ecu, vcu )
CALL pz_polarized( length, rs, ecp, vcp )
! ... first case: analytical derivative available
!
!$omp parallel do private( rs, ex_s, vx_s, ecu, vcu, ecp, vcp, dmcu, dmcp, aa, bb, cc )
DO ir = 1, length
fz = ( (1.0_DP + zeta(ir))**p43 + (1.0_DP - zeta(ir))**p43 - 2.0_DP ) &
/ (2.0_DP**p43 - 2.0_DP)
fz1 = p43 * ( (1.0_DP + zeta(ir))**third - (1.0_DP - zeta(ir))**third) &
/ (2.0_DP**p43 - 2.0_DP)
fz2 = p49 * ( (1.0_DP + zeta(ir))**m23 + (1.0_DP - zeta(ir))**m23 ) &
/ (2.0_DP**p43 - 2.0_DP)
!
IF (rhotot(ir) < small) CYCLE
zeta_s = (rho_in(ir,1) - rho_in(ir,2)) / rhotot(ir)
IF (ABS(zeta_s) > 1.0_DP) CYCLE
!
! ... exchange
!
rs = ( pi34 / (2.0_DP * rho_in(ir,1)) )**third
CALL slater( rs, ex_s, vx_s )
!
dmuxc(ir,1,1) = vx_s / (3.0_DP * rho_in(ir,1))
!
rs = ( pi34 / (2.0_DP * rho_in(ir,2)) )**third
CALL slater( rs, ex_s, vx_s )
!
dmuxc(ir,2,2) = vx_s / (3.0_DP * rho_in(ir,2))
!
! ... correlation
!
rs = (pi34 / rhotot(ir))**third
!
CALL pz( rs, 1, ecu, vcu )
CALL pz_polarized( rs, ecp, vcp )
!
fz = ( (1.0_DP + zeta_s)**p43 + (1.0_DP - zeta_s)**p43 - 2.0_DP ) &
/ (2.0_DP**p43 - 2.0_DP)
fz1 = p43 * ( (1.0_DP + zeta_s)**third - (1.0_DP - zeta_s)**third) &
/ (2.0_DP**p43 - 2.0_DP)
fz2 = p49 * ( (1.0_DP + zeta_s)**m23 + (1.0_DP - zeta_s)**m23) &
/ (2.0_DP**p43 - 2.0_DP)
!
iflg = 2
IF (rs(ir) < 1.0_DP) iflg = 1
IF (rs < 1.0_DP) iflg = 1
!
dmcu = dpz( rs(ir), iflg )
dmcp = dpz_polarized( rs(ir), iflg )
dmcu = dpz( rs, iflg )
dmcp = dpz_polarized( rs, iflg )
!
aa = dmcu + fz * (dmcp - dmcu)
bb = 2.0_DP * fz1 * (vcp(ir) - vcu(ir) - (ecp(ir) - ecu(ir)) ) / rhotot(ir)
cc = fz2 * (ecp(ir) - ecu(ir)) / rhotot(ir)
bb = 2.0_DP * fz1 * (vcp - vcu - (ecp - ecu) ) / rhotot(ir)
cc = fz2 * (ecp - ecu) / rhotot(ir)
!
dmuxc(ir,1,1) = dmuxc(ir,1,1) + aa + (1.0_DP - zeta(ir)) * bb + &
(1.0_DP - zeta(ir))**2 * cc
dmuxc(ir,2,1) = dmuxc(ir,2,1) + aa + (-zeta(ir)) * bb + &
(zeta(ir)**2 - 1.0_DP) * cc
dmuxc(ir,1,1) = dmuxc(ir,1,1) + aa + (1.0_DP - zeta_s) * bb + &
(1.0_DP - zeta_s)**2 * cc
dmuxc(ir,2,1) = dmuxc(ir,2,1) + aa + (-zeta_s) * bb + &
(zeta_s**2 - 1.0_DP) * cc
dmuxc(ir,1,2) = dmuxc(ir,2,1)
dmuxc(ir,2,2) = dmuxc(ir,2,2) + aa - (1.0_DP + zeta(ir)) * bb + &
(1.0_DP + zeta(ir))**2 * cc
dmuxc(ir,2,2) = dmuxc(ir,2,2) + aa - (1.0_DP + zeta_s) * bb + &
(1.0_DP + zeta_s)**2 * cc
ENDDO
!
DEALLOCATE( rs )
DEALLOCATE( vcu, vcp, vx )
DEALLOCATE( ecu, ecp )
!$omp end parallel do
!
ELSE
!
ALLOCATE( vx(length,2), vc(length,2) )
ALLOCATE( aux2(length) )
ALLOCATE( zeta_eff(length), ds(length) )
!
! ... here ds=drho
ds(:) = MIN( 1.E-6_DP, 1.E-4_DP * rhotot(:) )
ALLOCATE( vx(4*length,2) , vc(4*length,2), vxc(2*length,2) )
ALLOCATE( rhoaux(4*length), zetaux(4*length) )
ALLOCATE( aux1(4*length) , aux2(4*length) )
ALLOCATE( dr(length), dz(length) )
!
CALL xc_spin( length, rhotot+ds, zeta, aux1, aux2, vx, vc )
i1 = 1 ; f1 = length ! four blocks: [ rho+dr , zeta ]
i2 = f1+1 ; f2 = 2*length ! [ rho-dr , zeta ]
i3 = f2+1 ; f3 = 3*length ! [ rho , zeta+dz ]
i4 = f3+1 ; f4 = 4*length ! [ rho , zeta-dz ]
!
dmuxc(:,1,1) = vx(:,1) + vc(:,1)
dmuxc(:,2,2) = vx(:,2) + vc(:,2)
!
CALL xc_spin( length, rhotot-ds, zeta, aux1, aux2, vx, vc )
dz(:) = 1.E-6_DP ! dz(:) = MIN( 1.d-6, 1.d-4*ABS(zeta(:)) )
!
dmuxc(:,1,1) = (dmuxc(:,1,1) - vx(:,1) - vc(:,1)) / (2.0_DP * ds(:))
dmuxc(:,2,1) = dmuxc(:,1,1)
dmuxc(:,2,2) = (dmuxc(:,2,2) - vx(:,2) - vc(:,2)) / (2.0_DP * ds(:))
dmuxc(:,1,2) = dmuxc(:,2,2)
! ... THRESHOLD STUFF
DO ir = 1, length
zeta_s=zeta_trash
IF (rhotot(ir) <= small) THEN
rhotot(ir)=rho_trash ; null_v(ir) = 0.0_DP
ENDIF
IF (rhotot(ir) > small) zeta_s = (rho_in(ir,1) - rho_in(ir,2)) / rhotot(ir)
zeta(ir) = zeta_s
! ... If zeta is too close to +-1, the derivative is computed at a slightly
! smaller zeta
zeta_eff(ir) = SIGN( MIN( ABS(zeta_s), (1.0_DP-2.0_DP*dz(ir)) ), zeta_s )
IF (ABS(zeta_s) > 1.0_DP) null_v(ir) = 0.0_DP
ENDDO
!
! ... now ds=dzeta
ds(:) = 1.E-6_DP
! ds() = min (1.d-6, 1.d-4 * abs(zeta(:)) )
dr(:) = MIN( 1.E-6_DP, 1.E-4_DP * rhotot(:) )
!
! ... If zeta is too close to +-1, the derivative is computed at a slightly
! smaller zeta
rhoaux(i1:f1) = rhotot + dr ; zetaux(i1:f1) = zeta
rhoaux(i2:f2) = rhotot - dr ; zetaux(i2:f2) = zeta
rhoaux(i3:f3) = rhotot ; zetaux(i3:f3) = zeta_eff + dz
rhoaux(i4:f4) = rhotot ; zetaux(i4:f4) = zeta_eff - dz
!
zeta_eff(:) = SIGN( MIN( ABS(zeta(:)), (1.0_DP-2.0_DP*ds(:)) ), zeta(:) )
!
CALL xc_spin( length, rhotot, zeta_eff+ds, aux1, aux2, vx, vc )
CALL xc_lsda( length*4, rhoaux, zetaux, aux1, aux2, vx, vc )
!
aux1 = 1.0_DP / rhotot / (2.0_DP*ds)
!
dmuxc(:,1,1) = ( vx(i1:f1,1) + vc(i1:f1,1) - vx(i2:f2,1) - vc(i2:f2,1) ) / (2.0_DP*dr)
dmuxc(:,2,2) = ( vx(i1:f1,2) + vc(i1:f1,2) - vx(i2:f2,2) - vc(i2:f2,2) ) / (2.0_DP*dr)
!
aux1(i1:f1) = 1.0_DP / rhotot(:) / (2.0_DP*dz(:))
aux1(i2:f2) = aux1(i1:f1)
DO is = 1, 2
vx(:,is) = (vx(:,is) + vc(:,is)) * aux1(:) ! vx as workspace here..
vxc(i1:f2,is) = ( vx(i3:f4,is) + vc(i3:f4,is) ) * aux1(i1:f2)
ENDDO
dmuxc(:,1,1) = dmuxc(:,1,1) + vx(:,1) * (1.0_DP-zeta(:))
dmuxc(:,1,2) = dmuxc(:,1,2) + vx(:,2) * (1.0_DP-zeta(:))
dmuxc(:,2,1) = dmuxc(:,2,1) - vx(:,1) * (1.0_DP+zeta(:))
dmuxc(:,2,2) = dmuxc(:,2,2) - vx(:,2) * (1.0_DP+zeta(:))
!
CALL xc_spin( length, rhotot, zeta_eff-ds, aux1, aux2, vx, vc )
dmuxc(:,2,1) = dmuxc(:,1,1) - (vxc(i1:f1,1) - vxc(i2:f2,1)) * (1.0_DP+zeta)
dmuxc(:,1,2) = dmuxc(:,2,2) + (vxc(i1:f1,2) - vxc(i2:f2,2)) * (1.0_DP-zeta)
dmuxc(:,1,1) = dmuxc(:,1,1) + (vxc(i1:f1,1) - vxc(i2:f2,1)) * (1.0_DP-zeta)
dmuxc(:,2,2) = dmuxc(:,2,2) - (vxc(i1:f1,2) - vxc(i2:f2,2)) * (1.0_DP+zeta)
!
aux1 = 1.0_DP / rhotot / (2.0_DP*ds)
DO is=1,2
vx(:,is) = (vx(:,is) + vc(:,is)) * aux1(:) ! ..and here
ENDDO
dmuxc(:,1,1) = dmuxc(:,1,1) - vx(:,1) * (1.0_DP-zeta(:))
dmuxc(:,1,2) = dmuxc(:,1,2) - vx(:,2) * (1.0_DP-zeta(:))
dmuxc(:,2,1) = dmuxc(:,2,1) + vx(:,1) * (1.0_DP+zeta(:))
dmuxc(:,2,2) = dmuxc(:,2,2) + vx(:,2) * (1.0_DP+zeta(:))
!
DEALLOCATE( vx, vc )
DEALLOCATE( aux2 )
DEALLOCATE( zeta_eff, ds )
DEALLOCATE( vx, vc, vxc )
DEALLOCATE( rhoaux, zetaux )
DEALLOCATE( aux1, aux2 )
DEALLOCATE( dr, dz )
!
ENDIF
!
@ -286,15 +443,16 @@ SUBROUTINE dmxc_spin( length, rho, dmuxc )
!
RETURN
!
END SUBROUTINE dmxc_spin
END SUBROUTINE dmxc_lsda
!
!
!-----------------------------------------------------------------------
SUBROUTINE dmxc_nc( length, rho_in, m, dmuxc )
!-----------------------------------------------------------------------
!! Computes the derivative of the xc potential with respect to the
!! local density and magnetization in the non-collinear case.
!
USE xc_lda_lsda, ONLY: xc_spin
USE xc_lda_lsda, ONLY: xc_lsda
USE funct, ONLY: init_lda_xc
USE kinds, ONLY: DP
!
@ -311,20 +469,22 @@ SUBROUTINE dmxc_nc( length, rho_in, m, dmuxc )
!
! ... local variables
!
REAL(DP), DIMENSION(length) :: rho, amag, zeta
REAL(DP), DIMENSION(length) :: zeta_eff, ds, null_v, null_m
REAL(DP), DIMENSION(length) :: vs, aux1, aux2
REAL(DP), DIMENSION(length,2) :: vx, vxm, vxp
REAL(DP), DIMENSION(length,2) :: vc, vcm, vcp
REAL(DP), DIMENSION(length) :: dvxc_rho, dbx_rho, dby_rho, dbz_rho
REAL(DP), DIMENSION(length) :: rhotot, amag, zeta, zeta_eff, dr, dz
REAL(DP), ALLOCATABLE, DIMENSION(:) :: rhoaux, zetaux
REAL(DP), DIMENSION(length) :: vs
INTEGER, DIMENSION(length) :: null_v
REAL(DP), ALLOCATABLE, DIMENSION(:) :: aux1, aux2
REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: vx, vc
REAL(DP), DIMENSION(length) :: dvxc_rho, dbx_rho, dby_rho, dbz_rho
!
REAL(DP) :: dvxc_mx, dvxc_my, dvxc_mz, &
dbx_mx, dbx_my, dbx_mz, &
dby_mx, dby_my, dby_mz, &
dbz_mx, dbz_my, dbz_mz
REAL(DP) :: rnull
REAL(DP) :: rnull, zeta_s
!
INTEGER :: i
INTEGER :: i1, i2, i3, i4, i5, i
INTEGER :: f1, f2, f3, f4, f5
!
REAL(DP), PARAMETER :: small = 1.E-30_DP, e2 = 2.0_DP, &
rho_trash = 0.5_DP, zeta_trash = 0.5_DP, &
@ -332,133 +492,129 @@ SUBROUTINE dmxc_nc( length, rho_in, m, dmuxc )
!
dmuxc = 0.0_DP
!
rho = rho_in
ALLOCATE( rhoaux(length*5), zetaux(length*5) )
ALLOCATE( aux1(length*5), aux2(length*5) )
ALLOCATE( vx(length*5,2), vc(length*5,2) )
!
rhotot = rho_in
zeta = zeta_trash
amag = amag_trash
null_v = 1.0_DP
null_m = 1.0_DP
null_v = 1
!
CALL init_lda_xc()
i1 = 1 ; f1 = length ! five blocks: [ rho , zeta ]
i2 = f1+1 ; f2 = 2*length ! [ rho+dr , zeta ]
i3 = f2+1 ; f3 = 3*length ! [ rho-dr , zeta ]
i4 = f3+1 ; f4 = 4*length ! [ rho , zeta+dz ]
i5 = f4+1 ; f5 = 5*length ! [ rho , zeta-dz ]
!
WHERE (rho_in <= small)
rho = rho_trash
null_v = 0.0_DP
ELSEWHERE
amag = SQRT( m(:,1)**2 + m(:,2)**2 + m(:,3)**2 )
zeta = amag / rho
END WHERE
!
WHERE (ABS(zeta) > 1.0_DP)
zeta = zeta_trash
null_v = 0.0_DP
END WHERE
!
CALL xc_spin( length, rho, zeta, aux1, aux2, vx, vc )
!
vs = 0.5_DP*( vx(:,1)+vc(:,1)-vx(:,2)-vc(:,2) )
!
! ... here ds=drho
ds = MIN( 1.E-6_DP, 1.E-4_DP * rho )
!
CALL xc_spin( length, rho-ds, zeta, aux1, aux2, vxm, vcm )
!
CALL xc_spin( length, rho+ds, zeta, aux1, aux2, vxp, vcp )
!
dvxc_rho = ((vxp(:,1) + vcp(:,1) - vxm(:,1) - vcm(:,1)) + &
(vxp(:,2) + vcp(:,2) - vxm(:,2) - vcm(:,2))) / (4.0_DP*ds)
!
WHERE (amag < 1.E-10_DP)
rho = 0.5_DP
zeta = 0.5_DP
amag = 0.025_DP
null_m = 0.0_DP
END WHERE
!
!
aux2(:) = vxp(:,1) + vcp(:,1) - vxm(:,1) - vcm(:,1) - &
( vxp(:,2) + vcp(:,2) - vxm(:,2) - vcm(:,2) )
!
dbx_rho(:) = aux2 * m(:,1) / (4.0_DP*ds*amag)
dby_rho(:) = aux2 * m(:,2) / (4.0_DP*ds*amag)
dbz_rho(:) = aux2 * m(:,3) / (4.0_DP*ds*amag)
!
! ... now ds=dzeta
ds = 1.0E-6_DP
! ds = min (1.d-6, 1.d-4 * abs (zeta) )
!
! ... If zeta is too close to +-1, the derivative is computed at a slightly
! smaller zeta
dz = 1.0E-6_DP !dz = MIN( 1.d-6, 1.d-4*ABS(zeta) )
!
DO i = 1, length
zeta_eff(i) = SIGN( MIN( ABS( zeta(i) ), ( 1.0_DP - 2.0_DP*ds(i) ) ) , zeta(i) )
zeta_s = zeta_trash
IF (rhotot(i) <= small) THEN
rhotot(i) = rho_trash ; null_v(i) = 0.0_DP
ENDIF
amag(i) = SQRT( m(i,1)**2 + m(i,2)**2 + m(i,3)**2 )
IF (rhotot(i) > small) zeta_s = amag(i) / rhotot(i)
zeta(i) = zeta_s
zeta_eff(i) = SIGN( MIN( ABS(zeta_s), (1.0_DP-2.0_DP*dz(i)) ), zeta_s )
IF (ABS(zeta_s) > 1.0_DP) null_v(i) = 0
ENDDO
!
CALL xc_spin( length, rho, zeta_eff-ds, aux1, aux2, vxm, vcm )
dr = MIN( 1.E-6_DP, 1.E-4_DP * rhotot )
!
rhoaux(i1:f1) = rhotot ; zetaux(i1:f1) = zeta
rhoaux(i2:f2) = rhotot + dr ; zetaux(i2:f2) = zeta
rhoaux(i3:f3) = rhotot - dr ; zetaux(i3:f3) = zeta
rhoaux(i4:f4) = rhotot ; zetaux(i4:f4) = zeta_eff + dz
rhoaux(i5:f5) = rhotot ; zetaux(i5:f5) = zeta_eff - dz
!
CALL xc_spin( length, rho, zeta_eff+ds, aux1, aux2, vxp, vcp )
!
! ... The variables are rho and m, so zeta depends on rho
CALL xc_lsda( length*5, rhoaux, zetaux, aux1, aux2, vx, vc )
!
aux1(:) = vxp(:,1) + vcp(:,1) - vxm(:,1) - vcm(:,1) + &
vxp(:,2) + vcp(:,2) - vxm(:,2) - vcm(:,2)
aux2(:) = vxp(:,1) + vcp(:,1) - vxm(:,1) - vcm(:,1) - &
( vxp(:,2) + vcp(:,2) - vxm(:,2) - vcm(:,2) )
!
dvxc_rho = dvxc_rho - aux1 * zeta/rho / (4.0_DP*ds) * null_m
dbx_rho(:) = dbx_rho - aux2 * m(:,1) * zeta/rho / (4.0_DP*ds*amag)
dby_rho(:) = dby_rho - aux2 * m(:,2) * zeta/rho / (4.0_DP*ds*amag)
dbz_rho(:) = dbz_rho - aux2 * m(:,3) * zeta/rho / (4.0_DP*ds*amag)
vs(:) = 0.5_DP*( vx(i1:f1,1)+vc(i1:f1,1)-vx(i1:f1,2)-vc(i1:f1,2) )
!
dmuxc(:,1,1) = dvxc_rho * null_v
dmuxc(:,2,1) = dbx_rho * null_v
dmuxc(:,3,1) = dby_rho * null_v
dmuxc(:,4,1) = dbz_rho * null_v
dvxc_rho(:) = ((vx(i2:f2,1) + vc(i2:f2,1) - vx(i3:f3,1) - vc(i3:f3,1)) + &
(vx(i2:f2,2) + vc(i2:f2,2) - vx(i3:f3,2) - vc(i3:f3,2))) / (4.0_DP*dr)
!
aux2(1:length) = vx(i2:f2,1) + vc(i2:f2,1) - vx(i3:f3,1) - vc(i3:f3,1) - &
( vx(i2:f2,2) + vc(i2:f2,2) - vx(i3:f3,2) - vc(i3:f3,2) )
!
! ... Here the derivatives with respect to m
dbx_rho(:) = aux2(1:length) * m(:,1) / (4.0_DP*dr*amag)
dby_rho(:) = aux2(1:length) * m(:,2) / (4.0_DP*dr*amag)
dbz_rho(:) = aux2(1:length) * m(:,3) / (4.0_DP*dr*amag)
!
aux1(1:length) = vx(i4:f4,1) + vc(i4:f4,1) - vx(i5:f5,1) - vc(i5:f5,1) + &
vx(i4:f4,2) + vc(i4:f4,2) - vx(i5:f5,2) - vc(i5:f5,2)
aux2(1:length) = vx(i4:f4,1) + vc(i4:f4,1) - vx(i5:f5,1) - vc(i5:f5,1) - &
( vx(i4:f4,2) + vc(i4:f4,2) - vx(i5:f5,2) - vc(i5:f5,2) )
!
DO i = 1, length
!
IF ( null_v(i) == 0 ) THEN
dmuxc(i,:,:) = 0.0_DP
CYCLE
ENDIF
!
IF (amag(i) <= 1.E-10_DP) THEN
dmuxc(i,1,1) = dvxc_rho(i)
CYCLE
ENDIF
!
dvxc_rho(i) = dvxc_rho(i) - aux1(i) * zeta(i)/rhotot(i) / (4.0_DP*dz(i))
dbx_rho(i) = dbx_rho(i) - aux2(i) * m(i,1) * zeta(i)/rhotot(i) / (4.0_DP*dz(i)*amag(i))
dby_rho(i) = dby_rho(i) - aux2(i) * m(i,2) * zeta(i)/rhotot(i) / (4.0_DP*dz(i)*amag(i))
dbz_rho(i) = dbz_rho(i) - aux2(i) * m(i,3) * zeta(i)/rhotot(i) / (4.0_DP*dz(i)*amag(i))
!
dmuxc(i,1,1) = dvxc_rho(i)
dmuxc(i,2,1) = dbx_rho(i)
dmuxc(i,3,1) = dby_rho(i)
dmuxc(i,4,1) = dbz_rho(i)
!
! ... Here the derivatives with respect to m
!
rnull = null_v(i)
!
dvxc_mx = aux1(i) * m(i,1) / rho(i) / (4.0_DP*ds(i)*amag(i))
dvxc_my = aux1(i) * m(i,2) / rho(i) / (4.0_DP*ds(i)*amag(i))
dvxc_mz = aux1(i) * m(i,3) / rho(i) / (4.0_DP*ds(i)*amag(i))
dvxc_mx = aux1(i) * m(i,1) / rhotot(i) / (4.0_DP*dz(i)*amag(i))
dvxc_my = aux1(i) * m(i,2) / rhotot(i) / (4.0_DP*dz(i)*amag(i))
dvxc_mz = aux1(i) * m(i,3) / rhotot(i) / (4.0_DP*dz(i)*amag(i))
!
dbx_mx = (aux2(i) * m(i,1) * m(i,1) * amag(i)/rho(i) / (4.0_DP*ds(i)) + &
dbx_mx = (aux2(i) * m(i,1) * m(i,1) * amag(i)/rhotot(i) / (4.0_DP*dz(i)) + &
vs(i) * (m(i,2)**2+m(i,3)**2)) / amag(i)**3
dbx_my = (aux2(i) * m(i,1) * m(i,2) * amag(i)/rho(i) / (4.0_DP*ds(i)) - &
dbx_my = (aux2(i) * m(i,1) * m(i,2) * amag(i)/rhotot(i) / (4.0_DP*dz(i)) - &
vs(i) * m(i,1) * m(i,2) ) / amag(i)**3
dbx_mz = (aux2(i) * m(i,1) * m(i,3) * amag(i)/rho(i) / (4.0_DP*ds(i)) - &
dbx_mz = (aux2(i) * m(i,1) * m(i,3) * amag(i)/rhotot(i) / (4.0_DP*dz(i)) - &
vs(i) * m(i,1) * m(i,3) ) / amag(i)**3
!
dby_mx = dbx_my
dby_my = (aux2(i) * m(i,2) * m(i,2) * amag(i)/rho(i) / (4.0_DP*ds(i)) + &
dby_my = (aux2(i) * m(i,2) * m(i,2) * amag(i)/rhotot(i) / (4.0_DP*dz(i)) + &
vs(i) * (m(i,1)**2 + m(i,3)**2)) / amag(i)**3
dby_mz = (aux2(i) * m(i,2) * m(i,3) * amag(i)/rho(i) / (4.0_DP*ds(i)) - &
dby_mz = (aux2(i) * m(i,2) * m(i,3) * amag(i)/rhotot(i) / (4.0_DP*dz(i)) - &
vs(i) * m(i,2) * m(i,3)) / amag(i)**3
!
dbz_mx = dbx_mz
dbz_my = dby_mz
dbz_mz = (aux2(i) * m(i,3) * m(i,3) * amag(i)/rho(i) / (4.0_DP*ds(i)) + &
dbz_mz = (aux2(i) * m(i,3) * m(i,3) * amag(i)/rhotot(i) / (4.0_DP*dz(i)) + &
vs(i)*(m(i,1)**2 + m(i,2)**2)) / amag(i)**3
!
! ... assigns values to dmuxc and sets to zero trash points
!
dmuxc(i,1,2) = dvxc_mx * rnull
dmuxc(i,1,3) = dvxc_my * rnull
dmuxc(i,1,4) = dvxc_mz * rnull
dmuxc(i,1,2) = dvxc_mx
dmuxc(i,1,3) = dvxc_my
dmuxc(i,1,4) = dvxc_mz
!
dmuxc(i,2,2) = dbx_mx * rnull
dmuxc(i,2,3) = dbx_my * rnull
dmuxc(i,2,4) = dbx_mz * rnull
dmuxc(i,2,2) = dbx_mx
dmuxc(i,2,3) = dbx_my
dmuxc(i,2,4) = dbx_mz
!
dmuxc(i,3,2) = dby_mx * rnull
dmuxc(i,3,3) = dby_my * rnull
dmuxc(i,3,4) = dby_mz * rnull
dmuxc(i,3,2) = dby_mx
dmuxc(i,3,3) = dby_my
dmuxc(i,3,4) = dby_mz
!
dmuxc(i,4,2) = dbz_mx * rnull
dmuxc(i,4,3) = dbz_my * rnull
dmuxc(i,4,4) = dbz_mz * rnull
dmuxc(i,4,2) = dbz_mx
dmuxc(i,4,3) = dbz_my
dmuxc(i,4,4) = dbz_mz
!
ENDDO
!
@ -466,6 +622,96 @@ SUBROUTINE dmxc_nc( length, rho_in, m, dmuxc )
!
dmuxc = e2 * dmuxc
!
DEALLOCATE( rhoaux, zetaux)
DEALLOCATE( aux1, aux2 )
DEALLOCATE( vx, vc )
!
RETURN
!
END SUBROUTINE dmxc_nc
!
!
!-----------------------------------------------------------------------
FUNCTION dpz( rs, iflg )
!-----------------------------------------------------------------------
!! Derivative of the correlation potential with respect to local density
!! Perdew and Zunger parameterization of the Ceperley-Alder functional.
!
USE kinds, ONLY: DP
USE constants, ONLY: pi, fpi
!
IMPLICIT NONE
!
REAL(DP), INTENT(IN) :: rs
INTEGER, INTENT(IN) :: iflg
REAL(DP) :: dpz
!
! ... local variables
! a,b,c,d,gc,b1,b2 are the parameters defining the functional
!
REAL(DP), PARAMETER :: a = 0.0311d0, b = -0.048d0, c = 0.0020d0, &
d = -0.0116d0, gc = -0.1423d0, b1 = 1.0529d0, b2 = 0.3334d0,&
a1 = 7.0d0 * b1 / 6.d0, a2 = 4.d0 * b2 / 3.d0
REAL(DP) :: x, den, dmx, dmrs
!
IF (iflg == 1) THEN
dmrs = a / rs + 2.d0 / 3.d0 * c * (LOG(rs) + 1.d0) + &
(2.d0 * d-c) / 3.d0
ELSE
x = SQRT(rs)
den = 1.d0 + x * (b1 + x * b2)
dmx = gc * ( (a1 + 2.d0 * a2 * x) * den - 2.d0 * (b1 + 2.d0 * &
b2 * x) * (1.d0 + x * (a1 + x * a2) ) ) / den**3
dmrs = 0.5d0 * dmx / x
ENDIF
!
dpz = - fpi * rs**4.d0 / 9.d0 * dmrs
!
RETURN
!
END FUNCTION dpz
!
!
!-----------------------------------------------------------------------
FUNCTION dpz_polarized( rs, iflg )
!-----------------------------------------------------------------------
!! Derivative of the correlation potential with respect to local density
!! Perdew and Zunger parameterization of the Ceperley-Alder functional. |
!! Spin-polarized case.
!
USE kinds, ONLY: DP
USE constants, ONLY: pi, fpi
!
IMPLICIT NONE
!
REAL(DP), INTENT(IN) :: rs
INTEGER, INTENT(IN) :: iflg
REAL(DP) :: dpz_polarized
!
! ... local variables
!
! a,b,c,d,gc,b1,b2 are the parameters defining the functional
!
REAL(DP), PARAMETER :: a=0.01555_DP, b=-0.0269_DP, c=0.0007_DP, &
d=-0.0048_DP, gc=-0.0843_DP, b1=1.3981_DP,&
b2=0.2611_DP, a1=7.0_DP*b1/6._DP, a2=4._DP*b2/3._DP
REAL(DP) :: x, den, dmx, dmrs
!
!
IF (iflg == 1) THEN
dmrs = a/rs + 2._DP/3._DP * c * (LOG(rs) + 1._DP) + &
(2._DP*d - c)/3._DP
ELSE
x = SQRT(rs)
den = 1._DP + x * (b1 + x*b2)
dmx = gc * ( (a1 + 2._DP * a2 * x) * den - 2._DP * (b1 + 2._DP * &
b2 * x) * (1._DP + x * (a1 + x*a2) ) ) / den**3
dmrs = 0.5d0 * dmx/x
ENDIF
!
dpz_polarized = - fpi * rs**4._DP / 9._DP * dmrs
!
!
RETURN
!
END FUNCTION dpz_polarized

105
Modules/funct.f90
View File

@ -75,6 +75,7 @@ module funct
! PRIVATE variables defining the DFT functional
!
PRIVATE :: dft, iexch, icorr, igcx, igcc, imeta, inlc
PRIVATE :: qe_or_libxc
PRIVATE :: discard_input_dft
PRIVATE :: isgradient, ismeta, ishybrid
PRIVATE :: exx_fraction, exx_started
@ -314,6 +315,9 @@ module funct
!
integer, parameter:: notset = -1
!
! switches to decide between qe (1) and libxc (2) routines
integer :: qe_or_libxc(2)
!
! internal indices for exchange-correlation
! iexch: type of exchange
! icorr: type of correlation
@ -407,6 +411,14 @@ CONTAINS
do l = 1, len
dftout (l:l) = capital (dft_(l:l) )
enddo
!
!
! PROVISIONAL: it should be put as an environment variable or something
qe_or_libxc(:)=0
#if defined(__LIBXC)
qe_or_libxc(1)=1
qe_or_libxc(2)=1
#endif
!
! ----------------------------------------------
! FIRST WE CHECK ALL THE SHORT NAMES
@ -416,6 +428,7 @@ CONTAINS
! special cases : PZ (LDA is equivalent to PZ)
IF (('PZ' .EQ. TRIM(dftout) ).OR.('LDA' .EQ. TRIM(dftout) )) THEN
dft_defined = set_dft_values(1,1,0,0,0,0)
! speciale cases : PW ( LDA with PW correlation )
ELSE IF ( 'PW' .EQ. TRIM(dftout)) THEN
dft_defined = set_dft_values(1,4,0,0,0,0)
@ -671,9 +684,9 @@ CONTAINS
igcc = matching (dftout, ngcc,gradc)
imeta = matching (dftout,nmeta, meta)
inlc = matching (dftout, ncnl, nonlocc)
endif
! ----------------------------------------------------------------
! Last check
! No more defaults, the code exits if the dft is not defined
@ -805,24 +818,86 @@ CONTAINS
end subroutine set_auxiliary_flags
!
!-----------------------------------------------------------------------
logical function set_dft_values (i1,i2,i3,i4,i5,i6)
logical function set_dft_values(i1,i2,i3,i4,i5,i6)
!-----------------------------------------------------------------------
!
implicit none
integer :: i1,i2,i3,i4,i5,i6
!
iexch=i1
icorr=i2
igcx =i3
igcc =i4
inlc =i5
imeta=i6
!
set_dft_values = .true.
!
return
!
end function set_dft_values
!
!
!--------------------------------------------------------------------------
FUNCTION qe_to_libxc_index( index_qe, term )
!-----------------------------------------------------------------------
!! PROVISIONAL: converts q-e indexes of functionals into libxc ones,
!! when possible.
!! In the next commit this will be done directly in 'funct.f90' and
!! including all the cases available.
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: index_qe
!! index of q-e functional term
CHARACTER(LEN=*), INTENT(IN) :: term
!! which term to consider
INTEGER :: qe_to_libxc_index
!! index converted in libxc notation
!
! ... exchange term
IF (term .EQ. 'exch_LDA') THEN
!
qe_or_libxc(1) = 1
!
SELECT CASE( index_qe )
CASE( 1 )
qe_to_libxc_index = 1
CASE DEFAULT
qe_to_libxc_index = -1
qe_or_libxc(1) = 0
END SELECT
!
! ... correlation term
ELSEIF (term .EQ. 'corr_LDA') THEN
!
qe_or_libxc(2) = 1
!
SELECT CASE( index_qe )
CASE( 1 ) !pz
qe_to_libxc_index = 9
CASE( 5 ) !wigner
qe_to_libxc_index = 2
CASE( 2 ) !vwn
qe_to_libxc_index = 7
CASE( 4 ) !pw
qe_to_libxc_index = 12
CASE DEFAULT
qe_to_libxc_index = -1
qe_or_libxc(2) = 0
END SELECT
!
ELSE
!
CALL errore( 'qe_to_libxc_index', 'ERROR: option not available', 3 )
!
ENDIF
!
RETURN
!
END FUNCTION qe_to_libxc_index
!
!
!-----------------------------------------------------------------------
subroutine enforce_input_dft (dft_, nomsg)
!
@ -2844,16 +2919,21 @@ end subroutine tau_xc_array_spin
SUBROUTINE init_lda_xc()
!! Gets from inside parameters needed to initialize lda xc-drivers.
!
USE kinds, ONLY: DP
USE xc_lda_lsda, ONLY: iexch_l, icorr_l, &
exx_started_l, is_there_finite_size_corr, &
exx_fraction_l, finite_size_cell_volume_l
USE kinds, ONLY: DP
USE xc_lda_lsda, ONLY: libxc_switches, iexch_l, icorr_l, &
exx_started_l, is_there_finite_size_corr, &
exx_fraction_l, finite_size_cell_volume_l
!
IMPLICIT NONE
!
! use qe or libxc for the different terms (0: qe, 1: libxc)
libxc_switches(:) = qe_or_libxc(:)
!
! exchange-correlation indexes
iexch_l = get_iexch()
IF (libxc_switches(1)==1) iexch_l = qe_to_libxc_index( iexch, 'exch_LDA' )
icorr_l = get_icorr()
IF (libxc_switches(2)==1) icorr_l = qe_to_libxc_index( icorr, 'corr_LDA' )
!
! hybrid exchange vars
exx_started_l = exx_is_active()
@ -2861,7 +2941,8 @@ SUBROUTINE init_lda_xc()
IF ( exx_started ) exx_fraction_l = get_exx_fraction()
!
! finite size correction vars
CALL get_finite_size_cell_volume( is_there_finite_size_corr, finite_size_cell_volume_l )
CALL get_finite_size_cell_volume( is_there_finite_size_corr, &
finite_size_cell_volume_l )
!
RETURN
!

90
Modules/functionals.f90
View File

@ -447,8 +447,7 @@ subroutine ggac( rho, grho, sc, v1c, v2c )
!
USE kinds, ONLY : DP
implicit none
integer, parameter :: length=1 !^^^ PROVISIONAL (xc-lib)
real(DP), dimension(length) :: rs, ec, vc
real(DP) :: rs, ec, vc
real(DP) :: rho, grho, sc, v1c, v2c
real(DP) :: al, pa, pb, pc, pd, cx, cxc0, cc0
parameter (al = 0.09d0, pa = 0.023266d0, pb = 7.389d-6, pc = &
@ -465,36 +464,36 @@ subroutine ggac( rho, grho, sc, v1c, v2c )
real(DP) :: h0, dh0, ddh0, ee, cn, dcn, cna, dcna, cnb, dcnb, h1, &
dh1, ddh1
!
rs(1) = pi34 / rho**third
rs2 = rs(1) * rs(1)
rs3 = rs(1) * rs2
rs = pi34 / rho**third
rs2 = rs * rs
rs3 = rs * rs2
!
call pw( length, rs, 1, ec, vc )
call pw( rs, 1, ec, vc )
!
kf = xkf / rs(1)
kf = xkf / rs
ks = xks * sqrt(kf)
t = sqrt(grho) / (2.d0 * ks * rho)
expe = exp( - 2.d0 * al * ec(1) / (be * be) )
expe = exp( - 2.d0 * al * ec / (be * be) )
af = 2.d0 * al / be * (1.d0 / (expe-1.d0) )
bf = expe * (vc(1) - ec(1))
bf = expe * (vc - ec)
y = af * t * t
xy = (1.d0 + y) / (1.d0 + y + y * y)
qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
s1 = 1.d0 + 2.d0 * al / be * t * t * xy
h0 = be * be / (2.d0 * al) * log(s1)
dh0 = be * t * t / s1 * ( - 7.d0 / 3.d0 * xy - qy * (af * bf / &
be-7.d0 / 3.d0) )
be-7.d0 / 3.d0) )
ddh0 = be / (2.d0 * ks * ks * rho) * (xy - qy) / s1
ee = - 100.d0 * (ks / kf * t) **2
cna = cxc0 + pa * rs(1) + pb * rs2
dcna = pa * rs(1) + 2.d0 * pb * rs2
cnb = 1.d0 + pc * rs(1) + pd * rs2 + 1.d4 * pb * rs3
dcnb = pc * rs(1) + 2.d0 * pd * rs2 + 3.d4 * pb * rs3
cna = cxc0 + pa * rs + pb * rs2
dcna = pa * rs + 2.d0 * pb * rs2
cnb = 1.d0 + pc * rs + pd * rs2 + 1.d4 * pb * rs3
dcnb = pc * rs + 2.d0 * pd * rs2 + 3.d4 * pb * rs3
cn = cna / cnb - cx
dcn = dcna / cnb - cna * dcnb / (cnb * cnb)
h1 = nu * (cn - cc0 - 3.d0 / 7.d0 * cx) * t * t * exp(ee)
dh1 = - third * (h1 * (7.d0 + 8.d0 * ee) + nu * t * t * exp(ee) &
* dcn)
* dcn)
ddh1 = 2.d0 * h1 * (1.d0 + ee) * rho / grho
sc = rho * (h0 + h1)
v1c = h0 + h1 + dh0 + dh1
@ -774,29 +773,28 @@ subroutine pbec(rho, grho, iflag, sc, v1c, v2c)
real(DP), parameter :: third = 1.d0 / 3.d0, pi34 = 0.6203504908994d0
real(DP), parameter :: xkf = 1.919158292677513d0, xks = 1.128379167095513d0
! pi34=(3/4pi)^(1/3), xkf=(9 pi/4)^(1/3), xks= sqrt(4/pi)
integer, parameter :: length=1 !^^^ PROVISIONAL
real(DP), dimension(length) :: rs, ec, vc
real(DP) :: rs, ec, vc
!
real(DP) :: kf, ks, t, expe, af, bf, y, xy, qy
real(DP) :: s1, h0, dh0, ddh0, sc2D, v1c2D, v2c2D
!
rs(1) = pi34 / rho**third
rs = pi34 / rho**third
!
call pw( length, rs, 1, ec, vc )
call pw( rs, 1, ec, vc )
!
kf = xkf / rs(1)
ks = xks * sqrt (kf)
t = sqrt (grho) / (2.d0 * ks * rho)
expe = exp ( - ec(1) / ga)
kf = xkf / rs
ks = xks * sqrt(kf)
t = sqrt(grho) / (2.d0 * ks * rho)
expe = exp( - ec / ga)
af = be(iflag) / ga * (1.d0 / (expe-1.d0) )
bf = expe * (vc(1) - ec(1))
bf = expe * (vc - ec)
y = af * t * t
xy = (1.d0 + y) / (1.d0 + y + y * y)
qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
s1 = 1.d0 + be(iflag) / ga * t * t * xy
h0 = ga * log (s1)
h0 = ga * log(s1)
dh0 = be(iflag) * t * t / s1 * ( - 7.d0 / 3.d0 * xy - qy * (af * bf / &
be(iflag)-7.d0 / 3.d0) )
be(iflag)-7.d0 / 3.d0) )
ddh0 = be(iflag) / (2.d0 * ks * ks * rho) * (xy - qy) / s1
sc = rho * h0
v1c = h0 + dh0
@ -1352,46 +1350,6 @@ subroutine wcx (rho, grho, sx, v1x, v2x)
return
end subroutine wcx
!
!-----------------------------------------------------------------------
function dpz(rs, iflg)
!-----------------------------------------------------------------------
! derivative of the correlation potential with respect to local density
! Perdew and Zunger parameterization of the Ceperley-Alder functional
!
use kinds, only: DP
USE constants, ONLY: pi, fpi
!
implicit none
!
real(DP), intent(in) :: rs
integer, intent(in) :: iflg
real(DP) :: dpz
!
! local variables
! a,b,c,d,gc,b1,b2 are the parameters defining the functional
!
real(DP), parameter :: a = 0.0311d0, b = -0.048d0, c = 0.0020d0, &
d = -0.0116d0, gc = -0.1423d0, b1 = 1.0529d0, b2 = 0.3334d0,&
a1 = 7.0d0 * b1 / 6.d0, a2 = 4.d0 * b2 / 3.d0
real(DP) :: x, den, dmx, dmrs
!
!
if (iflg == 1) then
dmrs = a / rs + 2.d0 / 3.d0 * c * (log (rs) + 1.d0) + &
(2.d0 * d-c) / 3.d0
else
x = sqrt (rs)
den = 1.d0 + x * (b1 + x * b2)
dmx = gc * ( (a1 + 2.d0 * a2 * x) * den - 2.d0 * (b1 + 2.d0 * &
b2 * x) * (1.d0 + x * (a1 + x * a2) ) ) / den**3
dmrs = 0.5d0 * dmx / x
endif
!
dpz = - fpi * rs**4.d0 / 9.d0 * dmrs
return
!
end function dpz
!----------------------------------------------------------------------
!
! HSE (wPBE) stabbing starts HERE
!

110
Modules/lsda_functionals.f90
View File

@ -101,9 +101,8 @@ subroutine ggac_spin(rho, zeta, grho, sc, v1cup, v1cdw, v2c)
implicit none
real(DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c
!
integer, parameter :: length=1 !^^^ PROVISIONAL
real(dp), dimension(length) :: rs , zeta_v, ec
real(dp), dimension(length,2) :: vc
real(dp) :: rs, ec
real(dp) :: vc(2)
!
real(DP) :: al, pa, pb, pc, pd, cx, cxc0, cc0
parameter (al = 0.09d0, pa = 0.023266d0, pb = 7.389d-6, pc = &
@ -121,16 +120,14 @@ subroutine ggac_spin(rho, zeta, grho, sc, v1cup, v1cdw, v2c)
ddh1, fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, dh0zup, &
dh0zdw, dh1zup, dh1zdw
!
rs(1) = pi34 / rho**third
rs2 = rs(1) * rs(1)
rs3 = rs(1) * rs2
rs = pi34 / rho**third
rs2 = rs * rs
rs3 = rs * rs2
!
zeta_v(1) = zeta
call pw_spin( rs, zeta, ec, vc )
!
call pw_spin( length, rs, zeta_v, ec, vc )
!
kf = xkf / rs(1)
ks = xks * sqrt (kf)
kf = xkf / rs
ks = xks * sqrt(kf)
fz = 0.5d0 * ( (1.d0 + zeta) ** (2.d0 / 3.d0) + (1.d0 - zeta) ** ( &
2.d0 / 3.d0) )
fz2 = fz * fz
@ -138,35 +135,35 @@ subroutine ggac_spin(rho, zeta, grho, sc, v1cup, v1cdw, v2c)
fz4 = fz3 * fz
dfz = ( (1.d0 + zeta) ** ( - 1.d0 / 3.d0) - (1.d0 - zeta) ** ( - &
1.d0 / 3.d0) ) / 3.d0
t = sqrt (grho) / (2.d0 * fz * ks * rho)
expe = exp ( - 2.d0 * al * ec(1) / (fz3 * be * be) )
t = sqrt(grho) / (2.d0 * fz * ks * rho)
expe = exp( - 2.d0 * al * ec / (fz3 * be * be) )
af = 2.d0 * al / be * (1.d0 / (expe-1.d0) )
bfup = expe * (vc(1,1) - ec(1)) / fz3
bfdw = expe * (vc(1,2) - ec(1)) / fz3
bfup = expe * (vc(1) - ec) / fz3
bfdw = expe * (vc(2) - ec) / fz3
y = af * t * t
xy = (1.d0 + y) / (1.d0 + y + y * y)
qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
qy = y * y * (2.d0 + y) / (1.d0 + y + y * y)**2
s1 = 1.d0 + 2.d0 * al / be * t * t * xy
h0 = fz3 * be * be / (2.d0 * al) * log (s1)
h0 = fz3 * be * be / (2.d0 * al) * log(s1)
dh0up = be * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
(af * bfup / be-7.d0 / 3.d0) )
dh0dw = be * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
(af * bfdw / be-7.d0 / 3.d0) )
dh0zup = (3.d0 * h0 / fz - be * t * t * fz2 / s1 * (2.d0 * xy - &
qy * (3.d0 * af * expe * ec(1) / fz3 / be+2.d0) ) ) * dfz * (1.d0 - &
qy * (3.d0 * af * expe * ec / fz3 / be+2.d0) ) ) * dfz * (1.d0 - &
zeta)
dh0zdw = - (3.d0 * h0 / fz - be * t * t * fz3 / s1 * (2.d0 * xy - &
qy * (3.d0 * af * expe * ec(1) / fz3 / be+2.d0) ) ) * dfz * (1.d0 + &
qy * (3.d0 * af * expe * ec / fz3 / be+2.d0) ) ) * dfz * (1.d0 + &
zeta)
ddh0 = be * fz / (2.d0 * ks * ks * rho) * (xy - qy) / s1
ee = - 100.d0 * fz4 * (ks / kf * t) **2
cna = cxc0 + pa * rs(1) + pb * rs2
dcna = pa * rs(1) + 2.d0 * pb * rs2
cnb = 1.d0 + pc * rs(1) + pd * rs2 + 1.d4 * pb * rs3
dcnb = pc * rs(1) + 2.d0 * pd * rs2 + 3.d4 * pb * rs3
cna = cxc0 + pa * rs + pb * rs2
dcna = pa * rs + 2.d0 * pb * rs2
cnb = 1.d0 + pc * rs + pd * rs2 + 1.d4 * pb * rs3
dcnb = pc * rs + 2.d0 * pd * rs2 + 3.d4 * pb * rs3
cn = cna / cnb - cx
dcn = dcna / cnb - cna * dcnb / (cnb * cnb)
h1 = nu * (cn - cc0 - 3.d0 / 7.d0 * cx) * fz3 * t * t * exp (ee)
h1 = nu * (cn - cc0 - 3.d0 / 7.d0 * cx) * fz3 * t * t * exp(ee)
dh1 = - third * (h1 * (7.d0 + 8.d0 * ee) + fz3 * nu * t * t * exp &
(ee) * dcn)
ddh1 = 2.d0 * h1 * (1.d0 + ee) * rho / grho
@ -192,9 +189,8 @@ subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
integer, intent(in) :: iflag
real(DP) :: rho, zeta, grho, sc, v1cup, v1cdw, v2c
!
integer, parameter :: length=1 !^^^PROVISIONAL
real(dp), dimension(length) :: rs , zeta_v, ec
real(dp), dimension(length,2) :: vc
real(dp) :: rs, ec
real(dp) :: vc(2)
!
real(DP) :: ga, be(2)
parameter (ga = 0.031091d0)
@ -208,12 +204,11 @@ subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
real(DP) :: fz, fz2, fz3, fz4, dfz, bfup, bfdw, dh0up, dh0dw, &
dh0zup, dh0zdw
!
rs(1) = pi34 / rho**third
zeta_v(1) = zeta
rs = pi34 / rho**third
!
call pw_spin( length, rs, zeta_v, ec, vc )
call pw_spin( rs, zeta, ec, vc )
!
kf = xkf / rs(1)
kf = xkf / rs
ks = xks * sqrt(kf)
fz = 0.5d0 * ( (1.d0 + zeta) ** (2.d0 / 3.d0) + (1.d0 - zeta) ** ( &
2.d0 / 3.d0) )
@ -223,10 +218,10 @@ subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
dfz = ( (1.d0 + zeta) ** ( - 1.d0 / 3.d0) - (1.d0 - zeta) ** ( - &
1.d0 / 3.d0) ) / 3.d0
t = sqrt(grho) / (2.d0 * fz * ks * rho)
expe = exp( - ec(1) / (fz3 * ga) )
expe = exp( - ec / (fz3 * ga) )
af = be(iflag) / ga * (1.d0 / (expe-1.d0) )
bfup = expe * (vc(1,1) - ec(1)) / fz3
bfdw = expe * (vc(1,2) - ec(1)) / fz3
bfup = expe * (vc(1) - ec) / fz3
bfdw = expe * (vc(2) - ec) / fz3
y = af * t * t
xy = (1.d0 + y) / (1.d0 + y + y * y)
qy = y * y * (2.d0 + y) / (1.d0 + y + y * y) **2
@ -237,10 +232,10 @@ subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
dh0dw = be(iflag) * t * t * fz3 / s1 * ( - 7.d0 / 3.d0 * xy - qy * &
(af * bfdw / be(iflag)-7.d0 / 3.d0) )
dh0zup = (3.d0 * h0 / fz - be(iflag) * t * t * fz2 / s1 * (2.d0 * xy - &
qy * (3.d0 * af * expe * ec(1) / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 - zeta)
qy * (3.d0 * af * expe * ec / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 - zeta)
dh0zdw = - (3.d0 * h0 / fz - be(iflag) * t * t * fz2 / s1 * (2.d0 * xy - &
qy * (3.d0 * af * expe * ec(1) / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 + zeta)
qy * (3.d0 * af * expe * ec / fz3 / be(iflag)+2.d0) ) ) * dfz * (1.d0 + zeta)
!
ddh0 = be(iflag) * fz / (2.d0 * ks * ks * rho) * (xy - qy) / s1
sc = rho * h0
v1cup = h0 + dh0up + dh0zup
@ -249,44 +244,3 @@ subroutine pbec_spin (rho, zeta, grho, iflag, sc, v1cup, v1cdw, v2c)
return
end subroutine pbec_spin
!
!-----------------------------------------------------------------------
function dpz_polarized (rs, iflg)
!-----------------------------------------------------------------------
! derivative of the correlation potential with respect to local density
! Perdew and Zunger parameterization of the Ceperley-Alder functional
! spin-polarized case
!
USE kinds, only : DP
USE constants, ONLY : pi, fpi
!
implicit none
!
real(DP), intent (in) :: rs
integer, intent(in) :: iflg
real(DP) :: dpz_polarized
!
! local variables
! a,b,c,d,gc,b1,b2 are the parameters defining the functional
!
real(DP), parameter :: a = 0.01555d0, b = -0.0269d0, c = 0.0007d0, &
d = -0.0048d0, gc = -0.0843d0, b1 = 1.3981d0, b2 = 0.2611d0,&
a1 = 7.0d0 * b1 / 6.d0, a2 = 4.d0 * b2 / 3.d0
real(DP) :: x, den, dmx, dmrs
!
!
if (iflg == 1) then
dmrs = a / rs + 2.d0 / 3.d0 * c * (log (rs) + 1.d0) + &
(2.d0 * d-c) / 3.d0
else
x = sqrt (rs)
den = 1.d0 + x * (b1 + x * b2)
dmx = gc * ( (a1 + 2.d0 * a2 * x) * den - 2.d0 * (b1 + 2.d0 * &
b2 * x) * (1.d0 + x * (a1 + x * a2) ) ) / den**3
dmrs = 0.5d0 * dmx / x
endif
!
dpz_polarized = - fpi * rs**4.d0 / 9.d0 * dmrs
return
!
end function dpz_polarized

114
Modules/metagga.f90
View File

@ -119,8 +119,7 @@ subroutine metax(rho,grho2,tau,ex,v1x,v2x,v3x)
! OUTPUT
real(DP) :: ex,v1x,v2x,v3x
! LOCAL
integer, parameter :: length=1 !^^^ PROVISIONAL
real(DP), dimension(length) :: rs, vx_unif, ex_unif
real(DP) :: rs, vx_unif, ex_unif
! ex_unif: lda \epsilon_x(rho)
! ec_unif: lda \epsilon_c(rho)
real(DP) :: small, pi34, third
@ -139,11 +138,11 @@ subroutine metax(rho,grho2,tau,ex,v1x,v2x,v3x)
v3x=0.0_DP
return
endif
rs(1) = pi34/rho**third
call slater( length, rs, ex_unif, vx_unif )
rs = pi34/rho**third
call slater( rs, ex_unif, vx_unif )
call metaFX(rho,grho2,tau,fx,f1x,f2x,f3x)
ex =rho*ex_unif(1)
v1x=vx_unif(1)*fx + ex*f1x
ex =rho*ex_unif
v1x=vx_unif*fx + ex*f1x
v2x=ex*f2x
v3x=ex*f3x
ex =ex*fx
@ -172,10 +171,9 @@ subroutine metac(rho,grho2,tau,ec,v1c,v2c,v3c)
real(DP) :: v1c_pbe, v2c_pbe, ec_pbe
real(DP) :: v1c_sum, v2c_sum, ec_sum
!
integer, parameter :: length=1 !^^^ PROVISIONAL
real(DP), dimension(length) :: rs, zeta
real(DP), dimension(length) :: ec_unif, vc_unif
real(DP), dimension(length,2) :: vc_unif_s
real(DP) :: rs, zeta
real(DP) :: ec_unif, vc_unif
real(DP) :: vc_unif_s(2)
!
real(DP) :: dd,cab,cabone
real(DP) :: rhoup,grhoup,dummy
@ -197,10 +195,9 @@ subroutine metac(rho,grho2,tau,ec,v1c,v2c,v3c)
grhoup=0.5_DP*SQRT(grho2)
if(rhoup.gt.small) then
!
rs(1)=(pi34/rhoup)**third
zeta(1)=1.0_DP
call pw_spin( length, rs, zeta, ec_unif, vc_unif_s )
dummy = vc_unif_s(1,2)
rs=(pi34/rhoup)**third
call pw_spin( rs, 1.0_DP, ec_unif, vc_unif_s )
dummy = vc_unif_s(2)
!
if(abs(grhoup).gt.small) then
!1.0_DP-small to avoid pow_e of 0 in pbec_spin
@ -211,8 +208,8 @@ subroutine metac(rho,grho2,tau,ec,v1c,v2c,v3c)
v1c_sum=0.0_DP
v2c_sum=0.0_DP
endif
ec_sum = ec_sum/rhoup + ec_unif(1)
v1c_sum = (v1c_sum + vc_unif_s(1,1)-ec_sum)/rho !rho, not rhoup
ec_sum = ec_sum/rhoup + ec_unif
v1c_sum = (v1c_sum + vc_unif_s(1)-ec_sum)/rho !rho, not rhoup
v2c_sum = v2c_sum/(2.0_DP*rho)
else
ec_sum=0.0_DP
@ -220,16 +217,16 @@ subroutine metac(rho,grho2,tau,ec,v1c,v2c,v3c)
v2c_sum=0.0_DP
endif
!
rs(1) = (pi34/rho)**third
call pw( length, rs, 1, ec_unif, vc_unif )
rs = (pi34/rho)**third
call pw( rs, 1, ec_unif, vc_unif )
!
! PBE correlation energy and potential
! ec_pbe=rho*H, not rho*(epsion_c_uinf + H)
! v1c_pbe=D (rho*H) /D rho
! v2c_pbe= for rho, 2 for
call pbec(rho,grho2,1,ec_pbe,v1c_pbe,v2c_pbe)
ec_pbe=ec_pbe/rho+ec_unif(1)
v1c_pbe=(v1c_pbe+vc_unif(1)-ec_pbe)/rho
ec_pbe=ec_pbe/rho+ec_unif
v1c_pbe=(v1c_pbe+vc_unif-ec_pbe)/rho
v2c_pbe=v2c_pbe/rho
!
if(ec_sum .lt. ec_pbe) then
@ -499,9 +496,8 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
real(DP) :: rhoup, rhodw,tauw,grhovec(3),grho2,grho,&
grhoup2,grhodw2
!
integer, parameter :: length=1 !^^^ PROVISIONAL
real(DP), dimension(length) :: rs, zeta_v, ec_u
real(DP), dimension(length,2) :: vc_u
real(DP) :: rs, zeta_v, ec_u
real(DP) :: vc_u(2)
!
!grhovec vector gradient of rho
!grho mod of gradient of rho
@ -545,9 +541,8 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
if(rho > small) then
v2_tmp=0.0_DP
!
rs(1) = (pi34/rho)**third
zeta_v(1) = zeta
call pw_spin( length, rs, zeta_v, ec_u, vc_u )
rs = (pi34/rho)**third
call pw_spin( rs, zeta, ec_u, vc_u )
!
if((abs(grho) > small) .and. (zeta <= 1.0_DP)) then
call pbec_spin(rho,zeta,grho2,1,&
@ -558,10 +553,10 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
v1dw_pbe=0.0_DP
v2up_pbe=0.0_DP
endif
ec_pbe = ec_pbe/rho+ec_u(1)
ec_pbe = ec_pbe/rho+ec_u
! v1xx_pbe = D_epsilon_c/ D_rho_xx :xx= up, dw
v1up_pbe = (v1up_pbe+vc_u(1,1)-ec_pbe)/rho
v1dw_pbe = (v1dw_pbe+vc_u(1,2)-ec_pbe)/rho
v1up_pbe = (v1up_pbe+vc_u(1)-ec_pbe)/rho
v1dw_pbe = (v1dw_pbe+vc_u(2)-ec_pbe)/rho
! v2xx_pbe = (D_Ec / D grho)/rho = (D_Ec/ D |grho| /|grho|)*grho/rho
v2up_pbe = v2_tmp/rho*grhovec
! v2dw === v2up for PBE
@ -578,8 +573,7 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
v2_tmp=0.0_DP
!
rs = (pi34/rhoup)**third
zeta_v(1) = 1.0_DP
call pw_spin( length, rs, zeta_v, ec_u, vc_u )
call pw_spin( rs, 1.0_DP, ec_u, vc_u )
!
if(sqrt(grhoup2) > small) then
call pbec_spin(rhoup,1.0_DP-small,grhoup2,1,&
@ -589,8 +583,8 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
v1up_0=0.0_DP
v2up_0=0.0_DP
endif
ecup_0 = ecup_0/rhoup + ec_u(1)
v1up_0 = (v1up_0 + vc_u(1,1)-ecup_0)/rhoup
ecup_0 = ecup_0/rhoup + ec_u
v1up_0 = (v1up_0 + vc_u(1)-ecup_0)/rhoup
v2up_0 = v2_tmp/rhoup*grhoup
else
ecup_0 = 0.0_DP
@ -616,9 +610,8 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
if(rhodw > small) then
v2_tmp=0.0_DP
!
rs(1) = (pi34/rhodw)**third
zeta_v(1) = -1.0_DP
call pw_spin( length, rs, zeta_v, ec_u, vc_u )
rs = (pi34/rhodw)**third
call pw_spin( rs, -1.0_DP, ec_u, vc_u )
!
if(sqrt(grhodw2) > small) then
call pbec_spin(rhodw,-1.0_DP+small,grhodw2,1,&
@ -628,8 +621,8 @@ subroutine metac_spin(rho,zeta,grhoup,grhodw, &
v1dw_0=0.0_DP
v2dw_0=0.0_DP
endif
ecdw_0 = ecdw_0/rhodw + ec_u(1)
v1dw_0 = (v1dw_0 + vc_u(1,2)-ecdw_0)/rhodw
ecdw_0 = ecdw_0/rhodw + ec_u
v1dw_0 = (v1dw_0 + vc_u(2)-ecdw_0)/rhodw
v2dw_0 = v2_tmp/rhodw*grhodw
else
ecdw_0 = 0.0_DP
@ -1110,12 +1103,11 @@ subroutine m06lc(rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c_up, v2c_up, v3c
implicit none
!-------------------------------------------------------------------------
real(dp), intent(in) :: rhoa, rhob, grho2a, grho2b, taua, taub
real(dp), intent(out) :: ec, v1c_up, v2c_up, v3c_up, v1c_dw, v2c_dw, v3c_dw
real(dp), intent(in) :: rhoa, rhob, grho2a, grho2b, taua, taub
real(dp), intent(out) :: ec, v1c_up, v2c_up, v3c_up, v1c_dw, v2c_dw, v3c_dw
!
integer, parameter :: length=1 !^^^ PROVISIONAL (xc-lib)
real(DP), dimension(length) :: rs, zeta
real(DP), dimension(length,2) :: vc_v
real(DP) :: rs, zeta
real(DP) :: vc_v(2)
!
!
real(dp), parameter :: zero = 0._dp, one = 1.0_dp, two=2.0_dp, three = 3.0_dp, &
@ -1136,7 +1128,7 @@ subroutine m06lc(rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c_up, v2c_up, v3c
!
! functions and variables
!
real(dp),dimension(length) :: ec_pw_a, ec_pw_b, ec_pw_ab
real(dp) :: ec_pw_a, ec_pw_b, ec_pw_ab
!
real(dp) :: vv, vc_pw_a, vc_pw_b, vc_pw_ab, vc_pw_up, vc_pw_dw, Ecaa, Ecbb, Ecab, &
& Ec_UEG_ab, Ec_UEG_aa, Ec_UEG_bb, decab_drhoa, decab_drhob, &
@ -1224,16 +1216,16 @@ subroutine m06lc(rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c_up, v2c_up, v3c
ec_pw_a = zero
vc_pw_a = zero
!
rs(1) = rsa
zeta(1) = one
call pw_spin( length, rs, zeta, ec_pw_a, vc_v )
vc_pw_a = vc_v(1,1)
vv = vc_v(1,2)
rs = rsa
zeta = one
call pw_spin( rs, zeta, ec_pw_a, vc_v )
vc_pw_a = vc_v(1)
vv = vc_v(2)
!
call gvt4 (xs2a, zsa, ds0, ds1, ds2, ds3, ds4, ds5, alpha_s, hga, dhga_dxs2a, dhga_dzsa)
call gfunc (cs, gama_s, xs2a, gsa, dgsa_dxs2a)
!
Ec_UEG_aa = rhoa*ec_pw_a(1)
Ec_UEG_aa = rhoa*ec_pw_a
num = (dgsa_dxs2a + dhga_dxs2a)*Dsa + (gsa + hga)*dDsa_dxs2a
!
!
@ -1276,16 +1268,16 @@ subroutine m06lc(rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c_up, v2c_up, v3c
dDsb_dxs2b = - one/(four * (zsb + CF))
dDsb_dzsb = xs2b/(four * (zsb + CF)**2)
!
rs(1) = rsb
zeta(1) = one
call pw_spin( length, rs, zeta, ec_pw_b, vc_v )
vc_pw_b = vc_v(1,1)
vv = vc_v(1,2)
zeta = one
rs = rsb
call pw_spin( rs, zeta, ec_pw_b, vc_v )
vc_pw_b = vc_v(1)
vv = vc_v(2)
!
call gvt4 (xs2b, zsb, ds0, ds1, ds2, ds3, ds4, ds5, alpha_s, hgb, dhgb_dxs2b, dhgb_dzsb)
call gfunc (cs, gama_s, xs2b, gsb, dgsb_dxs2b)
Ec_UEG_bb = rhob*ec_pw_b(1)
Ec_UEG_bb = rhob*ec_pw_b
num = (dgsb_dxs2b + dhgb_dxs2b)*Dsb + (gsb + hgb)*dDsb_dxs2b
!
!
@ -1318,20 +1310,20 @@ subroutine m06lc(rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c_up, v2c_up, v3c
xs2ab = xs2a + xs2b
zsab = zsa + zsb
rho = rhoa + rhob
zeta(1) = (rhoa - rhob)/rho
rs(1) = (pi34/rho)**f13
zeta = (rhoa - rhob)/rho
rs = (pi34/rho)**f13
call gvt4(xs2ab, zsab, dab0, dab1, dab2, dab3, dab4, dab5, alpha_ab, hgab, dhgab_dxs2ab, dhgab_dzsab)
!
call pw_spin( length, rs, zeta, ec_pw_ab, vc_v )
vc_pw_up = vc_v(1,1) ; vc_pw_dw=vc_v(1,2)
call pw_spin( rs, zeta, ec_pw_ab, vc_v )
vc_pw_up = vc_v(1) ; vc_pw_dw=vc_v(2)
!
call gfunc(cab, gama_ab, xs2ab, gsab, dgsab_dxs2ab)
!
decab_drhoa = vc_pw_up - vc_pw_a
decab_drhob = vc_pw_dw - vc_pw_b
Ec_UEG_ab = ec_pw_ab(1)*rho - ec_pw_a(1)*rhoa - ec_pw_b(1)*rhob
Ec_UEG_ab = ec_pw_ab*rho - ec_pw_a*rhoa - ec_pw_b*rhob
!
!
Ecab = Ec_UEG_ab * (gsab + hgab)

869
Modules/xc_lda_lsda_drivers.f90
View File

@ -16,12 +16,16 @@ PRIVATE
SAVE
!
! LDA and LSDA exchange-correlation drivers
PUBLIC :: xc, xc_spin, select_functionals
PUBLIC :: xc, xc_lda, xc_lsda, select_lda_functionals
!
PUBLIC :: libxc_switches
PUBLIC :: iexch_l, icorr_l
PUBLIC :: exx_started_l, exx_fraction_l
PUBLIC :: is_there_finite_size_corr, finite_size_cell_volume_l
!
! use qe or libxc for the different terms (0: qe, 1: libxc)
INTEGER :: libxc_switches(2)
!
! indexes defining xc functionals
INTEGER :: iexch_l, icorr_l
!
@ -36,7 +40,7 @@ REAL(DP) :: exx_fraction_l, finite_size_cell_volume_l
!----------------------------------------------------------------------------
!----- Select functionals by the corresponding indexes ----------------------
!----------------------------------------------------------------------------
SUBROUTINE select_functionals( iexch, icorr, exx_fraction, finite_size_cell_volume )
SUBROUTINE select_lda_functionals( iexch, icorr, exx_fraction, finite_size_cell_volume )
!
IMPLICIT NONE
!
@ -48,409 +52,632 @@ SUBROUTINE select_functionals( iexch, icorr, exx_fraction, finite_size_cell_volu
icorr_l = icorr
!
! hybrid exchange vars
exx_started_l = .false.
exx_started_l = .FALSE.
exx_fraction_l = 0._DP
IF ( present(exx_fraction) ) THEN
exx_started_l = .true.
IF ( PRESENT(exx_fraction) ) THEN
exx_started_l = .TRUE.
exx_fraction_l = exx_fraction
ENDIF
!
! finite size correction vars
is_there_finite_size_corr = .false.
is_there_finite_size_corr = .FALSE.
finite_size_cell_volume_l = -1.0_DP
IF ( present(finite_size_cell_volume) ) THEN
is_there_finite_size_corr = .true.
IF ( PRESENT(finite_size_cell_volume) ) THEN
is_there_finite_size_corr = .TRUE.
finite_size_cell_volume_l = finite_size_cell_volume
ENDIF
!
!
RETURN
!
END SUBROUTINE
END SUBROUTINE select_lda_functionals
!
!
!----------------------------------------------------------------------------
!------- LDA-LSDA DRIVERS --------------------------------------------------
!----------------------------------------------------------------------------
!
!----------------------------------------------------------------------------
SUBROUTINE xc( length, rho, ex, ec, vx, vc )
!--------------------------------------------------------------------------
! lda exchange and correlation functionals - Hartree a.u.
!---------------------------------------------------------------------
SUBROUTINE xc( length, sr_d, sv_d, rho_in, ex_out, ec_out, vx_out, vc_out )
!---------------------------------------------------------------------
!! Wrapper routine. Calls xc-driver routines from internal libraries
!! of q-e or from the external libxc, depending on the input choice.
!
! exchange : Slater, relativistic Slater
! correlation: Ceperley-Alder (Perdew-Zunger parameters)
! Vosko-Wilk-Nusair
! Lee-Yang-Parr
! Perdew-Wang
! Wigner
! Hedin-Lundqvist
! Ortiz-Ballone (Perdew-Zunger formula)
! Ortiz-Ballone (Perdew-Wang formula)
! Gunnarsson-Lundqvist
!
! input : rho=rho(r)
! definitions: E_x = \int E_x(rho) dr, E_x(rho) = rho\epsilon_c(rho)
! same for correlation
! output: ex = \epsilon_x(rho) ( NOT E_x(rho) )
! vx = dE_x(rho)/drho ( NOT d\epsilon_x(rho)/drho )
! ec, vc as above for correlation
!! NOTE: look at 'PP/src/benchmark_libxc.f90' to see the differences
!! between q-e and libxc.
!
#if defined(__LIBXC)
USE xc_f90_types_m
USE xc_f90_lib_m
#endif
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN), DIMENSION(length) :: rho
REAL(DP), INTENT(OUT), DIMENSION(length) :: ec, vc, ex, vx
!! length of the I/O arrays
INTEGER, INTENT(IN) :: sr_d
!! spin dimension of rho
INTEGER, INTENT(IN) :: sv_d
!! spin dimension of v
REAL(DP), INTENT(IN) :: rho_in(length,sr_d)
!! Charge density
REAL(DP), INTENT(OUT) :: ex_out(length)
!! \(\epsilon_x(rho)\) ( NOT \(E_x(\text{rho})\) )
REAL(DP), INTENT(OUT) :: vx_out(length,sv_d)
!! \(dE_x(\text{rho})/d\text{rho} ( NOT d\epsilon_x(\text{rho})/d\text{rho} )
REAL(DP), INTENT(OUT) :: ec_out(length)
!! \(\epsilon_c(rho)\) ( NOT \(E_c(\text{rho})\) )
REAL(DP), INTENT(OUT) :: vc_out(length,sv_d)
!! \(dE_c(\text{rho})/d\text{rho} ( NOT d\epsilon_c(\text{rho})/d\text{rho} )
!
! local vars
! ... local variables
!
REAL(DP), ALLOCATABLE, DIMENSION(:) :: ec_, vc_
REAL(DP), PARAMETER :: small = 1.E-10_DP, third = 1.0_DP / 3.0_DP, &
pi34 = 0.6203504908994_DP ! pi34=(3/4pi)^(1/3)
REAL(DP) :: rs(length)
#if defined(__LIBXC)
TYPE(xc_f90_pointer_t) :: xc_func
TYPE(xc_f90_pointer_t) :: xc_info1, xc_info2
REAL(DP) :: amag
REAL(DP), ALLOCATABLE :: rho_lxc(:)
REAL(DP), ALLOCATABLE :: vx_lxc(:), vc_lxc(:)
#endif
!
REAL(DP), ALLOCATABLE :: arho(:), zeta(:)
!
INTEGER :: ir
REAL(DP), PARAMETER :: small = 1.E-10_DP
!
!
WHERE ( rho > small )
rs = pi34 / rho**third
ELSEWHERE
rs = 1.0_DP
END WHERE
#if defined(__LIBXC)
!
IF ( icorr_l >= 12 .AND. icorr_l <= 14 ) THEN
ALLOCATE( ec_(length) )
ALLOCATE( vc_(length) )
IF (SUM(libxc_switches) /= 0) THEN
!
ALLOCATE( rho_lxc(length*sv_d) )
ALLOCATE( vx_lxc(length*sv_d), vc_lxc(length*sv_d) )
!
! ... set libxc input
SELECT CASE( sr_d )
CASE( 1 )
!
rho_lxc(:) = rho_in(:,1)
!
CASE( 2 )
!
DO ir = 1, length
rho_lxc(2*ir-1) = (rho_in(ir,1) + rho_in(ir,2)) * 0.5_DP
rho_lxc(2*ir) = (rho_in(ir,1) - rho_in(ir,2)) * 0.5_DP
ENDDO
!
CASE( 4 )
!
DO ir = 1, length
amag = SQRT( SUM(rho_in(ir,2:4)**2) )
rho_lxc(2*ir-1) = (rho_in(ir,1) + amag) * 0.5_DP
rho_lxc(2*ir) = (rho_in(ir,1) - amag) * 0.5_DP
ENDDO
!
CASE DEFAULT
!
CALL errore( 'xc_LDA', 'Wrong number of spin dimensions', 1 )
!
END SELECT
!
ENDIF
!
!
! ... EXCHANGE
!
SELECT CASE( iexch_l )
CASE( 1 ) ! 'sla'
!
CALL slater( length, rs, ex, vx )
!
CASE( 2 ) ! 'sl1'
!
CALL slater1( length, rs, ex, vx )
!
CASE( 3 ) ! 'rxc'
!
CALL slater_rxc( length, rs, ex, vx )
!
CASE( 4, 5 ) ! 'oep','hf'
!
IF ( exx_started_l ) THEN
ex = 0.0_DP
vx = 0.0_DP
ELSE
CALL slater( length, rs, ex, vx )
ENDIF
!
CASE( 6, 7 ) ! 'pb0x' or 'DF-cx-0', or 'DF2-0',
! ! 'B3LYP'
CALL slater( length, rs, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE( 8 ) ! 'sla+kzk'
!
IF (.NOT. is_there_finite_size_corr) CALL errore ('XC',&
'finite size corrected exchange used w/o initialization',1)
CALL slaterKZK( length, rs, ex, vx, finite_size_cell_volume_l )
!
CASE( 9 ) ! 'X3LYP'
!
CALL slater( length, rs, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE DEFAULT
!
ex = 0.0_DP
vx = 0.0_DP
!
END SELECT
!
IF ( libxc_switches(1)==1 ) THEN
CALL xc_f90_func_init( xc_func, xc_info1, iexch_l, sv_d )
CALL xc_f90_lda_exc_vxc( xc_func, length, rho_lxc(1), ex_out(1), vx_lxc(1) )
CALL xc_f90_func_end( xc_func )
ENDIF
!
! ... CORRELATION
IF ( libxc_switches(2)==1 ) THEN
CALL xc_f90_func_init( xc_func, xc_info2, icorr_l, sv_d )
CALL xc_f90_lda_exc_vxc( xc_func, length, rho_lxc(1), ec_out(1), vc_lxc(1) )
CALL xc_f90_func_end( xc_func )
ENDIF
!
SELECT CASE( icorr_l )
!
IF ( (libxc_switches(1)==0) .OR. (libxc_switches(2)==0) ) THEN
!
SELECT CASE( sr_d )
CASE( 1 )
!
CALL xc_lda( length, ABS(rho_in(:,1)), ex_out, ec_out, vx_out(:,1), vc_out(:,1) )
!
CASE( 2 )
!
ALLOCATE( arho(length), zeta(length) )
arho = ABS(rho_in(:,1))
WHERE (arho > small) zeta(:) = rho_in(:,2) / arho(:)
CALL xc_lsda( length, arho, zeta, ex_out, ec_out, vx_out, vc_out )
DEALLOCATE( arho, zeta )
!
CASE( 4 )
!
ALLOCATE( arho(length), zeta(length) )
arho = ABS( rho_in(:,1) )
WHERE (arho > small) zeta(:) = SQRT( rho_in(:,2)**2 + rho_in(:,3)**2 + &
rho_in(:,4)**2 ) / arho(:) ! amag/arho
CALL xc_lsda( length, arho, zeta, ex_out, ec_out, vx_out, vc_out )
DEALLOCATE( arho, zeta )
!
CASE DEFAULT
!
CALL errore( 'xc_LDA', 'Wrong ns input', 2 )
!
END SELECT
!
ENDIF
!
! ... fill output arrays
!
IF (sv_d == 1) THEN
IF (libxc_switches(1)==1) vx_out(:,1) = vx_lxc(:)
IF (libxc_switches(2)==1) vc_out(:,1) = vc_lxc(:)
ELSE
IF (libxc_switches(1)==1) THEN
DO ir = 1, length
vx_out(ir,1) = vx_lxc(2*ir-1)
vx_out(ir,2) = vx_lxc(2*ir)
ENDDO
ENDIF
IF (libxc_switches(2)==1) THEN
DO ir = 1, length
vc_out(ir,1) = vc_lxc(2*ir-1)
vc_out(ir,2) = vc_lxc(2*ir)
ENDDO
ENDIF
ENDIF
!
IF (SUM(libxc_switches) /= 0) THEN
DEALLOCATE( rho_lxc )
DEALLOCATE( vx_lxc, vc_lxc )
ENDIF
!
#else
!
SELECT CASE( sr_d )
CASE( 1 )
!
CALL pz( length, rs, 1, ec, vc )
CALL xc_lda( length, ABS(rho_in(:,1)), ex_out, ec_out, vx_out(:,1), vc_out(:,1) )
!
CASE( 2 )
!
CALL vwn( length, rs, ec, vc )
ALLOCATE( arho(length), zeta(length) )
!
CASE( 3 )
arho = ABS(rho_in(:,1))
WHERE (arho > small) zeta(:) = rho_in(:,2) / arho(:)
!
CALL lyp( length, rs, ec, vc )
CALL xc_lsda( length, arho, zeta, ex_out, ec_out, vx_out, vc_out )
!
DEALLOCATE( arho, zeta )
!
CASE( 4 )
!
CALL pw( length, rs, 1, ec, vc )
ALLOCATE( arho(length), zeta(length) )
!
CASE( 5 )
arho = ABS( rho_in(:,1) )
WHERE (arho > small) zeta(:) = SQRT( rho_in(:,2)**2 + rho_in(:,3)**2 + &
rho_in(:,4)**2 ) / arho(:) ! amag/arho
!
CALL wignerc( length, rs, ec, vc )
CALL xc_lsda( length, arho, zeta, ex_out, ec_out, vx_out, vc_out )
!
CASE( 6 )
!
CALL hl( length, rs, ec, vc )
!
CASE( 7 )
!
CALL pz( length, rs, 2, ec, vc )
!
CASE( 8 )
!
CALL pw( length, rs, 2, ec, vc )
!
CASE( 9 )
!
CALL gl( length, rs, ec, vc )
!
CASE( 10 )
!
IF (.NOT. is_there_finite_size_corr) CALL errore ('XC',&
'finite size corrected exchange used w/o initialization',2)
CALL pzKZK ( length, rs, ec, vc, finite_size_cell_volume_l )
!
CASE( 11 )
!
CALL vwn1_rpa( length, rs, ec, vc )
!
CASE( 12 ) ! 'B3LYP'
!
CALL vwn( length, rs, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lyp( length, rs, ec_, vc_ )
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 13 ) ! 'B3LYP-V1R'
!
CALL vwn1_rpa ( length, rs, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lyp( length, rs, ec_, vc_ )
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 14 ) ! 'X3LYP'
!
CALL vwn1_rpa( length, rs, ec, vc )
ec = 0.129_DP * ec
vc = 0.129_DP * vc
!
CALL lyp( length, rs, ec_, vc_ )
ec = ec + 0.871_DP * ec_
vc = vc + 0.871_DP * vc_
DEALLOCATE( arho, zeta )
!
CASE DEFAULT
!
ec = 0.0_DP
vc = 0.0_DP
CALL errore( 'xc_LDA', 'Wrong ns input', 2 )
!
END SELECT
!
!
WHERE ( rho <= small )
ex = 0.0_DP
ec = 0.0_DP
vx = 0.0_DP
vc = 0.0_DP
END WHERE
!
IF ( icorr_l >= 12 .AND. icorr_l <= 14 ) THEN
DEALLOCATE( ec_ )
DEALLOCATE( vc_ )
ENDIF
#endif
!
!
RETURN
!
END SUBROUTINE xc
!
!----------------------------------------------------------------------------
!------- LDA-LSDA DRIVERS --------------------------------------------------
!----------------------------------------------------------------------------
!
!
SUBROUTINE xc_spin( length, rho, zeta, ex, ec, vx, vc )
!---------------------------------------------------------------------
! lsd exchange and correlation functionals - Hartree a.u.
!----------------------------------------------------------------------------
SUBROUTINE xc_lda( length, rho_in, ex_out, ec_out, vx_out, vc_out )
!--------------------------------------------------------------------------
!! LDA exchange and correlation functionals - Hartree a.u.
!
! exchange : Slater (alpha=2/3)
! correlation: Ceperley & Alder (Perdew-Zunger parameters)
! Perdew & Wang
!! * Exchange:
!! * Slater;
!! * relativistic Slater.
!! * Correlation:
!! * Ceperley-Alder (Perdew-Zunger parameters);
!! * Vosko-Wilk-Nusair;
!! * Lee-Yang-Parr;
!! * Perdew-Wang;
!! * Wigner;
!! * Hedin-Lundqvist;
!! * Ortiz-Ballone (Perdew-Zunger formula);
!! * Ortiz-Ballone (Perdew-Wang formula);
!! * Gunnarsson-Lundqvist.
!
! input : rho = rhoup(r)+rhodw(r)
! zeta=(rhoup(r)-rhodw(r))/rho
!! NOTE:
!! $$ E_x = \int E_x(\text{rho}) dr, E_x(\text{rho}) =
!! \text{rho}\epsilon_c(\text{rho})\ . $$
!! Same for correlation.
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
REAL(DP), INTENT(IN), DIMENSION(length) :: rho, zeta
REAL(DP), INTENT(OUT), DIMENSION(length) :: ex, ec !^^^
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: vx, vc
!! length of the I/O arrays
REAL(DP), INTENT(IN), DIMENSION(length) :: rho_in
!! Charge density
REAL(DP), INTENT(OUT), DIMENSION(length) :: ex_out
!! \(\epsilon_x(rho)\) ( NOT \(E_x(\text{rho})\) )
REAL(DP), INTENT(OUT), DIMENSION(length) :: vx_out
!! \(dE_x(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_x(\text{rho})/d\text{rho}\) )
REAL(DP), INTENT(OUT), DIMENSION(length) :: ec_out
!! \(\epsilon_c(rho)\) ( NOT \(E_c(\text{rho})\) )
REAL(DP), INTENT(OUT), DIMENSION(length) :: vc_out
!! \(dE_c(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_c(\text{rho})/d\text{rho}\) )
!
! local vars
REAL(DP), DIMENSION(length) :: rs
REAL(DP), ALLOCATABLE, DIMENSION(:) :: ec_
REAL(DP), ALLOCATABLE, DIMENSION(:,:) :: vc_
! ... local variables
!
REAL(DP), PARAMETER :: small= 1.E-10_DP, third = 1.0_DP/3.0_DP, &
pi34= 0.6203504908994_DP ! pi34=(3/4pi)^(1/3)
!
WHERE ( rho > small )
rs = pi34 / rho**third
ELSEWHERE
rs = 1.0_DP
END WHERE
!
IF ( icorr_l >= 12 .AND. icorr_l <= 14 ) THEN
ALLOCATE( ec_(length) )
ALLOCATE( vc_(length,2) )
ENDIF
INTEGER :: ir
REAL(DP) :: rho, rs
REAL(DP) :: ex, ec, ec_
REAL(DP) :: vx, vc, vc_
REAL(DP), PARAMETER :: small = 1.E-10_DP, third = 1.0_DP/3.0_DP, &
pi34 = 0.6203504908994_DP, e2 = 2.0_DP
! pi34 = (3/4pi)^(1/3)
!
!
! ... EXCHANGE
!
SELECT CASE( iexch_l )
CASE( 1 ) ! 'sla'
!$omp parallel do private( rho, rs, ex, ec, ec_, vx, vc, vc_ )
DO ir = 1, length
!
CALL slater_spin( length, rho, zeta, ex, vx )
rho = ABS(rho_in(ir))
!
CASE( 2 ) ! 'sl1'
! ... RHO THRESHOLD
!
CALL slater1_spin( length, rho, zeta, ex, vx )
IF ( rho > small ) THEN
rs = pi34 / rho**third
ELSE
ex_out(ir) = 0.0_DP ; ec_out(ir) = 0.0_DP
vx_out(ir) = 0.0_DP ; vc_out(ir) = 0.0_DP
CYCLE
ENDIF
!
CASE( 3 ) ! 'rxc'
! ... EXCHANGE
!
CALL slater_rxc_spin( length, rho, zeta, ex, vx )
!
CASE( 4, 5 ) ! 'oep','hf'
!
IF ( exx_started_l ) THEN
SELECT CASE( iexch_l )
CASE( 1 ) ! 'sla'
!
CALL slater( rs, ex, vx )
!
CASE( 2 ) ! 'sl1'
!
CALL slater1( rs, ex, vx )
!
CASE( 3 ) ! 'rxc'
!
CALL slater_rxc( rs, ex, vx )
!
CASE( 4, 5 ) ! 'oep','hf'
!
IF ( exx_started_l ) THEN
ex = 0.0_DP
vx = 0.0_DP
ELSE
CALL slater( rs, ex, vx )
ENDIF
!
CASE( 6, 7 ) ! 'pb0x' or 'DF-cx-0', or 'DF2-0',
! ! 'B3LYP'
CALL slater( rs, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE( 8 ) ! 'sla+kzk'
!
IF (.NOT. is_there_finite_size_corr) CALL errore( 'XC',&
'finite size corrected exchange used w/o initialization', 2 )
CALL slaterKZK( rs, ex, vx, finite_size_cell_volume_l )
!
CASE( 9 ) ! 'X3LYP'
!
CALL slater( rs, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE DEFAULT
!
ex = 0.0_DP
vx = 0.0_DP
ELSE
CALL slater_spin( length, rho, zeta, ex, vx )
ENDIF
!
END SELECT
!
CASE( 6 ) ! 'pb0x'
!
CALL slater_spin( length, rho, zeta, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
! ... CORRELATION
!
CASE( 7 ) ! 'B3LYP'
SELECT CASE( icorr_l )
CASE( 1 )
!
CALL pz( rs, 1, ec, vc )
!
CASE( 2 )
!
CALL vwn( rs, ec, vc )
!
CASE( 3 )
!
CALL lyp( rs, ec, vc )
!
CASE( 4 )
!
CALL pw( rs, 1, ec, vc )
!
CASE( 5 )
!
CALL wignerc( rs, ec, vc )
!
CASE( 6 )
!
CALL hl( rs, ec, vc )
!
CASE( 7 )
!
CALL pz( rs, 2, ec, vc )
!
CASE( 8 )
!
CALL pw( rs, 2, ec, vc )
!
CASE( 9 )
!
CALL gl( rs, ec, vc )
!
CASE( 10 )
!
IF (.NOT. is_there_finite_size_corr) CALL errore( 'XC',&
'finite size corrected exchange used w/o initialization', 3 )
CALL pzKZK( rs, ec, vc, finite_size_cell_volume_l )
!
CASE( 11 )
!
CALL vwn1_rpa( rs, ec, vc )
!
CASE( 12 ) ! 'B3LYP'
!
CALL vwn( rs, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lyp( rs, ec_, vc_ )
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 13 ) ! 'B3LYP-V1R'
!
CALL vwn1_rpa ( rs, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lyp( rs, ec_, vc_ )
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 14 ) ! 'X3LYP'
!
CALL vwn1_rpa( rs, ec, vc )
ec = 0.129_DP * ec
vc = 0.129_DP * vc
!
CALL lyp( rs, ec_, vc_ )
ec = ec + 0.871_DP * ec_
vc = vc + 0.871_DP * vc_
!
CASE DEFAULT
!
ec = 0.0_DP
vc = 0.0_DP
!
END SELECT
!
CALL slater_spin( length, rho, zeta, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
ex_out(ir) = ex ; ec_out(ir) = ec
vx_out(ir) = vx ; vc_out(ir) = vc
!
CASE( 9 ) ! 'X3LYP'
!
CALL slater_spin( length, rho, zeta, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE DEFAULT
!
ex = 0.0_DP
vx = 0.0_DP
!
END SELECT
!
!
! ... CORRELATION
!
SELECT CASE( icorr_l )
CASE( 0 )
!
ec = 0.0_DP
vc = 0.0_DP
!
CASE( 1 )
!
CALL pz_spin( length, rs, zeta, ec, vc )
!
CASE( 2 )
!
CALL vwn_spin( length, rs, zeta, ec, vc )
!
CASE( 3 )
!
CALL lsd_lyp( length, rho, zeta, ec, vc ) ! from CP/FPMD (more_functionals)
!
CASE( 4 )
!
CALL pw_spin( length, rs, zeta, ec, vc )
!
CASE( 12 ) ! 'B3LYP'
!
CALL vwn_spin( length, rs, zeta, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lsd_lyp( length, rho, zeta, ec_, vc_ ) ! from CP/FPMD (more_functionals)
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 13 ) ! 'B3LYP-V1R'
!
CALL vwn1_rpa_spin( length, rs, zeta, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lsd_lyp( length, rho, zeta, ec_, vc_ ) ! from CP/FPMD (more_functionals)
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 14 ) ! 'X3LYP
!
CALL vwn1_rpa_spin( length, rs, zeta, ec, vc )
ec = 0.129_DP * ec
vc = 0.129_DP * vc
!
CALL lsd_lyp( length, rho, zeta, ec_, vc_ ) ! from CP/FPMD (more_functionals)
ec = ec + 0.871_DP * ec_
vc = vc + 0.871_DP * vc_
!
CASE DEFAULT
!
CALL errore( 'lsda_functional (xc_spin)', 'not implemented', icorr_l )
!
END SELECT
!
!
WHERE ( rho <= small )
ex = 0.0_DP
ec = 0.0_DP
vx(:,1) = 0.0_DP
vx(:,2) = 0.0_DP
vc(:,1) = 0.0_DP
vc(:,2) = 0.0_DP
END WHERE
ENDDO
!$omp end parallel do
!
!
RETURN
!
END SUBROUTINE xc_spin
END SUBROUTINE xc_lda
!
!
!-----------------------------------------------------------------------------
SUBROUTINE xc_lsda( length, rho_in, zeta_in, ex_out, ec_out, vx_out, vc_out )
!-----------------------------------------------------------------------------
!! LSD exchange and correlation functionals - Hartree a.u.
!
!! * Exchange:
!! * Slater (alpha=2/3).
!! * Correlation:
!! * Ceperley & Alder (Perdew-Zunger parameters);
!! * Perdew & Wang.
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: length
!! length of the I/O arrays
REAL(DP), INTENT(IN), DIMENSION(length) :: rho_in
!! Total charge density
REAL(DP), INTENT(IN), DIMENSION(length) :: zeta_in
!! zeta = mag / rho_tot
REAL(DP), INTENT(OUT), DIMENSION(length) :: ex_out
!! \(\epsilon_x(rho)\) ( NOT \(E_x(\text{rho})\) )
REAL(DP), INTENT(OUT), DIMENSION(length) :: ec_out
!! \(\epsilon_c(rho)\) ( NOT \(E_c(\text{rho})\) )
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: vx_out
!! \(dE_x(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_x(\text{rho})/d\text{rho}\) )
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: vc_out
!! \(dE_c(\text{rho})/d\text{rho}\) ( NOT \(d\epsilon_c(\text{rho})/d\text{rho}\) )
!
! ... local variables
!
INTEGER :: ir
REAL(DP) :: rho, rs, zeta
REAL(DP) :: ex, ec, ec_
REAL(DP) :: vx(2), vc(2), vc_(2)
!
REAL(DP), PARAMETER :: small= 1.E-10_DP, third = 1.0_DP/3.0_DP, &
pi34 = 0.6203504908994_DP
! pi34 = (3/4pi)^(1/3)
!
!
!$omp parallel do private( rho, rs, zeta, ex, ec, ec_, vx, vc, vc_ )
DO ir = 1, length
!
zeta = zeta_in(ir)
IF (ABS(zeta) > 1.D0) zeta = SIGN( 1.D0, zeta )
!
rho = ABS(rho_in(ir))
!
IF ( rho > small ) THEN
rs = pi34 / rho**third
ELSE
ex_out(ir) = 0.0_DP ; vx_out(ir,:) = 0.0_DP
ec_out(ir) = 0.0_DP ; vc_out(ir,:) = 0.0_DP
CYCLE
ENDIF
!
!
! ... EXCHANGE
!
SELECT CASE( iexch_l )
CASE( 1 ) ! 'sla'
!
CALL slater_spin( rho, zeta, ex, vx )
!
CASE( 2 ) ! 'sl1'
!
CALL slater1_spin( rho, zeta, ex, vx )
!
CASE( 3 ) ! 'rxc'
!
CALL slater_rxc_spin( rho, zeta, ex, vx )
!
CASE( 4, 5 ) ! 'oep','hf'
!
IF ( exx_started_l ) THEN
ex = 0.0_DP
vx = 0.0_DP
ELSE
CALL slater_spin( rho, zeta, ex, vx )
ENDIF
!
CASE( 6 ) ! 'pb0x'
!
CALL slater_spin( rho, zeta, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE( 7 ) ! 'B3LYP'
!
CALL slater_spin( rho, zeta, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE( 9 ) ! 'X3LYP'
!
CALL slater_spin( rho, zeta, ex, vx )
IF ( exx_started_l ) THEN
ex = (1.0_DP - exx_fraction_l) * ex
vx = (1.0_DP - exx_fraction_l) * vx
ENDIF
!
CASE DEFAULT
!
ex = 0.0_DP
vx = 0.0_DP
!
END SELECT
!
!
! ... CORRELATION
!
SELECT CASE( icorr_l )
CASE( 0 )
!
ec = 0.0_DP
vc = 0.0_DP
!
CASE( 1 )
!
CALL pz_spin( rs, zeta, ec, vc )
!
CASE( 2 )
!
CALL vwn_spin( rs, zeta, ec, vc )
!
CASE( 3 )
!
CALL lsd_lyp( rho, zeta, ec, vc ) ! from CP/FPMD (more_functionals)
!
CASE( 4 )
!
CALL pw_spin( rs, zeta, ec, vc )
!
CASE( 12 ) ! 'B3LYP'
!
CALL vwn_spin( rs, zeta, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lsd_lyp( rho, zeta, ec_, vc_ ) ! from CP/FPMD (more_functionals)
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 13 ) ! 'B3LYP-V1R'
!
CALL vwn1_rpa_spin( rs, zeta, ec, vc )
ec = 0.19_DP * ec
vc = 0.19_DP * vc
!
CALL lsd_lyp( rho, zeta, ec_, vc_ ) ! from CP/FPMD (more_functionals)
ec = ec + 0.81_DP * ec_
vc = vc + 0.81_DP * vc_
!
CASE( 14 ) ! 'X3LYP
!
CALL vwn1_rpa_spin( rs, zeta, ec, vc )
ec = 0.129_DP * ec
vc = 0.129_DP * vc
!
CALL lsd_lyp( rho, zeta, ec_, vc_ ) ! from CP/FPMD (more_functionals)
ec = ec + 0.871_DP * ec_
vc = vc + 0.871_DP * vc_
!
CASE DEFAULT
!
CALL errore( 'lsda_functional (xc_spin)', 'not implemented', icorr_l )
!
END SELECT
!
ex_out(ir) = ex ; vx_out(ir,:) = vx(:)
ec_out(ir) = ec ; vc_out(ir,:) = vc(:)
!
ENDDO
!$omp end parallel do
!
!
RETURN
!
END SUBROUTINE xc_lsda
!
!
END MODULE xc_lda_lsda

42
Modules/xc_vdW_DF.f90
View File

@ -594,14 +594,12 @@ CONTAINS
integer, parameter :: m_cut = 12 ! How many terms to include in the sum
! of SOLER equation 5.
integer, parameter :: length = 1 !^^^ PROVISIONAL (xc-lib)
real(dp) :: rho ! Local variable for the density.
real(dp) :: r_s(length) ! Wigner–Seitz radius.
real(dp) :: r_s ! Wigner–Seitz radius.
real(dp) :: s ! Reduced gradient.
real(dp) :: q
real(dp) :: ec(length), dq0dr_v(length) !^^^ PROVISIONAL
real(dp) :: ec
real(dp) :: dq0_dq ! The derivative of the saturated
! q0 with respect to q.
@ -634,20 +632,19 @@ CONTAINS
! -----------------------------------------------------------------
! Calculate some intermediate values needed to find q.
r_s(1) = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
r_s = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
s = sqrt( grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2 ) &
/ (2.0D0 * kF(rho) * rho )
s = sqrt( grad_rho(1,i_grid)**2 + grad_rho(2,i_grid)**2 + grad_rho(3,i_grid)**2 ) / &
(2.0D0 * kF(rho) * rho )
! -----------------------------------------------------------------
! This is the q value defined in equations 11 and 12 of DION.
! Use pw() from flib/functionals.f90 to get qc = kf/eps_x * eps_c.
!
call pw(length, r_s, 1, ec, dq0dr_v)
dq0_drho(i_grid) = dq0dr_v(1)
call pw(r_s, 1, ec, dq0_drho(i_grid))
!
q = -4.0D0*pi/3.0D0 * ec(1) + kF(rho) * Fs(s)
q = -4.0D0*pi/3.0D0 * ec + kF(rho) * Fs(s)
! -----------------------------------------------------------------
@ -674,7 +671,7 @@ CONTAINS
! dP_dq0 term will be interpolated later.
dq0_drho(i_grid) = dq0_dq * rho * ( -4.0D0*pi/3.0D0 * &
(dq0_drho(i_grid) - ec(1))/rho + dqx_drho(rho, s) )
(dq0_drho(i_grid) - ec)/rho + dqx_drho(rho, s) )
dq0_dgradrho(i_grid) = dq0_dq * rho * kF(rho) * dFs_ds(s) * ds_dgradrho(rho)
end do
@ -737,14 +734,13 @@ CONTAINS
real(dp), intent(OUT) :: dq0_dgradrho_up(:), dq0_dgradrho_down(:) ! Output variables.
complex(dp), intent(inout) :: thetas(:,:) ! The thetas from SOLER.
integer, parameter :: length=1 !^^^ PROVISIONAL (xc-lib)
real(dp) :: rho, up, down ! Local copy of densities.
real(dp) :: zeta(length) ! Spin polarization.
real(dp) :: r_s(length) ! Wigner-Seitz radius.
real(dp) :: zeta ! Spin polarization.
real(dp) :: r_s ! Wigner-Seitz radius.
real(dp) :: q, qc, qx, qx_up, qx_down ! q for exchange and correlation.
real(dp) :: q0x_up, q0x_down ! Saturated q values.
real(dp) :: fac
real(dp) :: ec(length), vc_v(length,2) !^^^ PROVISIONAL
real(dp) :: ec, vc(2)
real(dp) :: dq0_dq, dq0x_up_dq, dq0x_down_dq ! Derivative of q0 w.r.t q.
real(dp) :: dqc_drho_up, dqc_drho_down ! Intermediate values.
real(dp) :: dqx_drho_up, dqx_drho_down ! Intermediate values.
@ -821,14 +817,14 @@ CONTAINS
! This is the q value defined in equations 11 and 12 of DION and
! equation 8 of THONHAUSER (also see text above that equation).
r_s(1) = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
zeta(1) = (up - down) / rho
IF ( ABS(zeta(1)) > 1.0D0 ) zeta(1) = SIGN(1.0D0, zeta(1))
call pw_spin( length, r_s, zeta, ec, vc_v )
dqc_drho_up = vc_v(1,1) ; dqc_drho_down = vc_v(1,2)
r_s = ( 3.0D0 / (4.0D0*pi*rho) )**(1.0D0/3.0D0)
zeta = (up - down) / rho
IF ( ABS(zeta) > 1.0D0 ) zeta = SIGN(1.0D0, zeta)
call pw_spin( r_s, zeta, ec, vc )
dqc_drho_up = vc(1) ; dqc_drho_down = vc(2)
!
qx = ( up * q0x_up + down * q0x_down ) / rho
qc = -4.0D0*pi/3.0D0 * ec(1)
qc = -4.0D0*pi/3.0D0 * ec
q = qx + qc
@ -870,8 +866,8 @@ CONTAINS
if (calc_qx_down) dqx_drho_up = dqx_drho_up - q0x_down*down/rho
if (calc_qx_up) dqx_drho_down = dqx_drho_down - q0x_up *up /rho
dqc_drho_up = -4.0D0*pi/3.0D0 * (dqc_drho_up - ec(1))
dqc_drho_down = -4.0D0*pi/3.0D0 * (dqc_drho_down - ec(1))
dqc_drho_up = -4.0D0*pi/3.0D0 * (dqc_drho_up - ec)
dqc_drho_down = -4.0D0*pi/3.0D0 * (dqc_drho_down - ec)
dq0_drho_up (i_grid) = dq0_dq * (dqc_drho_up + dqx_drho_up )
dq0_drho_down(i_grid) = dq0_dq * (dqc_drho_down + dqx_drho_down)

6
PHonon/Gamma/cg_setup.f90
View File

@ -18,7 +18,7 @@ SUBROUTINE cg_setup
USE uspp_param, ONLY: upf
USE wavefunctions, ONLY: evc
USE io_files, ONLY: prefix, iunpun, iunres, diropn
USE funct, ONLY: dft_is_gradient
USE funct, ONLY: dft_is_gradient, init_lda_xc
USE dfunct, ONLY: newd
USE fft_base, ONLY: dfftp
USE gvect, ONLY: g, ngm, eigts1, eigts2, eigts3
@ -85,6 +85,8 @@ SUBROUTINE cg_setup
!
! derivative of the xc potential - NOT IMPLEMENTED FOR LSDA
!
CALL init_lda_xc()
!
rhotot(:) = rho%of_r(:,1) + rho_core(:)
!
sign_r = 1.0_DP
@ -98,7 +100,7 @@ SUBROUTINE cg_setup
ENDIF
ENDDO
!
CALL dmxc( dfftp%nnr, rhotot, dmuxc )
CALL dmxc_lda( dfftp%nnr, rhotot, dmuxc )
!
dmuxc = dmuxc * sign_r
!

1
PHonon/Gamma/make.depend
View File

@ -138,6 +138,7 @@ phcg.o : ../../Modules/constants.o
phcg.o : ../../Modules/control_flags.o
phcg.o : ../../Modules/environment.o
phcg.o : ../../Modules/fft_base.o
phcg.o : ../../Modules/funct.o
phcg.o : ../../Modules/io_files.o
phcg.o : ../../Modules/io_global.o
phcg.o : ../../Modules/ions_base.o

5
PHonon/Gamma/phcg.f90
View File

@ -421,6 +421,7 @@ SUBROUTINE cg_neweps
USE ions_base, ONLY : nat, tau
USE fft_base, ONLY : dfftp
USE scf, ONLY : rho, rho_core
USE funct, ONLY : init_lda_xc
USE cgcom
!
IMPLICIT NONE
@ -433,6 +434,8 @@ SUBROUTINE cg_neweps
!
! new derivative of the xc potential - NOT IMPLEMENTED FOR LSDA
!
CALL init_lda_xc()
!
rhotot(:) = rho%of_r(:,1) + rho_core(:)
!
sign_r = 1.0_DP
@ -446,7 +449,7 @@ SUBROUTINE cg_neweps
ENDIF
ENDDO
!
CALL dmxc( dfftp%nnr, rhotot, dmuxc )
CALL dmxc_lda( dfftp%nnr, rhotot, dmuxc )
!
dmuxc = dmuxc * sign_r
!

405
PP/src/benchmark_libxc.f90
View File

@ -8,11 +8,11 @@
!
PROGRAM benchmark_libxc
!
!**--------------------------------------------------------------------------------**
!** REMEMBER to comment the libxc blocks in 'Modules/correlation_lda_lsda.f90' and **
!** 'Modules/exchange_lda_lsda.f90' in order to run consistent tests (they should **
!** be commented already however). **
!**--------------------------------------------------------------------------------**
!------------------------------------------------------------------------------------!
! REMEMBER to comment eventual libxc blocks in the functional routines in 'Modules' !
! folder in order to run consistent tests (in the present version they should be !
! already absent, however). !
!------------------------------------------------------------------------------------!
!
#if defined(__LIBXC)
!
@ -23,94 +23,351 @@ PROGRAM benchmark_libxc
!
IMPLICIT NONE
!
INTEGER, PARAMETER :: DP = selected_real_kind(14,200)
!-------- Common vars ----------------------
INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(14,200)
INTEGER, PARAMETER :: nnr = 6
CHARACTER(LEN=120) :: e_q, f_q
INTEGER :: ii, ns, np, quit, i_sub
LOGICAL :: LDA, POLARIZED, ENERGY_ONLY, DF_OK
REAL(DP), PARAMETER :: null = 0.0_DP, pi34 = 0.6203504908994_DP
!
!----------QE vars --------------------------
INTEGER :: iexch_qe, icorr_qe
REAL(DP) :: rs(nnr)
REAL(DP), ALLOCATABLE :: rho_qe(:,:)
REAL(DP), ALLOCATABLE :: rho_tot(:), zeta(:)
REAL(DP), ALLOCATABLE :: ex_qe(:), ec_qe(:)
REAL(DP), ALLOCATABLE :: vx_qe(:,:), vc_qe(:,:)
REAL(DP), ALLOCATABLE :: dmuxc(:,:,:)
!
!--------- LIBXC vars -----------------------
TYPE(xc_f90_pointer_t) :: xc_func
TYPE(xc_f90_pointer_t) :: xc_info1, xc_info2
!
! The dimension of the arrays
INTEGER, PARAMETER :: nnr=10
!
REAL(DP), DIMENSION(nnr) :: ex_lda, vx_lda, ec_lda, vc_lda
!-----------------------------
REAL(DP), DIMENSION(nnr) :: rho, ex, vx, ec, vc
CHARACTER(LEN=120) :: name1, name2
INTEGER :: ii,nnr
INTEGER :: iexch_qe, icorr_qe, iexch_xc, icorr_xc
INTEGER :: iexch_lxc, icorr_lxc
INTEGER :: pol_unpol
REAL(DP), ALLOCATABLE :: rho_lxc(:)
REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
REAL(DP), ALLOCATABLE :: vx_lxc(:), vc_lxc(:)
REAL(DP), ALLOCATABLE :: dex_lxc(:), dcr_lxc(:), df_lxc(:)
!
!
!
! **************************************************************************
! *------------------------------------------------------------------------*
! * libxc funct. indexes: http://bigdft.org/Wiki/index.php?title=XC_codes *
! * qe " " : see comments in Modules/funct.f90 *
! *------------------------------------------------------------------------*
! * *
! * ... some examples for LDA case: *
! * *
! * *
! * | q-e | libxc | *
! * |___________|__________| *
! * slater (x) | 1 | 1 | *
! * pz (c) | 1 | 9 | *
! * | | | *
! * wigner (c) | 5 | 2 | *
! * vwn (c) | 2 | 7 | *
! * pw (c) | 4 | 12 | *
! * ... | ... | ... | *
! * *
! **************************************************************************
!
!
!
WRITE (*,'(/,1x,a)', ADVANCE='no') "Derivative of xc?(y/n) "
READ(*,*) f_q
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)
RETURN
ENDIF
!
IF ( .NOT.DF_OK ) THEN
WRITE (*,'(/,1x,a)', ADVANCE='no') "Energy only (y/n)? "
READ(*,*) e_q
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)
RETURN
ENDIF
ENDIF
!
WRITE (*,'(/,1x,a)', ADVANCE='no') "Polarization switch (1 unpolarized, &
& 2 polarized): "
READ(*,*) ns
IF ( ns/=1 .AND. ns/=2 ) THEN
PRINT *, CHAR(10)//"ERROR: you can only choose 1 or 2"//CHAR(10)
RETURN
ENDIF
WRITE (*,'(/,1x,a)') "-- Functional indexes "
WRITE (*,'(/,1x,a)', ADVANCE='no') "iexch_libxc icorr_libxc: "
READ(*,*) iexch_lxc, icorr_lxc
WRITE (*,'(/,1x,a)', ADVANCE='no') "iexch_qe icorr_qe: "
READ(*,*) iexch_qe, icorr_qe
IF (ns == 2 .AND. icorr_qe/=0 .AND. icorr_qe/=1 .AND. icorr_qe/=2 .AND. &
icorr_qe/=4 .AND. icorr_qe/=8 .AND. icorr_qe/=3 .AND. &
icorr_qe/=7 .AND. icorr_qe/=13) THEN
PRINT *, CHAR(10)//" ERROR: icorr_qe not available at these conditions"//CHAR(10)
RETURN
ENDIF
!
!
POLARIZED = .FALSE.
IF (ns == 2) THEN
POLARIZED = .TRUE.
ENDIF
!
pol_unpol = XC_UNPOLARIZED
np = 1
!
IF ( ns == 2 ) THEN
pol_unpol = XC_POLARIZED
np = 3
ENDIF
!
!
! *******************************************************************************
! * Polarized case: *
! * *
! * qe => rho_qe(:,1) = up | libxc => rho_lxc(2n+1) = up *
! * rho_qe(:,2) = down | (dim=2*nnr) rho_lxc(2n+2) = down *
! * | *
! * dmxc(:,1,1) = uu | (dim=3*nnr) df_lxc(3n+1) = uu *
! * dmxc(:,1,2) = ud | df_lxc(3n+2) = ud *
! * dmxc(:,2,1) = du | df_lxc(3n+2) = du *
! * dmxc(:,2,2) = dd | df_lxc(3n+3) = dd *
! * *
! *******************************************************************************
!
!
! ------ Allocations --------
!
! ... qe
!
ALLOCATE( rho_qe(nnr,ns) )
IF ( POLARIZED ) ALLOCATE( rho_tot(nnr), zeta(nnr) )
ALLOCATE( ex_qe(nnr), ec_qe(nnr) )
ALLOCATE( vx_qe(nnr,ns), vc_qe(nnr,ns) )
IF ( DF_OK ) THEN
IF ( .NOT.POLARIZED ) ALLOCATE( dmuxc(nnr,1,1) )
IF ( POLARIZED ) ALLOCATE( dmuxc(nnr,2,2) )
ENDIF
!
! ... libxc
!
ALLOCATE( rho_lxc(nnr*ns) )
ALLOCATE( ex_lxc(nnr), ec_lxc(nnr) )
ALLOCATE( vx_lxc(nnr*ns), vc_lxc(nnr*ns) )
IF ( DF_OK ) THEN
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
!
!
!----- Initializations -----
!
! ... qe
!
rho_qe = 0.0_DP
IF ( POLARIZED ) THEN
rho_tot = 0.0_DP
zeta = 0.0_DP
ENDIF
!
! ... libcx
!
rho_lxc = 0.0_DP
!
!
! -------- Setting up an arbitrary input for both qe and libxc -----
!
! ... qe
!
DO ii = 1, nnr
rho(ii)=DBLE(ii)/DBLE(nnr)
!
rho_qe(ii,1) = DBLE(ii)/DBLE(nnr+2)
!
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)
!
ENDIF
!
ENDDO
!
ex_lda = 0._dp ; ec_lda = 0._dp
vx_lda = 0._dp ; vc_lda = 0._dp
ex = 0._dp ; ec = 0._dp
vx = 0._dp ; vc = 0._dp
! ... libxc
!
!**--------------------------------------------------------------------------------**
!** libxc funct. indexes: http://bigdft.org/Wiki/index.php?title=XC_codes **
!** qe " " : see comments in Modules/funct.f90 **
!**--------------------------------------------------------------------------------**
!** **
!** ... some examples: **
!** **
!** | q-e | libxc | **
!** |___________|______________| **
!** exchange: | 1 | 1 | **
!** pz | 1 | 9 | **
!** wigner | 5 | 2 | **
!** vwn | 2 | 7 | **
!** pw | 4 | 12 | **
!** **
!**--------------------------------------------------------------------------------**
DO ii = 1, nnr
!
IF ( .NOT. POLARIZED ) THEN
!
rho_lxc(ii) = rho_qe(ii,1)
!
ELSE
!
rho_lxc(2*ii-1) = rho_qe(ii,1)
rho_lxc(2*ii) = rho_qe(ii,2)
!
ENDIF
!
ENDDO
!
!
PRINT *, "Insert functional indexes"
PRINT *, " "
PRINT *, "iexch_xc, icorr_xc: "
READ(*,*) iexch_xc, icorr_xc
PRINT *, "iexch_qe, icorr_qe: "
READ(*,*) iexch_qe, icorr_qe
!-------- Calculation of energies and potential -------------------
!
!------ LIBXC ------
!
CALL xc_f90_func_init( xc_func, xc_info1, iexch_xc, XC_UNPOLARIZED )
CALL xc_f90_lda_exc_vxc( xc_func, nnr, rho(1), ex_lda(1), vx_lda(1) )
CALL xc_f90_func_end( xc_func )
!
CALL xc_f90_func_init( xc_func, xc_info2, icorr_xc, XC_UNPOLARIZED )
CALL xc_f90_lda_exc_vxc( xc_func, nnr, rho(1), ec_lda(1), vc_lda(1) )
CALL xc_f90_func_end( xc_func )
!
!----- QE ----------
!
CALL select_functionals( iexch_qe, icorr_qe )
CALL xc( nnr, rho, ex, ec, vx, vc )
!CALL xc_spin( nnr, rho, zeta, ex, ec, vx, vc )
!
CALL xc_f90_func_init( xc_func, xc_info1, iexch_lxc, pol_unpol ) ! ... EXCHANGE
CALL xc_f90_lda_exc_vxc( xc_func, nnr, rho_lxc(1), ex_lxc(1), vx_lxc(1) )
!
IF ( DF_OK ) CALL xc_f90_lda_fxc( xc_func, nnr, rho_lxc(1), dex_lxc(1) )
CALL xc_f90_func_end( xc_func )
!
CALL xc_f90_func_init( xc_func, xc_info2, icorr_lxc, pol_unpol ) ! ... CORRELATION
CALL xc_f90_lda_exc_vxc( xc_func, nnr, rho_lxc(1), ec_lxc(1), vc_lxc(1) )
!
IF ( DF_OK ) THEN
CALL xc_f90_lda_fxc( xc_func, nnr, rho_lxc(1), dcr_lxc(1) )
df_lxc = (dex_lxc + dcr_lxc) * 2.0_DP
ENDIF
CALL xc_f90_func_end( xc_func )
!
!----- QE ----------
!
CALL select_lda_functionals( iexch_qe, icorr_qe ) ! ... EXCHANGE and CORRELATION
!
IF ( DF_OK ) THEN
IF ( .NOT.POLARIZED ) THEN
CALL dmxc_lda( nnr, rho_qe(:,1), dmuxc(:,1,1) )
ELSE
CALL dmxc_lsda( nnr, rho_qe, dmuxc )
ENDIF
ENDIF
!
IF ( .NOT. POLARIZED ) THEN
CALL xc_lda( nnr, rho_qe(:,1), ex_qe, ec_qe, vx_qe(:,1), vc_qe(:,1) )
ELSE
CALL xc_lsda( nnr, rho_tot, zeta, ex_qe, ec_qe, vx_qe, vc_qe )
ENDIF
!
!------------------
!
CALL xc_f90_info_name(xc_info1,name1)
CALL xc_f90_info_name(xc_info2,name2)
PRINT *, " "
PRINT *, "libxc: ", trim(name1), " + ", trim(name2)
PRINT *, " "
DO ii = 1, nnr, nnr
PRINT *, "rho:",rho(ii)
PRINT *, "--- ex, vx ---"
PRINT *, "libxc:",ex_lda(ii),vx_lda(ii)
PRINT *, "libqe:",ex(ii),vx(ii)
PRINT *, "diffs:",ex_lda(ii)-ex(ii),vx_lda(ii)-vx(ii)
PRINT *, "--- ec, vc ---"
PRINT *, "libxc:",ec_lda(ii),vc_lda(ii)
PRINT *, "libqe:",ec(ii),vc(ii)
PRINT *, "diffs:",ec_lda(ii)-ec(ii),vc_lda(ii)-vc(ii)
PRINT *, "--- ---"
ENDDO
CALL xc_f90_info_name( xc_info1, name1 )
CALL xc_f90_info_name( xc_info2, name2 )
!
PRINT *, "=================================== "//CHAR(10)//" "
PRINT *, "libxc functionals:"//CHAR(10)//" "
PRINT *, "Exchange: ", TRIM(name1)
PRINT *, "Correlation: ", TRIM(name2)
PRINT *, " "
!
DO ii = 1, nnr, nnr-1
WRITE(*,*) ' '
WRITE(*,*) ' '
WRITE(*,909) ii, nnr
IF ( .NOT. POLARIZED ) THEN
WRITE (*, 401 ) rho_qe(ii,1)
ELSE
WRITE (*, 402 ) rho_qe(ii,1), rho_qe(ii,2)
ENDIF
PRINT *, " "
IF (.NOT. DF_OK) THEN
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)
!
IF (.NOT. ENERGY_ONLY) THEN
PRINT *, " "
PRINT *, "=== Exchange potential ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) vx_qe(ii,1)
WRITE (*,201) vx_lxc(ii)
PRINT *, " --- "
WRITE (*,301) vx_qe(ii,1)-vx_lxc(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) vx_qe(ii,1), vx_qe(ii,2)
WRITE (*,202) vx_lxc(2*ii-1), vx_lxc(2*ii)
PRINT *, " --- "
WRITE (*,302) vx_qe(ii,1)-vx_lxc(2*ii-1), vx_qe(ii,2)-vx_lxc(2*ii)
ENDIF
PRINT *, " "
PRINT *, "=== Correlation potential ==="
IF ( .NOT. POLARIZED ) THEN
WRITE (*,101) vc_qe(ii,1)
WRITE (*,201) vc_lxc(ii)
PRINT *, " --- "
WRITE (*,301) vc_qe(ii,1)-vc_lxc(ii)
ELSEIF ( POLARIZED ) THEN
WRITE (*,102) vc_qe(ii,1), vc_qe(ii,2)
WRITE (*,202) vc_lxc(2*ii-1), vc_lxc(2*ii)
PRINT *, " --- "
WRITE (*,302) vc_qe(ii,1)-vc_lxc(2*ii-1), vc_qe(ii,2)-vc_lxc(2*ii)
ENDIF
PRINT *, "--- ---"
ENDIF
ELSE
PRINT *, " "
PRINT *, "=== First derivative of xc functional: ==="
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)
ENDIF
ENDIF
ENDDO
!
101 FORMAT('qe: ',3x,F17.14)
102 FORMAT('qe: ',3x,F17.14,4x,F17.14)
103 FORMAT('qe: ',3x,F17.14,4x,F17.14,4x,F17.14)
104 FORMAT('qe: ',3x,F17.14,4x,F17.14,4x,F17.14,4x,F17.14)
!
201 FORMAT('libxc: ',F17.14)
202 FORMAT('libxc: ',F17.14,4x,F17.14)
203 FORMAT('libxc: ',F17.14,4x,F17.14,4x,F17.14)
!
301 FORMAT('diffs: ',F17.14)
302 FORMAT('diffs: ',F17.14,4x,F17.14)
303 FORMAT('diffs: ',F17.14,4x,F17.14,4x,F17.14)
!
401 FORMAT('rho: ',F17.14)
402 FORMAT('rho(up,down): ',F17.14,4x,F17.14)
!
909 FORMAT('grid-point: ',I4,' of',I4)
!
DEALLOCATE( rho_qe )
IF ( POLARIZED ) DEALLOCATE( rho_tot, zeta )
DEALLOCATE( ex_qe, ec_qe )
DEALLOCATE( vx_qe, vc_qe )
IF (DF_OK) DEALLOCATE( dmuxc )
!
DEALLOCATE( rho_lxc )
DEALLOCATE( ex_lxc, ec_lxc )
DEALLOCATE( vx_lxc, vc_lxc )
IF (DF_OK) DEALLOCATE( dex_lxc, dcr_lxc, df_lxc )
!
PRINT *, " "
!
RETURN
!
#else
!
PRINT *, "ERROR: library libxc not found. Check Makefile."
PRINT *, "ERROR: library libxc not included."
RETURN
!
#endif
!

20
PP/src/ppacf.f90
View File

@ -51,7 +51,7 @@ PROGRAM do_ppacf
USE funct, ONLY : gcxc,gcx_spin,gcc_spin,dft_is_nonlocc,nlc
USE funct, ONLY : get_iexch, get_icorr, get_igcx, get_igcc
USE funct, ONLY : set_exx_fraction,set_auxiliary_flags,enforce_input_dft, init_lda_xc
USE xc_lda_lsda, ONLY : xc, xc_spin
USE xc_lda_lsda, ONLY : xc_lda, xc_lsda
USE wvfct, ONLY : npw, npwx
USE environment, ONLY : environment_start, environment_end
USE kernel_table, ONLY : Nqs, vdw_table_name, kernel_file_name
@ -337,7 +337,7 @@ PROGRAM do_ppacf
rs = pi34 /arhox**third
CALL slater(1, rs, ex, vx(:,1)) ! \epsilon_x,\lambda[n]=\epsilon_x[n]
ELSE
CALL xc( 1, arhox, ex, ec, vx(:,1), vc(:,1) )
CALL xc_lda( 1, arhox, ex, ec, vx(:,1), vc(:,1) )
ENDIF
etx=etx+e2*ex(1)*rhox
etxlda=etxlda+e2*ex(1)*rhox
@ -356,15 +356,15 @@ PROGRAM do_ppacf
IF(icorr==4) THEN
CALL pwcc(rs(1),cc,ec(1),vc(1,1),ec_l)
ELSE
CALL xc(1,arhox/ccp3,expp,ecpp,vx(:,1),vc(:,1))
CALL xc(1,arhox/ccm3,exm,ecm,vx(:,1),vc(:,1))
CALL xc_lda(1,arhox/ccp3,expp,ecpp,vx(:,1),vc(:,1))
CALL xc_lda(1,arhox/ccm3,exm,ecm,vx(:,1),vc(:,1))
!
ec_l=(ccp2*ecpp(1)-ccm2*ecm(1))/dcc*0.5_DP
ENDIF
etcldalambda=etcldalambda+e2*ec_l*rhox
IF(icc == ncc) THEN
IF(icorr.NE.4) THEN
CALL xc(1,arhox,ex,ec,vx(:,1),vc(:,1))
CALL xc_lda(1,arhox,ex,ec,vx(:,1),vc(:,1))
ENDIF
tclda%of_r(ir,1)=e2*(ec(1)-ec_l)*rhox
ttclda=ttclda+e2*(ec(1)-ec_l)*rhox
@ -377,7 +377,7 @@ PROGRAM do_ppacf
etcgclambda=etcgclambda+e2*ecgc_l*segno
ENDIF
ENDIF
CALL xc(1,arhox,ex,ec,vx(:,1),vc(:,1))
CALL xc_lda(1,arhox,ex,ec,vx(:,1),vc(:,1))
!
etclda=etclda+e2*ec(1)*rhox
etc=etc+e2*ec(1)*rhox
@ -418,7 +418,7 @@ PROGRAM do_ppacf
CALL slater_spin(1, arhox, zeta, ex, vx)
!
ELSE
CALL xc_spin( 1, arhox, zeta, ex, ec, vx, vc )
CALL xc_lsda( 1, arhox, zeta, ex, ec, vx, vc )
!
ENDIF
etx=etx+e2*ex(1)*rhox
@ -440,9 +440,9 @@ PROGRAM do_ppacf
IF(icorr==4) THEN
CALL pwcc_spin (rs(1),cc, zeta(1), ec(1), vc(1,1), vc(1,2), ec_l)
ELSE
CALL xc_spin( 1, arhox/ccp3, zeta, expp, ecpp, vx, vc )
CALL xc_lsda( 1, arhox/ccp3, zeta, expp, ecpp, vx, vc )
!
CALL xc_spin( 1, arhox/ccm3, zeta, exm, ecm, vx, vc )
CALL xc_lsda( 1, arhox/ccm3, zeta, exm, ecm, vx, vc )
!
ec_l=(ccp2*ecpp(1)-ccm2*ecm(1))/dcc*0.5_DP
ENDIF
@ -458,7 +458,7 @@ PROGRAM do_ppacf
etcgclambda=etcgclambda+e2*ecgc_l
ENDIF
ENDIF
CALL xc_spin( 1, arhox, zeta, ex, ec, vx, vc )
CALL xc_lsda( 1, arhox, zeta, ex, ec, vx, vc )
!
etclda=etclda+e2*ec(1)*rhox
etc=etc+e2*ec(1)*rhox

15
PP/src/xc_vdW_scale_mod.f90
View File

@ -319,8 +319,6 @@ CONTAINS
real(dp) :: dqc_drho
integer :: i_grid, idx ! Indexing variables.
!
real(dp) :: rs_v(1), ec_v(1), dqcdr_v(1) !^^^ PROVISIONAL (xc-lib)
@ -359,10 +357,7 @@ CONTAINS
! This is the q value defined in equations 11 and 12 of DION.
! Use pw() from flib/functionals.f90 to get qc = kf/eps_x * eps_c.
rs_v(1)=cc*r_s
call pw(1, rs_v, 1, ec_v, dqcdr_v)
ec=ec_v(1)
dqc_drho=dqcdr_v(1)
call pw( cc*r_s, 1, ec, dqc_drho)
!
q = -4.0D0*pi/3.0D0 * ec + kF(rho) * Fs(s)/cc
if(lecnl_qx) then
@ -452,7 +447,7 @@ CONTAINS
integer :: i_grid, idx ! Indexing variables
logical :: calc_qx_up, calc_qx_down
!
real(dp) :: rs_v(1), zeta_v(1), ec_v(1), vc_v(1,2)
real(dp) :: vc_v(2) ! auxiliary array for pw_spin call
@ -525,10 +520,8 @@ CONTAINS
zeta = (up - down) / rho
IF (ABS(zeta) > 1.0D0 ) zeta = SIGN(1.0D0, zeta)
!
rs_v(1)=cc*r_s ; zeta_v(1)=zeta
call pw_spin( 1, rs_v, zeta_v, ec_v, vc_v )
ec=ec_v(1)
dqc_drho_up=vc_v(1,1) ; dqc_drho_down=vc_v(1,2)
call pw_spin( cc*r_s, zeta, ec, vc_v )
dqc_drho_up=vc_v(1) ; dqc_drho_down=vc_v(2)
!
qx = ( up * q0x_up + down * q0x_down ) / rho / cc
qc = -4.0D0*pi/3.0D0 * ec

56
PW/src/paw_onecenter.f90
View File

@ -408,8 +408,8 @@ SUBROUTINE PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy)
USE uspp_param, ONLY : upf
USE lsda_mod, ONLY : nspin
USE atom, ONLY : g => rgrid
USE funct, ONLY : dft_is_gradient, init_lda_xc !^^^
USE xc_lda_lsda, ONLY : xc, xc_spin !^^^
USE funct, ONLY : dft_is_gradient, init_lda_xc
USE xc_lda_lsda, ONLY : xc
USE constants, ONLY : fpi ! REMOVE
TYPE(paw_info), INTENT(IN) :: i ! atom's minimal info
@ -435,7 +435,7 @@ SUBROUTINE PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy)
INTEGER :: mytid, ntids
!
!^^^****************************************** !^^^
REAL(DP), ALLOCATABLE :: arho(:), zeta(:), amag(:)
REAL(DP), ALLOCATABLE :: arho(:,:), zeta(:), amag(:)
REAL(DP), ALLOCATABLE :: ex(:), ec(:)
REAL(DP), ALLOCATABLE :: vx(:,:), vc(:,:)
REAL(DP), PARAMETER :: eps = 1.e-30_dp
@ -473,7 +473,7 @@ SUBROUTINE PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy)
!
ALLOCATE( rho_rad(i%m,nspin_mag) )
!
ALLOCATE( arho(i%m) ) !^^^
ALLOCATE( arho(i%m,2) ) !^^^
ALLOCATE( zeta(i%m) )
ALLOCATE( amag(i%m) )
ALLOCATE( ex(i%m) )
@ -505,18 +505,12 @@ SUBROUTINE PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy)
IF ( nspin_mag ==4 ) THEN
IF (with_small_so.AND.i%ae==1) CALL add_small_mag(i,ix,rho_rad)
!
!
DO k=1,i%m
rho_loc(k,1:nspin) = rho_rad(k,1:nspin)*g(i%t)%rm2(k)
amag(k) = SQRT(rho_loc(k,2)**2+rho_loc(k,3)**2+rho_loc(k,4)**2)
arho(k) = ABS( rho_loc(k,1)+rho_core(k) )
IF ( arho(k) > eps12 ) THEN
zeta(k) = amag(k) / arho(k)
IF ( ABS( zeta(k) ) > 1.D0 ) zeta(k) = SIGN( 1.D0, zeta(k) )
ENDIF
rho_loc(k,1) = rho_loc(k,1) + rho_core(k)
ENDDO
!
CALL xc_spin( i%m, arho, zeta, ex, ec, vx, vc )
CALL xc( i%m, 4, 2, rho_loc, ex, ec, vx, vc )
!
DO k=1,i%m
IF (present(energy)) &
@ -550,33 +544,25 @@ SUBROUTINE PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy)
!
IF (lsd==0) THEN
!
arho = ABS( rho_loc(:,1)+rho_core ) !^^^ def arhox, eps, e2,
arho(:,1) = rho_loc(:,1) + rho_core
!
CALL xc( i%m, arho, ex, ec, vx(:,1), vc(:,1) )
CALL xc( i%m, 1, 1, arho(:,1), ex, ec, vx(:,1), vc(:,1) )
!
v_rad(:,ix,1) = e2*( vx(:,1) + vc(:,1) )
if (present(energy)) e_rad = e2*( ex + ec )
if (present(energy)) e_rad = e2*( ex(:) + ec(:) )
!
ELSE
!
zeta=0.0_DP
DO k = 1, i%m
arho(k) = ABS( rho_loc(k,1)+rho_loc(k,2)+rho_core(k) )
IF ( arho(k) > eps ) THEN
zeta(k) = (rho_loc(k,1)-rho_loc(k,2)) / arho(k)
IF ( ABS(zeta(k)) > 1.D0 ) zeta(k) = SIGN( 1.D0, zeta(k) )
ENDIF
ENDDO
arho(:,1) = rho_loc(:,1) + rho_loc(:,2) + rho_core(:)
arho(:,2) = rho_loc(:,1) - rho_loc(:,2)
!
CALL xc_spin( i%m, arho, zeta, ex, ec, vx, vc )
CALL xc( i%m, 2, 2, arho, ex, ec, vx, vc )
!
v_rad(:,ix,:) = e2*( vx(:,:) + vc(:,:) ) !^^^**
v_rad(:,ix,:) = e2*( vx(:,:) + vc(:,:) )
if (present(energy)) e_rad(:) = e2*( ex(:) + ec(:) )
!
ENDIF
!
!CALL evxc_t_vec(rho_loc, rho_core, lsd, i%m, v_rad(:,ix,:), e_rad) !^^^
!
IF (present(energy)) THEN
IF ( nspin_mag < 2 ) THEN
e_rad = e_rad * ( rho_rad(:,1) + rho_core*g(i%t)%r2 )
@ -1605,18 +1591,19 @@ SUBROUTINE PAW_dpotential(dbecsum, becsum, int3, npe)
CALL stop_clock('PAW_dpot')
END SUBROUTINE PAW_dpotential
!
!
SUBROUTINE PAW_dxc_potential( i, drho_lm, rho_lm, rho_core, v_lm )
!
!! This routine computes the change of the exchange and correlation
!! potential in the spherical basis. It receives as input the charge
!! density and its variation.
!
USE spin_orb, ONLY : domag
USE noncollin_module, ONLY : nspin_mag
USE lsda_mod, ONLY : nspin
USE atom, ONLY : g => rgrid
USE funct, ONLY : dft_is_gradient
USE funct, ONLY : dft_is_gradient, init_lda_xc
!
TYPE(paw_info), INTENT(IN) :: i ! atom's minimal info
REAL(DP), INTENT(IN) :: rho_lm(i%m,i%l**2,nspin_mag) ! charge density as
@ -1644,6 +1631,7 @@ SUBROUTINE PAW_dxc_potential( i, drho_lm, rho_lm, rho_core, v_lm )
ALLOCATE(v_rad(i%m,rad(i%t)%nx,nspin_mag))
ALLOCATE(dmuxc(i%m,nspin_mag,nspin_mag))
!
CALL init_lda_xc()
!
DO ix = ix_s, ix_e
!
@ -1662,13 +1650,14 @@ SUBROUTINE PAW_dxc_potential( i, drho_lm, rho_lm, rho_core, v_lm )
CASE( 4 )
!
rho_rad(:,1) = rho_rad(:,1) + rho_core(:)
CALL dmxc_nc( i%m, rho_rad(:,1),rho_rad(:,2:4), dmuxc )
CALL dmxc( i%m, 4, rho_rad, dmuxc )
!
CASE( 2 )
!
rho_rad(:,1) = rho_rad(:,1) + 0.5_DP*rho_core(:)
rho_rad(:,2) = rho_rad(:,2) + 0.5_DP*rho_core(:)
CALL dmxc_spin( i%m, rho_rad, dmuxc )
!
CALL dmxc( i%m, 2, rho_rad, dmuxc )
!
CASE DEFAULT
!
@ -1685,7 +1674,7 @@ SUBROUTINE PAW_dxc_potential( i, drho_lm, rho_lm, rho_core, v_lm )
ENDIF
ENDDO
!
CALL dmxc( i%m, rho_rad(:,1), dmuxc(:,1,1) )
CALL dmxc( i%m, 1, rho_rad(:,1), dmuxc )
!
v_rad(:,ix,1) = dmuxc(:,1,1)*sign_r(:)
!
@ -1730,6 +1719,7 @@ SUBROUTINE PAW_dxc_potential( i, drho_lm, rho_lm, rho_core, v_lm )
!
END SUBROUTINE PAW_dxc_potential
!
!
! add gradient correction to dvxc. Both unpolarized and
! spin polarized cases are supported.
!

132
PW/src/v_of_rho.f90
View File

@ -347,7 +347,7 @@ SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
USE cell_base, ONLY : omega
USE spin_orb, ONLY : domag
USE funct, ONLY : nlc, dft_is_nonlocc, init_lda_xc
USE xc_lda_lsda, ONLY : xc, xc_spin
USE xc_lda_lsda, ONLY : xc
USE scf, ONLY : scf_type
USE mp_bands, ONLY : intra_bgrp_comm
USE mp, ONLY : mp_sum
@ -371,7 +371,9 @@ SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
!
REAL(DP) :: rhoneg(2), vs
!
REAL(DP), ALLOCATABLE :: arhox(:), amag(:), zeta(:)
!REAL(DP), ALLOCATABLE :: arhox(:), amag(:), zeta(:)
REAL(DP) :: arho, amag
REAL(DP) :: rhoup2, rhodw2
REAL(DP), ALLOCATABLE :: ex(:), ec(:)
REAL(DP), ALLOCATABLE :: vx(:,:), vc(:,:)
! In order:
@ -392,123 +394,87 @@ SUBROUTINE v_xc( rho, rho_core, rhog_core, etxc, vtxc, v )
!
CALL start_clock( 'v_xc' )
!
ALLOCATE(ex(dfftp%nnr))
ALLOCATE(ec(dfftp%nnr))
ALLOCATE(vx(dfftp%nnr,nspin))
ALLOCATE(vc(dfftp%nnr,nspin))
ALLOCATE( ex(dfftp%nnr) )
ALLOCATE( ec(dfftp%nnr) )
ALLOCATE( vx(dfftp%nnr,nspin) )
ALLOCATE( vc(dfftp%nnr,nspin) )
!
CALL init_lda_xc()
!
!
etxc = 0.D0
vtxc = 0.D0
v(:,:) = 0.D0
rhoneg = 0.D0
!
!
rho%of_r(:,1) = rho%of_r(:,1) + rho_core(:)
!
IF ( nspin == 1 .OR. ( nspin == 4 .AND. .NOT. domag ) ) THEN
!
! ... spin-unpolarized case
!
rho%of_r(:,1) = rho%of_r(:,1) + rho_core(:)
CALL xc( dfftp%nnr, 1, 1, rho%of_r(:,1), ex, ec, vx(:,1), vc(:,1) )
!
CALL xc( dfftp%nnr, ABS(rho%of_r), ex, ec, vx(:,1), vc(:,1) )
!
v(:,1) = e2*( vx(:,1) + vc(:,1) )
etxc = SUM( e2*((ex(:) + ec(:))*rho%of_r(:,1)) )
rho%of_r(:,1)=rho%of_r(:,1)-rho_core(:)
vtxc = SUM( v(:,1)*rho%of_r(:,1) )
DO ir=1, dfftp%nnr
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
DO ir = 1, dfftp%nnr
v(ir,1) = e2*( vx(ir,1) + vc(ir,1) )
etxc = etxc + e2*( ex(ir) + ec(ir) )*rho%of_r(ir,1)
rho%of_r(ir,1) = rho%of_r(ir,1) - rho_core(ir)
vtxc = vtxc + v(ir,1)*rho%of_r(ir,1)
IF (rho%of_r(ir,1) < 0.D0) rhoneg(1) = rhoneg(1)-rho%of_r(ir,1)
ENDDO
!
!
ELSE IF ( nspin == 2 ) THEN
!
ELSEIF ( nspin == 2 ) THEN
! ... spin-polarized case
!
ALLOCATE(arhox(dfftp%nnr))
ALLOCATE(zeta(dfftp%nnr))
zeta=0.D0
CALL xc( dfftp%nnr, 2, 2, rho%of_r, ex, ec, vx, vc )
!
DO ir = 1, dfftp%nnr
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
rho%of_r(ir,1) = rho%of_r(ir,1) + rho_core(ir)
arhox(ir) = ABS(rho%of_r(ir,1))
IF ( arhox(ir) > vanishing_charge ) THEN
zeta(ir) = rho%of_r(ir,2) / arhox(ir)
IF ( ABS(zeta(ir)) > 1.D0 ) THEN
rhoneg(2) = rhoneg(2) + 1.D0 / omega
zeta(ir) = SIGN( 1.D0, zeta(ir) )
END IF
ENDIF
ENDDO
!
CALL xc_spin( dfftp%nnr, arhox, zeta, ex, ec, vx, vc )
!
DO ir = 1, dfftp%nnr
DO ir = 1, dfftp%nnr !OMP ?
v(ir,:) = e2*( vx(ir,:) + vc(ir,:) )
etxc = etxc + e2*( (ex(ir) + ec(ir))*rho%of_r(ir,1) )
rho%of_r(ir,1) = rho%of_r(ir,1) - rho_core(ir)
vtxc = vtxc + ( ( v(ir,1) + v(ir,2) )*rho%of_r(ir,1) + &
( v(ir,1) - v(ir,2) )*rho%of_r(ir,2) )
!
rhoup2 = rho%of_r(ir,1)+rho%of_r(ir,2)
rhodw2 = rho%of_r(ir,1)-rho%of_r(ir,2)
IF (rhoup2 < 0.d0) rhoneg(1) = rhoneg(1) + rhoup2
IF (rhodw2 < 0.d0) rhoneg(2) = rhoneg(2) + rhodw2
ENDDO
vtxc = 0.5d0 * vtxc
!
DEALLOCATE(arhox)
DEALLOCATE(zeta)
!
vtxc = 0.5d0 * vtxc
rhoneg = 0.5d0 * rhoneg
!
!
ELSE IF ( nspin == 4 ) THEN
!
! ... noncolinear case
!
ALLOCATE(arhox(dfftp%nnr))
ALLOCATE(amag(dfftp%nnr))
ALLOCATE(zeta(dfftp%nnr))
zeta=0.D0
CALL xc( dfftp%nnr, 4, 2, rho%of_r, ex, ec, vx, vc )
!
DO ir = 1, dfftp%nnr
amag(ir) = SQRT( rho%of_r(ir,2)**2 + rho%of_r(ir,3)**2 + rho%of_r(ir,4)**2 )
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
rho%of_r(ir,1) = rho%of_r(ir,1) + rho_core(ir)
arhox(ir) = ABS( rho%of_r(ir,1) )
IF ( arhox(ir) > vanishing_charge ) THEN
zeta(ir) = amag(ir) / arhox(ir)
IF ( ABS( zeta(ir) ) > 1.D0 ) THEN
rhoneg(2) = rhoneg(2) + 1.D0 / omega
zeta(ir) = SIGN( 1.D0, zeta(ir) )
END IF
DO ir = 1, dfftp%nnr !OMP ?
arho = ABS( rho%of_r(ir,1) )
IF ( arho < vanishing_charge ) CYCLE
vs = 0.5D0*( vx(ir,1) + vc(ir,1) - vx(ir,2) - vc(ir,2) )
v(ir,1) = e2*( 0.5D0*( vx(ir,1) + vc(ir,1) + vx(ir,2) + vc(ir,2) ) )
!
amag = SQRT( SUM( rho%of_r(ir,2:4)**2 ) )
IF ( amag > vanishing_mag ) THEN
v(ir,2:4) = e2 * vs * rho%of_r(ir,2:4) / amag
vtxc = vtxc + SUM( v(ir,2:4) * rho%of_r(ir,2:4) )
ENDIF
etxc = etxc + e2*( ex(ir) + ec(ir) ) * arho
!
rho%of_r(ir,1) = rho%of_r(ir,1) - rho_core(ir)
IF ( rho%of_r(ir,1) < 0.D0 ) rhoneg(1) = rhoneg(1) - rho%of_r(ir,1)
IF ( amag / arho > 1.D0 ) rhoneg(2) = rhoneg(2) + 1.D0/omega
vtxc = vtxc + v(ir,1) * rho%of_r(ir,1)
ENDDO
!
CALL xc_spin( dfftp%nnr, arhox, zeta, ex, ec, vx, vc )
!
DO ir = 1, dfftp%nnr
IF ( arhox(ir) > vanishing_charge ) THEN
vs = 0.5D0*( vx(ir,1) + vc(ir,1) - vx(ir,2) - vc(ir,2) )
v(ir,1) = e2*( 0.5D0*( vx(ir,1) + vc(ir,1) + vx(ir,2) + vc(ir,2) ) )
IF ( amag(ir) > vanishing_mag ) THEN
DO ipol = 2, 4
v(ir,ipol) = e2 * vs * rho%of_r(ir,ipol) / amag(ir)
vtxc = vtxc + v(ir,ipol) * rho%of_r(ir,ipol)
ENDDO
END IF
etxc = etxc + e2*( ex(ir) + ec(ir) ) * rho%of_r(ir,1)
rho%of_r(ir,1) = rho%of_r(ir,1) - rho_core(ir)
vtxc = vtxc + v(ir,1) * rho%of_r(ir,1)
ENDIF
ENDDO
!
DEALLOCATE(arhox)
DEALLOCATE(amag)
DEALLOCATE(zeta)
!
END IF
ENDIF
!
DEALLOCATE(ex)
DEALLOCATE(ec)
DEALLOCATE(vx)
DEALLOCATE(vc)
DEALLOCATE( ex, ec )
DEALLOCATE( vx, vc )
!
CALL mp_sum( rhoneg , intra_bgrp_comm )
!

7
atomic/src/c6_tfvw.f90
View File

@ -420,7 +420,7 @@ subroutine dvxc_dn(mesh, rho, dvxc)
! compute the derivative of xc-potential w.r.t local density.
! some routine in PH and flibs will be called
!
use funct, only : dft_is_gradient
use funct, only : dft_is_gradient, init_lda_xc
!
implicit none
!
@ -439,7 +439,10 @@ subroutine dvxc_dn(mesh, rho, dvxc)
call errore ('dvxc_dn', 'gradient correction to dvxc not yet implemented', 1)
!
! LDA only
CALL dmxc( mesh, rho, dvxc )
!
CALL init_lda_xc()
!
CALL dmxc_lda( mesh, rho, dvxc )
!
return
!

43
atomic/src/vxcgc.f90
View File

@ -15,15 +15,15 @@ subroutine vxc_t(lsd,rho,rhoc,exc,vxc)
!
use kinds, only : DP
use funct, only : init_lda_xc
use xc_lda_lsda, only: xc, xc_spin
use xc_lda_lsda, only: xc
implicit none
integer, intent(in) :: lsd ! 1 in the LSDA case, 0 otherwise
real(DP), intent(in) :: rho(2), rhoc ! the system density
real(DP), intent(out):: exc, vxc(2)
real(DP), intent(out):: exc(1), vxc(2) !^^^
!
integer, parameter :: length=1 !^^^ PROVISIONAL (xc-lib)
real(DP), dimension(length) :: arho, zeta, ex, ec
REAL(DP), dimension(length,2) :: vx, vc
integer, parameter :: length=1 !^^^ PROVISIONAL (xc-lib)
real(DP), dimension(length) :: ex, ec , arho
REAL(DP), dimension(length,2) :: rhoaux, vx, vc
!
real(DP), parameter :: e2=2.0_dp, eps=1.e-30_dp
!
@ -36,10 +36,10 @@ subroutine vxc_t(lsd,rho,rhoc,exc,vxc)
!
! LDA case
!
arho(1) = abs(rho(1) + rhoc)
if (arho(1) > eps) then
rhoaux(1,1) = abs(rho(1) + rhoc)
if (rhoaux(1,1) > eps) then
!
call xc( length, arho ,ex, ec, vx(:,1), vc(:,1) )
CALL xc( length, 1, 1, rhoaux, ex, ec, vx(:,1), vc(:,1) )
!
vxc(1) = e2 * ( vx(1,1) + vc(1,1) )
exc = e2 * ( ex(1) + ec(1) )
@ -49,24 +49,23 @@ subroutine vxc_t(lsd,rho,rhoc,exc,vxc)
! LSDA case
!
vxc(2)=0.0_dp
arho(1) = abs(rho(1)+rho(2)+rhoc)
if (arho(1) > eps) then
zeta(1) = (rho(1)-rho(2)) / arho(1)
! zeta has to stay between -1 and 1, but can get a little
! out the bound during the first iterations.
if (abs(zeta(1)) > 1.0_dp) zeta(1) = sign(1._dp, zeta(1))
!
CALL xc_spin( length, arho, zeta, ex, ec, vx, vc )
!
vxc(1) = e2 * ( vx(1,1) + vc(1,1) )
vxc(2) = e2 * ( vx(1,2) + vc(1,2) )
exc = e2 * ( ex(1) + ec(1) )
endif
!
rhoaux(1,1) = rho(1) + rho(2) + rhoc
rhoaux(1,2) = rho(1) - rho(2)
!
CALL xc( length, 2, 2, rhoaux, ex, ec, vx, vc )
!
vxc(1) = e2 * ( vx(1,1) + vc(1,1) )
vxc(2) = e2 * ( vx(1,2) + vc(1,2) )
exc = e2 * ( ex(1) + ec(1) )
!
endif
!
return
!
end subroutine vxc_t
!
!
!---------------------------------------------------------------
subroutine vxcgc ( ndm, mesh, nspin, r, r2, rho, rhoc, vgc, egc, &
tau, vtau, iflag)