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

XClib - libxc gga thresholds again

This commit is contained in:
fabrizio22 2021-10-13 12:02:47 +02:00
parent e5d3914db6
commit 176dc91d93
4 changed files with 92 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
XClib/xc_wrapper_d_gga.f90
View File

@ -58,9 +58,10 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
!
INTEGER :: k, length_lxc, length_dlxc
REAL(DP) :: rht, zeta
LOGICAL :: thr_dw_cond, thr_up_cond
REAL(DP), ALLOCATABLE :: sigma(:)
REAL(DP), PARAMETER :: small = 1.E-10_DP, rho_trash = 0.5_DP
REAL(DP), PARAMETER :: epsr=1.0d-6, epsg=1.0d-10
REAL(DP), PARAMETER :: epsr=1.0d-6, epsg=1.0d-6
!
IF ( ANY(.NOT.is_libxc(3:4)) ) THEN
rho_threshold_gga = small ; grho_threshold_gga = small
@ -117,7 +118,8 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
v2sigma2_x(length_dlxc*sp) )
! ... DERIVATIVE FOR EXCHANGE
v2rho2_x = 0._DP ; v2rhosigma_x = 0._DP ; v2sigma2_x = 0._DP
IF (igcx /= 0) THEN
IF (igcx /= 0) THEN
CALL xc_f03_func_set_dens_threshold( xc_func(3), epsr )
CALL xc_f03_gga_fxc( xc_func(3), lengthxc, rho_lbxc(1), sigma(1), v2rho2_x(1), &
v2rhosigma_x(1), v2sigma2_x(1) )
ENDIF
@ -130,6 +132,7 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
v2rho2_c = 0._DP ; v2rhosigma_c = 0._DP ; v2sigma2_c = 0._DP
IF (igcc /= 0) THEN
fkind = xc_f03_func_info_get_kind( xc_info(4) )
CALL xc_f03_func_set_dens_threshold( xc_func(4), epsr )
CALL xc_f03_gga_fxc( xc_func(4), lengthxc, rho_lbxc(1), sigma(1), v2rho2_c(1), &
v2rhosigma_c(1), v2sigma2_c(1) )
ENDIF
@ -199,7 +202,6 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
!
ELSEIF ( sp == 2 ) THEN
!
!IF ( ((.NOT.is_libxc(3)) .OR. (.NOT.is_libxc(4))) ) THEN
IF ( ((.NOT.is_libxc(3).AND.igcx/=0) .OR. (.NOT.is_libxc(4).AND.igcc/=0)) &
.AND. fkind/=XC_EXCHANGE_CORRELATION ) THEN
!
@ -234,35 +236,34 @@ SUBROUTINE dgcxc( length, sp, r_in, g_in, dvxc_rr, dvxc_sr, dvxc_ss )
!
IF ( is_libxc(3) ) THEN
!
DO k = 1, length
rht = r_in(k,1) + r_in(k,2)
IF (rht > epsr) THEN
dvxc_rr(k,1,1) = dvxc_rr(k,1,1) + e2 * v2rho2_x(3*k-2)
dvxc_rr(k,1,2) = dvxc_rr(k,1,2) + e2 * v2rho2_x(3*k-1)
dvxc_rr(k,2,1) = dvxc_rr(k,2,1) + e2 * v2rho2_x(3*k-1)
dvxc_rr(k,2,2) = dvxc_rr(k,2,2) + e2 * v2rho2_x(3*k)
ENDIF
dvxc_sr(k,1,1) = dvxc_sr(k,1,1) + e2 * v2rhosigma_x(6*k-5)*2._DP
DO k = 1, length
thr_up_cond = r_in(k,1)>epsr .AND. SQRT(ABS(sigma(3*k-2)))>epsg
thr_dw_cond = r_in(k,2)>epsr .AND. SQRT(ABS(sigma(3*k)))>epsg
IF ( .NOT.thr_up_cond .OR. .NOT.thr_dw_cond ) CYCLE
dvxc_rr(k,1,1) = dvxc_rr(k,1,1) + e2 * v2rho2_x(3*k-2)
dvxc_ss(k,1,1) = dvxc_ss(k,1,1) + e2 * v2sigma2_x(6*k-5)*4._DP
dvxc_sr(k,2,2) = dvxc_sr(k,2,2) + e2 * v2rhosigma_x(6*k)*2._DP
dvxc_rr(k,2,2) = dvxc_rr(k,2,2) + e2 * v2rho2_x(3*k)
dvxc_ss(k,2,2) = dvxc_ss(k,2,2) + e2 * v2sigma2_x(6*k)*4._DP
ENDDO
!
DEALLOCATE( v2rho2_x, v2rhosigma_x, v2sigma2_x )
!
ENDIF
!
!
IF ( is_libxc(4) ) THEN
!
DO k = 1, length
rht = r_in(k,1) + r_in(k,2)
IF (rht > epsr) THEN
dvxc_rr(k,1,1) = dvxc_rr(k,1,1) + e2 * v2rho2_c(3*k-2)
dvxc_rr(k,1,2) = dvxc_rr(k,1,2) + e2 * v2rho2_c(3*k-1)
dvxc_rr(k,2,1) = dvxc_rr(k,2,1) + e2 * v2rho2_c(3*k-1)
dvxc_rr(k,2,2) = dvxc_rr(k,2,2) + e2 * v2rho2_c(3*k)
ENDIF
dvxc_rr(k,1,2) = dvxc_rr(k,1,2) + e2 * v2rho2_x(3*k-1)
dvxc_sr(k,1,1) = dvxc_sr(k,1,1) + e2 * v2rhosigma_x(6*k-5)*2._DP
dvxc_rr(k,2,1) = dvxc_rr(k,2,1) + e2 * v2rho2_x(3*k-1)
dvxc_sr(k,2,2) = dvxc_sr(k,2,2) + e2 * v2rhosigma_x(6*k)*2._DP
ENDDO
!
DEALLOCATE( v2rho2_x, v2rhosigma_x, v2sigma2_x )
!
ENDIF
!
IF ( is_libxc(4) ) THEN
!
DO k = 1, length
thr_up_cond = r_in(k,1)>epsr .AND. SQRT(ABS(sigma(3*k-2)))>epsg
thr_dw_cond = r_in(k,2)>epsr .AND. SQRT(ABS(sigma(3*k)))>epsg
IF ( .NOT.thr_up_cond .OR. .NOT.thr_dw_cond ) CYCLE
dvxc_rr(k,1,1) = dvxc_rr(k,1,1) + e2 * v2rho2_c(3*k-2)
dvxc_rr(k,1,2) = dvxc_rr(k,1,2) + e2 * v2rho2_c(3*k-1)
dvxc_rr(k,2,1) = dvxc_rr(k,2,1) + e2 * v2rho2_c(3*k-1)
dvxc_rr(k,2,2) = dvxc_rr(k,2,2) + e2 * v2rho2_c(3*k)
dvxc_sr(k,1,1) = dvxc_sr(k,1,1) + e2 * v2rhosigma_c(6*k-5)*2.d0
dvxc_ss(k,1,1) = dvxc_ss(k,1,1) + e2 * v2sigma2_c(6*k)*4.d0
dvxc_sr(k,1,2) = dvxc_sr(k,1,2) + e2 * v2rhosigma_c(6*k-4)

25
XClib/xc_wrapper_d_lda_lsda.f90
View File

@ -94,7 +94,8 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
ALLOCATE( dmex_lxc(length_dlxc) )
! ... DERIVATIVE FOR EXCHANGE
dmex_lxc(:) = 0.0_DP
IF (iexch /= 0) THEN
IF (iexch /= 0) THEN
CALL xc_f03_func_set_dens_threshold( xc_func(1), rho_threshold_lda )
CALL xc_f03_lda_fxc( xc_func(1), lengthxc, rho_lxc(1), dmex_lxc(1) )
ENDIF
ENDIF
@ -105,6 +106,7 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
dmcr_lxc(:) = 0.0_DP
IF (icorr /= 0) THEN
fkind_x = xc_f03_func_info_get_kind( xc_info(2) )
CALL xc_f03_func_set_dens_threshold( xc_func(2), rho_threshold_lda )
CALL xc_f03_lda_fxc( xc_func(2), lengthxc, rho_lxc(1), dmcr_lxc(1) )
ENDIF
ENDIF
@ -127,13 +129,26 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
SELECT CASE( sr_d )
CASE( 1 )
!
IF ( is_libxc(1) ) dmuxc(:,1,1) = dmuxc(:,1,1) + dmex_lxc(:)*2.0_DP
IF ( is_libxc(2) ) dmuxc(:,1,1) = dmuxc(:,1,1) + dmcr_lxc(:)*2.0_DP
IF ( is_libxc(1) ) THEN
DO ir = 1, length
IF (rho_in(ir,1)<=rho_threshold_lda ) CYCLE
dmuxc(ir,1,1) = dmuxc(ir,1,1) + dmex_lxc(ir)*2.0_DP
ENDDO
ENDIF
!
IF ( is_libxc(2) ) THEN
DO ir = 1, length
IF (rho_in(ir,1)<=rho_threshold_lda ) CYCLE
dmuxc(ir,1,1) = dmuxc(ir,1,1) + dmcr_lxc(ir)*2.0_DP
ENDDO
ENDIF
!
CASE( 2 )
!
IF ( is_libxc(1) ) THEN
DO ir = 1, length
IF (rho_in(ir,1)<=rho_threshold_lda .OR. &
rho_in(ir,2)<=rho_threshold_lda) CYCLE
dmuxc(ir,1,1) = dmuxc(ir,1,1) + dmex_lxc(3*ir-2)*2.0_DP
dmuxc(ir,1,2) = dmuxc(ir,1,2) + dmex_lxc(3*ir-1)*2.0_DP
dmuxc(ir,2,1) = dmuxc(ir,2,1) + dmex_lxc(3*ir-1)*2.0_DP
@ -143,7 +158,9 @@ SUBROUTINE dmxc( length, sr_d, rho_in, dmuxc )
ENDIF
!
IF ( is_libxc(2) ) THEN
DO ir = 1, length
DO ir = 1, length
IF (rho_in(ir,1)<=rho_threshold_lda .OR. &
rho_in(ir,2)<=rho_threshold_lda) CYCLE
dmuxc(ir,1,1) = dmuxc(ir,1,1) + dmcr_lxc(3*ir-2)*2.0_DP
dmuxc(ir,1,2) = dmuxc(ir,1,2) + dmcr_lxc(3*ir-1)*2.0_DP
dmuxc(ir,2,1) = dmuxc(ir,2,1) + dmcr_lxc(3*ir-1)*2.0_DP

80
XClib/xc_wrapper_gga.f90
View File

@ -75,8 +75,7 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!
INTEGER :: k, is
REAL(DP) :: sgn(2)
REAL(DP) :: rho_up, rho_dw, grho2_up, grho2_dw
REAL(DP) :: xlda_up, xlda_dw, xgga_up, xgga_dw
REAL(DP) :: rho_up, rho_dw, grho_up, grho_dw
REAL(DP), PARAMETER :: small = 1.E-10_DP
!
!
@ -166,32 +165,28 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
IF ( rho_lxc(k) > rho_threshold_lda ) THEN
ec(k) = ec_lxc(k) * rho_lxc(k) * SIGN(1.0_DP, rho(k,1))
v1c(k,1) = vc_rho(k)
ENDIF
IF ( rho_lxc(k)>small .AND. SQRT(ABS(sigma(k)))>small ) &
v2c(k,1) = vc_sigma(k)*2.d0
IF ( rho_lxc(k) <= rho_threshold_lda ) CYCLE
ec(k) = ec_lxc(k) * rho_lxc(k) * SIGN(1.0_DP, rho(k,1))
v1c(k,1) = vc_rho(k)
IF ( rho_lxc(k) <= rho_threshold_gga .OR. &
SQRT(ABS(sigma(k))) <= grho_threshold_gga) CYCLE
v2c(k,1) = vc_sigma(k)*2.d0
ENDDO
ELSE
DO k = 1, length
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
grho_up = SQRT(ABS(sigma(3*k-2)))
grho_dw = SQRT(ABS(sigma(3*k)))
IF ( rho_up <= rho_threshold_lda .OR. rho_dw <= rho_threshold_lda ) CYCLE
ec(k) = ec_lxc(k) * (rho_up+rho_dw)
v1c(k,1) = vc_rho(2*k-1)
v1c(k,2) = vc_rho(2*k)
IF ( rho_up <= rho_threshold_gga .OR. rho_dw <= rho_threshold_gga .OR. &
grho_up<=grho_threshold_gga .OR. grho_dw<=grho_threshold_gga ) CYCLE
v2c(k,1) = vc_sigma(3*k-2)*2.d0
v2c_ud(k)= vc_sigma(3*k-1)
v2c(k,2) = vc_sigma(3*k)*2.d0
ENDDO
ENDIF
!
@ -256,8 +251,7 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!
IF ( is_libxc(3) ) THEN
!
!CALL xc_f03_func_set_dens_threshold( xc_func(3), rho_threshold_gga )
CALL xc_f03_func_set_dens_threshold( xc_func(3), small )
CALL xc_f03_func_set_dens_threshold( xc_func(3), grho_threshold_gga )
IF (libxc_flags(3,0)==1) THEN
CALL xc_f03_gga_exc_vxc( xc_func(3), lengthxc, rho_lxc(1), sigma(1), ex_lxc(1), vx_rho(1), vx_sigma(1) )
ELSE
@ -267,31 +261,27 @@ SUBROUTINE xc_gcx_( length, ns, rho, grho, ex, ec, v1x, v2x, v1c, v2c, v2c_ud )
!
IF (.NOT. POLARIZED) THEN
DO k = 1, length
xlda_up = 1.0_DP ; xgga_up = 1.0_DP
IF ( rho_lxc(k) <= rho_threshold_lda ) xlda_up = 0.0_DP
IF ( rho_lxc(k)<=small .OR. SQRT(ABS(sigma(k)))<=small ) xgga_up = 0.0_DP
ex(k) = ex_lxc(k) * rho_lxc(k) * SIGN(1.0_DP, rho(k,1)) *xlda_up
v1x(k,1) = vx_rho(k)*xlda_up
v2x(k,1) = vx_sigma(k)*2.d0*xgga_up
IF ( rho_lxc(k) <= rho_threshold_lda ) CYCLE
ex(k) = ex_lxc(k) * rho_lxc(k) * SIGN(1.0_DP, rho(k,1))
v1x(k,1) = vx_rho(k)
IF ( rho_lxc(k) <= rho_threshold_gga .OR. &
SQRT(ABS(sigma(k))) <= grho_threshold_gga) CYCLE
v2x(k,1) = vx_sigma(k)*2.d0
ENDDO
ELSE
DO k = 1, length
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 <= small ) 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
ex(k) = ex_lxc(k) * (rho_up*xlda_up+rho_dw*xlda_dw)
v1x(k,1) = vx_rho(2*k-1)*xlda_up
v1x(k,2) = vx_rho(2*k)*xlda_dw
v2x(k,1) = vx_sigma(3*k-2)*2.d0*xgga_up
v2x(k,2) = vx_sigma(3*k)*2.d0*xgga_dw
grho_up = SQRT(ABS(sigma(3*k-2)))
grho_dw = SQRT(ABS(sigma(3*k)))
IF ( rho_up <= rho_threshold_lda .OR. rho_dw <= rho_threshold_lda ) CYCLE
ex(k) = ex_lxc(k) * (rho_up+rho_dw)
v1x(k,1) = vx_rho(2*k-1)
v1x(k,2) = vx_rho(2*k)
IF ( rho_up <= rho_threshold_gga .OR. rho_dw <= rho_threshold_gga .OR. &
grho_up<=grho_threshold_gga .OR. grho_dw<=grho_threshold_gga ) CYCLE
v2x(k,1) = vx_sigma(3*k-2)*2.d0
v2x(k,2) = vx_sigma(3*k)*2.d0
ENDDO
ENDIF
!

10
XClib/xclib_test.f90
View File

@ -1592,8 +1592,8 @@ PROGRAM xclib_test
ENDIF
!
IF (mype==root .AND. TRIM(test)=='EXECUTE') THEN
CALL print_aver( what, xc_aver(:,1), aver(1), 'up' )
CALL print_aver( what, xc_aver(:,2), aver(2), 'down' )
CALL print_aver( what, xc_aver(:,1), aver(1), ' up' )
CALL print_aver( what, xc_aver(:,2), aver(2), ' down' )
ENDIF
!
ENDIF
@ -1659,9 +1659,9 @@ PROGRAM xclib_test
#endif
!
IF (mype==root .AND. TRIM(test)=='EXECUTE') THEN
CALL print_aver( what, dxc_aver(:,1), aver(1), 'up-up' )
CALL print_aver( what, dxc_aver(:,2), aver(2), 'up-down' )
CALL print_aver( what, dxc_aver(:,3), aver(3), 'down-down' )
CALL print_aver( what, dxc_aver(:,1), aver(1), ' up-up' )
CALL print_aver( what, dxc_aver(:,2), aver(2), ' up-down' )
CALL print_aver( what, dxc_aver(:,3), aver(3), ' down-down' )
ENDIF
!
ENDIF