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

XClib - q-e/libxc threshold alignment (attempt)

This commit is contained in:
fabrizio22 2021-10-01 17:02:04 +02:00
parent 843d85f7e6
commit 3dc87c62eb
5 changed files with 91 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
XClib/qe_drivers_mgga.f90
View File

@ -79,6 +79,9 @@ SUBROUTINE tau_xc( length, rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c
IF ( (arho<=rho_threshold_mgga).OR.(grho2(k)<=grho2_threshold_mgga).OR. &
(ABS(tau(k))<=rho_threshold_mgga) ) CYCLE
!
! ...libxc-like threshold management
!grho2(k) = MIN( grho2(k), (8.d0*rho(k)*tau(k))**2 )
!
SELECT CASE( imeta )
CASE( 1 )
!
@ -147,14 +150,15 @@ SUBROUTINE tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, v2c,
ex=0.0_DP ; v1x=0.0_DP ; v2x=0.0_DP ; v3x=0.0_DP
ec=0.0_DP ; v1c=0.0_DP ; v2c=0.0_DP ; v3c=0.0_DP
!
! FIXME: for SCAN, this will be calculated later
!
DO k = 1, length
!
rh = rho(k,1) + rho(k,2)
atau = tau(k,1) + tau(k,2) ! KE-density in Hartree
grho2(1) = SUM( grho(:,k,1)**2 )
grho2(2) = SUM( grho(:,k,2)**2 )
! ...libxc-like threshold management
!grho2(1) = MIN( SUM(grho(:,k,1)**2), (8.d0*rho(k,1)*tau(k,1))**2 )
!grho2(2) = MIN( SUM(grho(:,k,2)**2), (8.d0*rho(k,2)*tau(k,2))**2 )
grho2(1) = SUM(grho(:,k,1)**2)
grho2(2) = SUM(grho(:,k,2)**2)
ggrho2 = ( grho2(1) + grho2(2) ) * 4.0_DP
!
IF ( (rh <= rho_threshold_mgga).OR.(ggrho2 <= grho2_threshold_mgga).OR.&

14
XClib/qe_funct_mgga.f90
View File

@ -267,9 +267,7 @@ SUBROUTINE metac( rho, grho2, tau, ec, v1c, v2c, v3c ) !<GPU:
cf3 = cf3*rho
v2c = cf2*d2rev + cf3*2.0_DP/grho2
v3c = cf2*d3rev - cf3/tau
!
ec = rho*ec_rev*(1.0_DP+dd*ec_rev*z2*z) !-rho*ec_unif(1)
!v1c = v1c - vc_unif(1)
!
RETURN
!
@ -865,8 +863,8 @@ END SUBROUTINE metac_spin
!
! input: - rho
! - grho2=|\nabla rho|^2
! - tau = the kinetic energy density
! It is defined as summ_i( |nabla phi_i|**2 )
! - tau = the input kinetic energy density
! It is defined as 0.5summ_i( |nabla phi_i|**2 )
!
! definition: E_x = \int ex dr
!
@ -923,18 +921,20 @@ SUBROUTINE m06lxc( rho, grho2, tau, ex, ec, v1x, v2x, v3x, v1c, v2c, v3c )
!
taua = tau * two * 0.5_dp ! Taua, which is Tau_sigma is half Tau
taub = taua ! Tau is defined as summ_i( |nabla phi_i|**2 )
! in the M06L routine
! in the following M06L routines
!
CALL m06lx( rhoa, grho2a, taua, ex, v1x, v2x, v3x ) !<GPU:m06lx=>m06lx_d>
!
ex = two * ex ! Add the two components up + dw
!
v2x = 0.5_dp * v2x
v3x = 2.0_dp * v3x
!
CALL m06lc( rhoa, rhob, grho2a, grho2b, taua, taub, ec, v1c, v2c, v3c, & !<GPU:m06lc=>m06lc_d>
v1cb, v2cb, v3cb )
!
v2c = 0.5_dp * v2c
v3c = 2.0_dp * v3c
!
END SUBROUTINE m06lxc
!
@ -963,9 +963,13 @@ SUBROUTINE m06lxc_spin( rhoup, rhodw, grhoup2, grhodw2, tauup, taudw, &
CALL m06lx( rhodw, grhodw2, taub, exdw, v1xdw, v2xdw, v3xdw ) !<GPU:m06lx=>m06lx_d>
!
ex = exup + exdw
v3xup = 2.0_dp * v3xup
v3xdw = 2.0_dp * v3xdw
!
CALL m06lc( rhoup, rhodw, grhoup2, grhodw2, taua, taub, & !<GPU:m06lc=>m06lc_d>
ec, v1cup, v2cup, v3cup, v1cdw, v2cdw, v3cdw )
v3cup = 2.0_dp * v3cup
v3cdw = 2.0_dp * v3cdw
!
END SUBROUTINE m06lxc_spin
! !

30
XClib/xc_wrapper_gga.f90
View File

@ -174,19 +174,23 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
ENDDO
ELSE
DO k = 1, length
sgn(:) = 1.d0
!IF (rho_lxc(2*k-1)<rho_threshold_gga .OR. SQRT(ABS(sigma(3*k-2)))<grho_threshold_gga) sgn(1)=0.d0
!IF (rho_lxc(2*k) <rho_threshold_gga .OR. SQRT(ABS(sigma(3*k))) <grho_threshold_gga) sgn(2)=0.d0
IF (rho_lxc(2*k-1)<rho_threshold_gga) sgn(1)=0.d0
!.OR. SQRT(ABS(sigma(3*k-2)))<grho_threshold_gga) sgn(1)=0.d0
IF (rho_lxc(2*k) <rho_threshold_gga) sgn(2)=0.d0
!.OR. SQRT(ABS(sigma(3*k))) <grho_threshold_gga) sgn(2)=0.d0
ec(k) = ec_lxc(k) * (rho_lxc(2*k-1)*sgn(1)+rho_lxc(2*k)*sgn(2))
v1c(k,1) = vc_rho(2*k-1) * sgn(1)
v1c(k,2) = vc_rho(2*k) * sgn(2)
v2c(k,1) = vc_sigma(3*k-2)*2.d0 * sgn(1)
v2c_ud(k)= vc_sigma(3*k-1) * sgn(1)*sgn(2)
v2c(k,2) = vc_sigma(3*k)*2.d0 * sgn(2)
rho_up = rho_lxc(2*k-1)
rho_dw = rho_lxc(2*k)
grho2_up = sigma(3*k-2)
grho2_dw = sigma(3*k)
xlda_up = 1.0_DP ; xlda_dw = 1.0_DP
xgga_up = 1.0_DP ; xgga_dw = 1.0_DP
IF ( rho_up+rho_dw <= rho_threshold_gga ) CYCLE
IF ( rho_up <= rho_threshold_lda ) xlda_up = 0.0_DP
IF ( rho_dw <= rho_threshold_lda ) xlda_dw = 0.0_DP
IF ( rho_up<=small .OR. SQRT(ABS(grho2_up))<=small ) xgga_up = 0.0_DP
IF ( rho_dw<=small .OR. SQRT(ABS(grho2_dw))<=small ) xgga_dw = 0.0_DP
ec(k) = ec_lxc(k) * (rho_up*xlda_up+rho_dw*xlda_dw)
v1c(k,1) = vc_rho(2*k-1)*xlda_up
v1c(k,2) = vc_rho(2*k)*xlda_dw
v2c(k,1) = vc_sigma(3*k-2)*2.d0*xgga_up
v2c_ud(k)= vc_sigma(3*k-1)*xgga_up*xgga_dw
v2c(k,2) = vc_sigma(3*k)*2.d0*xgga_dw
ENDDO
ENDIF
!

82
XClib/xc_wrapper_mgga.f90
View File

@ -19,9 +19,9 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
#endif
!
USE kind_l, ONLY: DP
USE dft_setting_params, ONLY: imeta, imetac, is_libxc, rho_threshold_mgga, &
grho2_threshold_mgga, tau_threshold_mgga, scan_exx, &
exx_started, exx_fraction
USE dft_setting_params, ONLY: imetac, is_libxc, rho_threshold_mgga, &
grho2_threshold_mgga, tau_threshold_mgga, &
scan_exx, exx_started, exx_fraction
USE qe_drivers_mgga
!
IMPLICIT NONE
@ -59,16 +59,18 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
!
INTEGER :: is
REAL(DP), ALLOCATABLE :: grho2(:,:)
REAL(DP), PARAMETER :: small = 1.E-10_DP
!
#if defined(__LIBXC)
REAL(DP), ALLOCATABLE :: rho_lxc(:), sigma(:), tau_lxc(:)
REAL(DP), ALLOCATABLE :: ex_lxc(:), ec_lxc(:)
REAL(DP), ALLOCATABLE :: vx_rho(:), vx_sigma(:), vx_tau(:)
REAL(DP), ALLOCATABLE :: vc_rho(:), vc_sigma(:), vc_tau(:)
REAL(DP), ALLOCATABLE :: lapl_rho(:), vlapl_rho(:) ! not used in TPSS
REAL(DP), ALLOCATABLE :: lapl_rho(:), vlapl_rho(:) ! not used in QE
!
INTEGER :: k, ipol, pol_unpol, eflag
LOGICAL :: POLARIZED
INTEGER :: k, ipol, pol_unpol, eflag
LOGICAL :: POLARIZED
REAL(DP) :: rh, ggrho2, atau
#if (XC_MAJOR_VERSION > 4)
INTEGER(8) :: lengthxc
#else
@ -103,6 +105,7 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
sigma(k) = MAX( grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2, &
grho2_threshold_mgga )
tau_lxc(k) = MAX( tau(k,1), tau_threshold_mgga )
vlapl_rho(k) = 0._DP
ENDDO
!
ELSE
@ -110,16 +113,16 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
DO k = 1, length
rho_lxc(2*k-1) = ABS( rho(k,1) )
rho_lxc(2*k) = ABS( rho(k,2) )
!
sigma(3*k-2) = MAX( grho(1,k,1)**2 + grho(2,k,1)**2 + grho(3,k,1)**2, &
grho2_threshold_mgga )
sigma(3*k-1) = grho(1,k,1) * grho(1,k,2) + grho(2,k,1) * grho(2,k,2) + &
sigma(3*k-1) = grho(1,k,1) * grho(1,k,2) + grho(2,k,1) * grho(2,k,2) +&
grho(3,k,1) * grho(3,k,2)
sigma(3*k) = MAX( grho(1,k,2)**2 + grho(2,k,2)**2 + grho(3,k,2)**2, &
grho2_threshold_mgga )
!
tau_lxc(2*k-1) = MAX( tau(k,1), tau_threshold_mgga )
tau_lxc(2*k) = MAX( tau(k,2), tau_threshold_mgga )
tau_lxc(2*k-1) = MAX( tau(k,1), small )
tau_lxc(2*k) = MAX( tau(k,2), small )
vlapl_rho(2*k-1) = 0._DP
vlapl_rho(2*k) = 0._DP
ENDDO
!
ENDIF
@ -132,19 +135,17 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
grho2(:,is) = grho(1,:,is)**2 + grho(2,:,is)**2 + grho(3,:,is)**2
ENDDO
!
IF (ns == 1) THEN
CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, v1x(:,1), &
v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1), v3c(:,1) )
ELSEIF (ns == 2) THEN
CALL tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1c, &
v2c, v3c )
ENDIF
IF (ns == 1) CALL tau_xc( length, rho(:,1), grho2(:,1), tau(:,1), ex, ec, &
v1x(:,1), v2x(:,1), v3x(:,1), v1c(:,1), v2c(1,:,1),&
v3c(:,1) )
IF (ns == 2) CALL tau_xc_spin( length, rho, grho, tau, ex, ec, v1x, v2x, v3x,&
v1c, v2c, v3c )
!
DEALLOCATE( grho2 )
!
ENDIF
!
! META EXCHANGE
! ... META EXCHANGE
!
IF ( is_libxc(5) ) THEN
CALL xc_f03_func_set_dens_threshold( xc_func(5), rho_threshold_mgga )
@ -160,6 +161,9 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
IF ( ABS(rho_lxc(k))<=rho_threshold_mgga .OR. &
sigma(k)<=grho2_threshold_mgga .OR. &
ABS(tau_lxc(k))<=rho_threshold_mgga ) CYCLE
ex(k) = ex_lxc(k) * rho_lxc(k)
v1x(k,1) = vx_rho(k)
v2x(k,1) = vx_sigma(k) * 2.0_DP
@ -167,13 +171,22 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
ENDDO
ELSE
DO k = 1, length
IF (rho_lxc(2*k-1)+rho_lxc(2*k) <= rho_threshold_mgga) CYCLE
ex(k) = ex_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
v1x(k,1) = vx_rho(2*k-1)
v1x(k,2) = vx_rho(2*k)
v2x(k,1) = vx_sigma(3*k-2)*2.d0
v2x(k,2) = vx_sigma(3*k)*2.d0
v3x(k,1) = vx_tau(2*k-1)
v3x(k,2) = vx_tau(2*k)
IF ( ABS(rho_lxc(2*k-1))>rho_threshold_mgga .AND. &
sigma(3*k-2)>grho2_threshold_mgga .AND. &
ABS(tau_lxc(2*k-1))>tau_threshold_mgga ) THEN
v1x(k,1) = vx_rho(2*k-1)
v2x(k,1) = vx_sigma(3*k-2)*2.d0
v3x(k,1) = vx_tau(2*k-1)
ENDIF
IF ( ABS(rho_lxc(2*k))>rho_threshold_mgga .AND. &
sigma(3*k)>grho2_threshold_mgga .AND. &
ABS(tau_lxc(2*k))>tau_threshold_mgga ) THEN
v1x(k,2) = vx_rho(2*k)
v2x(k,2) = vx_sigma(3*k)*2.d0
v3x(k,2) = vx_tau(2*k) !m06l: *2 libxc o /2 qe????
ENDIF !-chiarisci il problema del fattore 2
ENDDO
ENDIF
!
@ -189,7 +202,7 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
!
ENDIF
!
! META CORRELATION
! ... META CORRELATION
!
IF ( is_libxc(6) ) THEN
!
@ -206,13 +219,22 @@ SUBROUTINE xc_metagcx( length, ns, np, rho, grho, tau, ex, ec, v1x, v2x, v3x, v1
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
ec(k) = ec_lxc(k) * rho_lxc(k)
v1c(k,1) = vc_rho(k)
v2c(1,k,1) = vc_sigma(k) * 2.0_DP
v3c(k,1) = vc_tau(k)
IF ( ABS(rho_lxc(k))<=rho_threshold_mgga .OR. &
sigma(k)<=grho2_threshold_mgga .OR. &
ABS(tau_lxc(k))<=rho_threshold_mgga ) CYCLE
ec(k) = ec_lxc(k) * rho_lxc(k)
v1c(k,1) = vc_rho(k)
v2c(1,k,1) = vc_sigma(k) * 2.0_DP
v3c(k,1) = vc_tau(k)
ENDDO
ELSE
DO k = 1, length
rh = rho_lxc(2*k-1) + rho_lxc(2*k)
atau = ABS(tau_lxc(2*k-1) + tau_lxc(2*k))
ggrho2 = (sigma(3*k-2) + sigma(3*k))*4.0_DP
IF ( rh <= rho_threshold_mgga .OR. &
ggrho2 <= grho2_threshold_mgga .OR. &
atau <= tau_threshold_mgga ) CYCLE
ec(k) = ec_lxc(k) * (rho_lxc(2*k-1)+rho_lxc(2*k))
v1c(k,1) = vc_rho(2*k-1)
v1c(k,2) = vc_rho(2*k)

15
XClib/xclib_test.f90
View File

@ -433,7 +433,7 @@ PROGRAM xclib_test
!
IF ( TRIM(dft)=='xxxx' .OR.TRIM(dft)=='NONE' .OR. &
TRIM(dft)=='META' .OR. & !TRIM(dft)=='TB09' .OR. &
TRIM(dft)=='SCA0' .OR.TRIM(dft)=='TPSS' .OR. &
TRIM(dft)=='SCA0' .OR. & !TRIM(dft)=='TPSS' .OR. &
TRIM(dft)=='SCAN0'.OR.TRIM(dft)=='PZ+META'.OR. &
TRIM(dft)=='PBE+META' ) THEN
CALL print_test_status( skipped )
@ -1469,10 +1469,10 @@ PROGRAM xclib_test
IF (what(1:1)=='V') thr = diff_thr_vmgga
ENDIF
!
!IF (mype==root .AND. test/='gen-benchmark') THEN
! WRITE(stdout,*) " "
! IF ( POLARIZED .AND. what(1:1)/='E' ) WRITE(stdout,*) TRIM(what)
!ENDIF
IF (mype==root .AND. TRIM(test)=='exe-benchmark') THEN
WRITE(stdout,*) " "
IF ( POLARIZED .AND. what(1:1)/='E' ) WRITE(stdout,*) TRIM(what)
ENDIF
!
xc_aver=0._DP
!
@ -1545,11 +1545,6 @@ PROGRAM xclib_test
IF (LDA .AND. .NOT.GGA) thr = diff_thr_dmuxc
IF ( GGA ) thr = diff_thr_dv
!
!IF (test=='exe-benchmark') THEN
! WRITE(stdout,*) " "
! WRITE(stdout,*) what
!ENDIF
!
dxc_aver=0._DP
!
dxc_aver(1,1) = SUM(dxc_qe(1:nnr,1,1))/DBLE(npoints)