mirror of https://gitlab.com/QEF/q-e.git
XClib-acc - GGA exch pol
This commit is contained in:
parent
ae37e1c9a4
commit
7727978d0a
|
@ -423,7 +423,7 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
|
|||
!! Up and down charge density
|
||||
REAL(DP), INTENT(IN), DIMENSION(length,2) :: grho2_in
|
||||
!! Up and down gradient of the charge
|
||||
REAL(DP), INTENT(OUT), DIMENSION(length) :: sx_tot
|
||||
REAL(DP), INTENT(OUT), DIMENSION(length) :: sx_tot
|
||||
!! Energy exchange GGA
|
||||
REAL(DP), INTENT(OUT), DIMENSION(length,2) :: v1x_out
|
||||
!! Exchange potential (density part)
|
||||
|
@ -432,17 +432,18 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
|
|||
!
|
||||
! ... local variables
|
||||
!
|
||||
INTEGER :: ir, is, iflag
|
||||
REAL(DP) :: rho(2), grho2(2)
|
||||
REAL(DP) :: v1x(2), v2x(2)
|
||||
REAL(DP) :: sx(2), rnull(2)
|
||||
REAL(DP) :: sxsr(2)
|
||||
REAL(DP) :: v1xsr(2), v2xsr(2)
|
||||
INTEGER :: ir, is, iflag
|
||||
REAL(DP) :: rho_up, rho_dw, grho2_up, grho2_dw
|
||||
REAL(DP) :: v1x_up, v1x_dw, v2x_up, v2x_dw
|
||||
REAL(DP) :: sx_up, sx_dw, rnull_up, rnull_dw
|
||||
REAL(DP) :: sxsr_up, sxsr_dw
|
||||
REAL(DP) :: v1xsr_up, v1xsr_dw, v2xsr_up, v2xsr_dw
|
||||
!
|
||||
REAL(DP), PARAMETER :: small=1.D-10
|
||||
REAL(DP), PARAMETER :: rho_trash=0.5_DP, grho2_trash=0.2_DP
|
||||
! temporary values assigned to rho and grho when they
|
||||
! are too small in order to avoid numerical problems.
|
||||
!
|
||||
#if defined(_OPENMP)
|
||||
INTEGER :: ntids
|
||||
INTEGER, EXTERNAL :: omp_get_num_threads
|
||||
|
@ -452,60 +453,78 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
|
|||
!
|
||||
sx_tot = 0.0_DP
|
||||
!
|
||||
!
|
||||
#if defined(_OPENACC)
|
||||
!$acc data copyin(rho_in, grho2_in), copyout(sx_tot, v1x_out, v2x_out)
|
||||
!$acc parallel loop
|
||||
#endif
|
||||
#if defined(__OPENMP) && !defined(_OPENACC)
|
||||
!$omp parallel if(ntids==1) default(none) &
|
||||
!$omp private( rho, grho2, rnull, sx, sxsr, v1x, v1xsr, &
|
||||
!$omp v2x, v2xsr, iflag ) &
|
||||
!$omp shared(rho_in, length, grho2_in, sx_tot, v1x_out, v2x_out,&
|
||||
!$omp igcx, exx_started, exx_fraction, screening_parameter, gau_parameter)
|
||||
!$omp private( rho_up, rho_dw, grho2_up, grho2_dw, rnull_up, rnull_dw, &
|
||||
!$omp sx_up, sx_dw, sxsr_up, sxsr_dw, v1xsr_up, v1xsr_dw, &
|
||||
!$omp v1x_up, v1x_dw, v2x_up, v2x_dw, v2xsr_up, v2xsr_dw, &
|
||||
!$omp iflag ) &
|
||||
!$omp shared( rho_in, length, grho2_in, sx_tot, v1x_out, v2x_out, &
|
||||
!$omp igcx, exx_started, exx_fraction, screening_parameter,&
|
||||
!$omp gau_parameter )
|
||||
!$omp do
|
||||
#endif
|
||||
!
|
||||
DO ir = 1, length
|
||||
!
|
||||
rho(:) = rho_in(ir,:)
|
||||
grho2(:) = grho2_in(ir,:)
|
||||
rnull(:) = 1.0_DP
|
||||
rho_up = rho_in(ir,1)
|
||||
rho_dw = rho_in(ir,2)
|
||||
grho2_up = grho2_in(ir,1)
|
||||
grho2_dw = grho2_in(ir,2)
|
||||
rnull_up = 1.0_DP
|
||||
rnull_dw = 1.0_DP
|
||||
!
|
||||
IF ( rho(1)+rho(2) <= small ) THEN
|
||||
IF ( rho_up+rho_dw <= small ) THEN
|
||||
sx_tot(ir) = 0.0_DP
|
||||
v1x_out(ir,:) = 0.0_DP ; v2x_out(ir,:) = 0.0_DP
|
||||
v1x_out(ir,1) = 0.0_DP
|
||||
v2x_out(ir,1) = 0.0_DP
|
||||
v1x_out(ir,2) = 0.0_DP
|
||||
v2x_out(ir,2) = 0.0_DP
|
||||
CYCLE
|
||||
ELSE
|
||||
DO is = 1, 2
|
||||
IF ( rho(is)<=small .OR. SQRT(ABS(grho2(is)))<=small ) THEN
|
||||
rho(is) = rho_trash
|
||||
grho2(is) = grho2_trash
|
||||
rnull(is) = 0.0_DP
|
||||
ENDIF
|
||||
ENDDO
|
||||
IF ( rho_up<=small .OR. SQRT(ABS(grho2_up))<=small ) THEN
|
||||
rho_up = rho_trash
|
||||
grho2_up = grho2_trash
|
||||
rnull_up = 0.0_DP
|
||||
ENDIF
|
||||
IF ( rho_dw<=small .OR. SQRT(ABS(grho2_dw))<=small ) THEN
|
||||
rho_dw = rho_trash
|
||||
grho2_dw = grho2_trash
|
||||
rnull_dw = 0.0_DP
|
||||
ENDIF
|
||||
ENDIF
|
||||
!
|
||||
!
|
||||
! ... exchange
|
||||
!
|
||||
SELECT CASE( igcx )
|
||||
CASE( 0 )
|
||||
!
|
||||
sx_tot(ir) = 0.0_DP
|
||||
v1x = 0.0_DP
|
||||
v2x = 0.0_DP
|
||||
v1x_up = 0.0_DP ; v1x_dw = 0.0_DP
|
||||
v2x_up = 0.0_DP ; v2x_dw = 0.0_DP
|
||||
!
|
||||
CASE( 1 )
|
||||
!
|
||||
CALL becke88_spin( rho(1), rho(2), grho2(1), grho2(2), sx(1), sx(2), &
|
||||
v1x(1), v1x(2), v2x(1), v2x(2) )
|
||||
CALL becke88_spin( rho_up, rho_dw, grho2_up, grho2_dw, sx_up, sx_dw, &
|
||||
v1x_up, v1x_dw, v2x_up, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = sx(1)*rnull(1) + sx(2)*rnull(2)
|
||||
sx_tot(ir) = sx_up*rnull_up + sx_dw*rnull_dw
|
||||
!
|
||||
CASE( 2 )
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL ggax( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL ggax( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL ggax( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL ggax( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 3, 4, 8, 10, 12, 20, 23, 24, 25, 44, 45 )
|
||||
! igcx=3: PBE, igcx=4: revised PBE, igcx=8: PBE0, igcx=10: PBEsol
|
||||
|
@ -521,272 +540,311 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
|
|||
IF ( igcx==44 ) iflag = 8
|
||||
IF ( igcx==45 ) iflag = 9
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL pbex( rho(1), grho2(1), iflag, sx(1), v1x(1), v2x(1) )
|
||||
CALL pbex( rho(2), grho2(2), iflag, sx(2), v1x(2), v2x(2) )
|
||||
CALL pbex( rho_up, grho2_up, iflag, sx_up, v1x_up, v2x_up )
|
||||
CALL pbex( rho_dw, grho2_dw, iflag, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( igcx == 8 .AND. exx_started ) THEN
|
||||
!
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_dw = (1.0_DP - exx_fraction) * v1x_dw
|
||||
v2x_up = (1.0_DP - exx_fraction) * v2x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
!
|
||||
ELSEIF ( igcx == 12 .AND. exx_started ) THEN
|
||||
!
|
||||
CALL pbexsr( rho(1), grho2(1), sxsr(1), v1xsr(1), &
|
||||
v2xsr(1), screening_parameter )
|
||||
CALL pbexsr( rho(2), grho2(2), sxsr(2), v1xsr(2), &
|
||||
v2xsr(2), screening_parameter )
|
||||
CALL pbexsr( rho_up, grho2_up, sxsr_up, v1xsr_up, &
|
||||
v2xsr_up, screening_parameter )
|
||||
CALL pbexsr( rho_dw, grho2_dw, sxsr_dw, v1xsr_dw, &
|
||||
v2xsr_dw, screening_parameter )
|
||||
!
|
||||
sx_tot(ir) = sx_tot(ir) - exx_fraction*0.5_DP*( sxsr(1)*rnull(1) + &
|
||||
sxsr(2)*rnull(2) )
|
||||
v1x = v1x - exx_fraction * v1xsr
|
||||
v2x = v2x - exx_fraction * v2xsr * 2.0_DP
|
||||
sx_tot(ir) = sx_tot(ir) - exx_fraction*0.5_DP*( sxsr_up*rnull_up + &
|
||||
sxsr_dw*rnull_dw )
|
||||
v1x_up = v1x_up - exx_fraction * v1xsr_up
|
||||
v1x_dw = v1x_dw - exx_fraction * v1xsr_dw
|
||||
v2x_up = v2x_up - exx_fraction * v2xsr_up * 2.0_DP
|
||||
v2x_dw = v2x_dw - exx_fraction * v2xsr_dw * 2.0_DP
|
||||
!
|
||||
ELSEIF ( igcx == 20 .AND. exx_started ) THEN
|
||||
! gau-pbe
|
||||
!CALL pbexgau_lsd( rho, grho2, sxsr, v1xsr, v2xsr, gau_parameter )
|
||||
CALL pbexgau( rho(1),grho2(1), sxsr(1), v1xsr(1),v2xsr(1), gau_parameter )
|
||||
CALL pbexgau( rho(2),grho2(2), sxsr(2), v1xsr(2),v2xsr(2), gau_parameter )
|
||||
CALL pbexgau( rho_up,grho2_up, sxsr_up, v1xsr_up,v2xsr_up, gau_parameter )
|
||||
CALL pbexgau( rho_dw,grho2_dw, sxsr_dw, v1xsr_dw,v2xsr_dw, gau_parameter )
|
||||
!
|
||||
sx_tot(ir) = sx_tot(ir) - exx_fraction*0.5_DP * ( sxsr(1)*rnull(1) + &
|
||||
sxsr(2)*rnull(2) )
|
||||
v1x = v1x - exx_fraction * v1xsr
|
||||
v2x = v2x - exx_fraction * v2xsr * 2.0_DP
|
||||
sx_tot(ir) = sx_tot(ir) - exx_fraction*0.5_DP * ( sxsr_up*rnull_up + &
|
||||
sxsr_dw*rnull_dw )
|
||||
v1x_up = v1x_up - exx_fraction * v1xsr_up
|
||||
v1x_dw = v1x_dw - exx_fraction * v1xsr_dw
|
||||
v2x_up = v2x_up - exx_fraction * v2xsr_up * 2.0_DP
|
||||
v2x_dw = v2x_dw - exx_fraction * v2xsr_dw * 2.0_DP
|
||||
!
|
||||
ENDIF
|
||||
!
|
||||
CASE( 9 ) ! B3LYP
|
||||
!
|
||||
CALL becke88_spin( rho(1), rho(2), grho2(1), grho2(2), sx(1), sx(2), &
|
||||
v1x(1), v1x(2), v2x(1), v2x(2) )
|
||||
CALL becke88_spin( rho_up, rho_dw, grho2_up, grho2_dw, sx_up, sx_dw, &
|
||||
v1x_up, v1x_dw, v2x_up, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = sx(1)*rnull(1) + sx(2)*rnull(2)
|
||||
sx_tot(ir) = sx_up*rnull_up + sx_dw*rnull_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot(ir) = 0.72_DP * sx_tot(ir)
|
||||
v1x = 0.72_DP * v1x
|
||||
v2x = 0.72_DP * v2x
|
||||
v1x_up = 0.72_DP * v1x_up ; v1x_dw = 0.72_DP * v1x_dw
|
||||
v2x_up = 0.72_DP * v2x_up ; v2x_dw = 0.72_DP * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 11 ) ! 'Wu-Cohen'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL wcx( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL wcx( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL wcx( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL wcx( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 13 ) ! 'revised PW86 for vdw-df2'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL rPW86( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL rPW86( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL rPW86( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL rPW86( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 16 ) ! 'c09x for vdw-df-c09.'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL c09x( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL c09x( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL c09x( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL c09x( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 21 ) ! 'PW86'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL pw86( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL pw86( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL pw86( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL pw86( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 22 ) ! 'B86B'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL becke86b( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL becke86b( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL becke86b( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL becke86b( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 26, 46 ) ! 'B86R for rev-vdW-DF2'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
IF ( igcx==26 ) iflag = 3 ! B86R for rev-vdW-DF2
|
||||
IF ( igcx==46 ) iflag = 4 ! W32X for vdW-DF3-opt2
|
||||
CALL b86b( rho(1), grho2(1), iflag, sx(1), v1x(1), v2x(1) )
|
||||
CALL b86b( rho(2), grho2(2), iflag, sx(2), v1x(2), v2x(2) )
|
||||
CALL b86b( rho_up, grho2_up, iflag, sx_up, v1x_up, v2x_up )
|
||||
CALL b86b( rho_dw, grho2_dw, iflag, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 27 ) ! 'cx13 for vdw-df-cx'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL cx13( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL cx13( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL cx13( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL cx13( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
CASE( 28 ) ! X3LYP
|
||||
!
|
||||
CALL becke88_spin( rho(1), rho(2), grho2(1), grho2(2), sx(1), sx(2), &
|
||||
v1x(1), v1x(2), v2x(1), v2x(2) )
|
||||
CALL becke88_spin( rho_up, rho_dw, grho2_up, grho2_dw, sx_up, sx_dw, &
|
||||
v1x_up, v1x_dw, v2x_up, v2x_dw )
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up
|
||||
rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up
|
||||
grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL pbex( rho(1), grho2(1), 1, sxsr(1), v1xsr(1), v2xsr(1) )
|
||||
CALL pbex( rho(2), grho2(2), 1, sxsr(2), v1xsr(2), v2xsr(2) )
|
||||
CALL pbex( rho_up, grho2_up, 1, sxsr_up, v1xsr_up, v2xsr_up )
|
||||
CALL pbex( rho_dw, grho2_dw, 1, sxsr_dw, v1xsr_dw, v2xsr_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP*( sxsr(1)*rnull(1) + sxsr(2)*rnull(2) )*0.235_DP + &
|
||||
( sx(1)*rnull(1) + sx(2)*rnull(2) )*0.765_DP
|
||||
v1x = v1xsr * 0.235_DP + v1x * 0.765_DP
|
||||
v2x = v2xsr * 0.235_DP * 2.0_DP + v2x * 0.765_DP
|
||||
sx_tot(ir) = 0.5_DP*( sxsr_up*rnull_up + sxsr_dw*rnull_dw )*0.235_DP + &
|
||||
( sx_up*rnull_up + sx_dw*rnull_dw )*0.765_DP
|
||||
v1x_up = v1xsr_up * 0.235_DP + v1x_up * 0.765_DP
|
||||
v1x_dw = v1xsr_dw * 0.235_DP + v1x_dw * 0.765_DP
|
||||
v2x_up = v2xsr_up * 0.235_DP * 2.0_DP + v2x_up * 0.765_DP
|
||||
v2x_dw = v2xsr_dw * 0.235_DP * 2.0_DP + v2x_dw * 0.765_DP
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot(ir) = 0.709_DP * sx_tot(ir)
|
||||
v1x = 0.709_DP * v1x
|
||||
v2x = 0.709_DP * v2x
|
||||
v1x_up = 0.709_DP * v1x_up
|
||||
v1x_dw = 0.709_DP * v1x_dw
|
||||
v2x_up = 0.709_DP * v2x_up
|
||||
v2x_dw = 0.709_DP * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 29, 31 ) ! 'cx0 for vdw-df-cx0' or `cx0p for vdW-DF-cx0p'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL cx13( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL cx13( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL cx13( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL cx13( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 30 ) ! 'R860' = 'rPW86-0' for vdw-df2-0'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL rPW86( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL rPW86( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL rPW86( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL rPW86( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_dw = (1.0_DP - exx_fraction) * v1x_dw
|
||||
v2x_up = (1.0_DP - exx_fraction) * v2x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 38 ) ! 'br0 for vdw-df2-BR0' etc
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL b86b( rho(1), grho2(1), 3, sx(1), v1x(1), v2x(1) )
|
||||
CALL b86b( rho(2), grho2(2), 3, sx(2), v1x(2), v2x(2) )
|
||||
CALL b86b( rho_up, grho2_up, 3, sx_up, v1x_up, v2x_up )
|
||||
CALL b86b( rho_dw, grho2_dw, 3, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_dw = (1.0_DP - exx_fraction) * v1x_dw
|
||||
v2x_up = (1.0_DP - exx_fraction) * v2x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 40 ) ! 'c090 for vdw-df-c090' etc
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL c09x( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL c09x( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL c09x( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL c09x( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_dw = (1.0_DP - exx_fraction) * v1x_dw
|
||||
v2x_up = (1.0_DP - exx_fraction) * v2x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 41 ) ! B86X for B86BPBEX hybrid
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL becke86b( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL becke86b( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL becke86b( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL becke86b( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot = (1.0_DP - exx_fraction) * sx_tot
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_dw = (1.0_DP - exx_fraction) * v1x_dw
|
||||
v2x_up = (1.0_DP - exx_fraction) * v2x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 42 ) ! B88X for BHANDHLYP
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw
|
||||
grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
!
|
||||
CALL becke88( rho(1), grho2(1), sx(1), v1x(1), v2x(1) )
|
||||
CALL becke88( rho(2), grho2(2), sx(2), v1x(2), v2x(2) )
|
||||
CALL becke88( rho_up, grho2_up, sx_up, v1x_up, v2x_up )
|
||||
CALL becke88( rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw )
|
||||
!
|
||||
sx_tot = 0.5_DP * ( sx(1)*rnull(1) + sx(2)*rnull(2) )
|
||||
v2x = 2.0_DP * v2x
|
||||
sx_tot(ir) = 0.5_DP * ( sx_up*rnull_up + sx_dw*rnull_dw )
|
||||
v2x_up = 2.0_DP * v2x_up
|
||||
v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
IF ( exx_started ) THEN
|
||||
sx_tot = (1.0_DP - exx_fraction) * sx_tot
|
||||
v1x = (1.0_DP - exx_fraction) * v1x
|
||||
v2x = (1.0_DP - exx_fraction) * v2x
|
||||
sx_tot(ir) = (1.0_DP - exx_fraction) * sx_tot(ir)
|
||||
v1x_up = (1.0_DP - exx_fraction) * v1x_up
|
||||
v1x_dw = (1.0_DP - exx_fraction) * v1x_dw
|
||||
v2x_up = (1.0_DP - exx_fraction) * v2x_up
|
||||
v2x_dw = (1.0_DP - exx_fraction) * v2x_dw
|
||||
ENDIF
|
||||
!
|
||||
CASE( 43 ) ! 'beefx'
|
||||
!
|
||||
rho = 2.0_DP * rho
|
||||
grho2 = 4.0_DP * grho2
|
||||
!
|
||||
CALL beefx(rho(1), grho2(1), sx(1), v1x(1), v2x(1), 0)
|
||||
CALL beefx(rho(2), grho2(2), sx(2), v1x(2), v2x(2), 0)
|
||||
!
|
||||
sx_tot(ir) = 0.5_DP * (sx(1)*rnull(1) + sx(2)*rnull(2))
|
||||
v2x = 2.0_DP * v2x
|
||||
! CASE( 43 ) ! 'beefx'
|
||||
! !
|
||||
! rho_up = 2.0_DP * rho_up ; rho_dw = 2.0_DP * rho_dw !*****TEMPORARY OUT -- openACC TEST
|
||||
! grho2_up = 4.0_DP * grho2_up ; grho2_dw = 4.0_DP * grho2_dw
|
||||
! !
|
||||
! CALL beefx(rho_up, grho2_up, sx_up, v1x_up, v2x_up, 0)
|
||||
! CALL beefx(rho_dw, grho2_dw, sx_dw, v1x_dw, v2x_dw, 0)
|
||||
! !
|
||||
! sx_tot(ir) = 0.5_DP * (sx_up*rnull_up + sx_dw*rnull_dw)
|
||||
! v2x_up = 2.0_DP * v2x_up
|
||||
! v2x_dw = 2.0_DP * v2x_dw
|
||||
!
|
||||
! case igcx == 5 (HCTH) and 6 (OPTX) not implemented
|
||||
! case igcx == 7 (meta-GGA) must be treated in a separate call to another
|
||||
|
@ -794,18 +852,25 @@ SUBROUTINE gcx_spin( length, rho_in, grho2_in, sx_tot, v1x_out, v2x_out )
|
|||
!
|
||||
CASE DEFAULT
|
||||
!
|
||||
sx = 0.0_DP
|
||||
v1x = 0.0_DP
|
||||
v2x = 0.0_DP
|
||||
sx_tot(ir) = 0.0_DP
|
||||
v1x_up = 0.0_DP ; v1x_dw = 0.0_DP
|
||||
v2x_up = 0.0_DP ; v2x_dw = 0.0_DP
|
||||
!
|
||||
END SELECT
|
||||
!
|
||||
v1x_out(ir,:) = v1x(:) * rnull(:)
|
||||
v2x_out(ir,:) = v2x(:) * rnull(:)
|
||||
v1x_out(ir,1) = v1x_up * rnull_up
|
||||
v1x_out(ir,2) = v1x_dw * rnull_dw
|
||||
v2x_out(ir,1) = v2x_up * rnull_up
|
||||
v2x_out(ir,2) = v2x_dw * rnull_dw
|
||||
!
|
||||
ENDDO
|
||||
#if defined(_OPENACC)
|
||||
!$acc end data
|
||||
#endif
|
||||
#if defined(__OPENMP) && !defined(_OPENACC)
|
||||
!$omp end do
|
||||
!$omp end parallel
|
||||
#endif
|
||||
!
|
||||
!
|
||||
RETURN
|
||||
|
|
|
@ -1219,7 +1219,8 @@ END SUBROUTINE cx13
|
|||
! ===========> SPIN <===========
|
||||
!
|
||||
!-----------------------------------------------------------------------
|
||||
SUBROUTINE becke88_spin( rho_up, rho_dw, grho_up, grho_dw, sx_up, sx_dw, v1x_up, v1x_dw, v2x_up, v2x_dw ) !<GPU:DEVICE>
|
||||
SUBROUTINE becke88_spin( rho_up, rho_dw, grho_up, grho_dw, sx_up, sx_dw, v1x_up, v1x_dw, v2x_up, v2x_dw )
|
||||
!$acc routine (becke88_spin) seq
|
||||
!-----------------------------------------------------------------------
|
||||
!! Becke exchange: A.D. Becke, PRA 38, 3098 (1988) - Spin polarized case
|
||||
!
|
||||
|
|
Loading…
Reference in New Issue