mirror of https://gitlab.com/QEF/q-e.git
beta smoothing: fix filter definition and enforce estimated erors in R and Q spaces. format printed summary.
This commit is contained in:
parent
8ca7f15e77
commit
992a4cf448
|
@ -13,7 +13,7 @@ SUBROUTINE init_us_b0
|
|||
USE kinds, ONLY : DP
|
||||
USE gvecw, ONLY : ecutwfc
|
||||
USE io_global, ONLY : stdout
|
||||
USE constants, ONLY : fpi, sqrt2, eps8
|
||||
USE constants, ONLY : fpi
|
||||
USE atom, ONLY : rgrid
|
||||
USE ions_base, ONLY : ntyp => nsp
|
||||
USE us, ONLY : dq
|
||||
|
@ -23,15 +23,14 @@ SUBROUTINE init_us_b0
|
|||
!
|
||||
IMPLICIT NONE
|
||||
!
|
||||
LOGICAL, PARAMETER :: tprint=.FALSE. ! Whether the beta_l(r) and
|
||||
! its relatives are printed or not.
|
||||
INTEGER, PARAMETER :: nn=16 ! Smoothing parameter, order of the polynomial
|
||||
! inverse gaussian approximant.
|
||||
REAL(DP), PARAMETER :: a=22.d0 ! Smoothing parameter, exponent of the gaussian
|
||||
! decaying factor.
|
||||
REAL(DP) :: rcut, drcut ! beta function cutoff radius and estimated increase of
|
||||
! it due to the filtering
|
||||
REAL(DP), PARAMETER :: eps = 5.d-9
|
||||
LOGICAL, PARAMETER :: tprint=.FALSE. ! Whether the beta_l(r) and its relatives are printed or not.
|
||||
INTEGER, PARAMETER :: nn=20 ! Smoothing parameter, order of the polynomial in the
|
||||
! inverse gaussian approximant.
|
||||
REAL(DP), PARAMETER:: a=1.125d0*(2*nn+1) ! Smoothing parameter, exponent of the gaussian
|
||||
! decaying factor, chosen so that the flex is at x=2/3.
|
||||
REAL(DP) :: rcut, drcut ! beta function cutoff radius and its estimated increase due to the filtering
|
||||
|
||||
REAL(DP), PARAMETER :: eps = 1.d-9 ! error tollerance for intergrals, norms etc.
|
||||
!
|
||||
INTEGER :: nqx
|
||||
REAL(DP), ALLOCATABLE :: tab0(:,:), tab(:,:), beta(:,:), betas(:,:)
|
||||
|
@ -39,15 +38,15 @@ SUBROUTINE init_us_b0
|
|||
!
|
||||
INTEGER :: nt, nb, l, ir, iq, startq, lastq, ndm
|
||||
! various counters
|
||||
REAL(DP), ALLOCATABLE :: aux(:), besr(:)
|
||||
REAL(DP), ALLOCATABLE :: aux(:), besr(:), ffrr(:)
|
||||
! various work space
|
||||
REAL(DP) :: q, qi, qmax
|
||||
REAL(DP) :: qi, qmax, vqint
|
||||
! the modulus of g for each shell
|
||||
! q-point grid for interpolation
|
||||
REAL(DP) :: vqint
|
||||
! interpolated value
|
||||
INTEGER :: kktest
|
||||
REAL(DP) :: test, test_save
|
||||
REAL(DP) :: test, safe_test, error_estimate
|
||||
LOGICAL :: first
|
||||
!
|
||||
CHARACTER(LEN=4) :: filename
|
||||
!
|
||||
|
@ -55,14 +54,19 @@ SUBROUTINE init_us_b0
|
|||
!
|
||||
! ... Initialization of variables
|
||||
!
|
||||
drcut = ABS(LOG(eps8))/SQRT(ecutwfc)/2.d0
|
||||
qmax = 3.d0 * SQRT(ecutwfc)
|
||||
nqx = INT( qmax / dq + 4) ! Think about what happens in a variable cell calculations
|
||||
drcut = 3.D0 * ABS(LOG(eps)) / SQRT(ecutwfc)/ 2.D0 ! estimate (to be restricted later)
|
||||
qmax = 3.D0 * SQRT(ecutwfc) ! maximum q that avoids aliasing, higher Fourier components must be negligible
|
||||
nqx = INT( qmax / dq + 4) ! Think about what happens in a variable cell calculations
|
||||
!
|
||||
ALLOCATE( tab0(nqx,nbetam), tab(nqx,nbetam) )
|
||||
ALLOCATE( power_r(nbetam), power_q(nbetam) )
|
||||
!
|
||||
IF (tprint) WRITE( stdout, * ) 'upf(nt)%kkbeta ', upf(1:ntyp)%kkbeta
|
||||
WRITE( stdout, '(/5X,"PSEUDOPOTENTIAL Beta Smoothing: ==============================================")' )
|
||||
IF (tprint) THEN
|
||||
WRITE( stdout,'(/5X,"table grid dimensions: NQX =",I7," NBETAM =", I4)'), nqx, nbetam
|
||||
WRITE( stdout,'(11X,"initial rcut indices:",10i7)') upf(1:ntyp)%kkbeta
|
||||
WRITE( stdout,'(19X,"rcut (ityp) :",10f7.3)'), (rgrid(nt)%r(upf(nt)%kkbeta),nt=1,ntyp)
|
||||
END IF
|
||||
!
|
||||
DO nt=1,ntyp
|
||||
rcut = rgrid(nt)%r(upf(nt)%kkbeta)
|
||||
|
@ -73,25 +77,20 @@ SUBROUTINE init_us_b0
|
|||
!
|
||||
ndm = MAXVAL( upf(:)%kkbeta )
|
||||
ALLOCATE( beta(ndm,nbetam), betas(ndm,nbetam) )
|
||||
ALLOCATE( aux(ndm),besr(ndm) )
|
||||
ALLOCATE( aux(ndm), besr(ndm), ffrr(ndm) )
|
||||
!
|
||||
WRITE(6,'(a)') 'Pseudopotential projectors are smoothed with filter( 1.2d0 * q / qmax, a, nn )'
|
||||
WRITE(6,'(a,f6.2,a,i4,4(a,f11.8))') 'FILTER : a=',a,', nn=',nn,', filter(1.1d0)=', filter(1.1d0,a,nn), &
|
||||
', filter(1.1d0/3)=', filter(1.1d0/3,a,nn),&
|
||||
', filter(1.2d0)=', filter(1.2d0,a,nn), &
|
||||
', filter(1.2d0/3)=', filter(1.2d0/3,a,nn)
|
||||
IF (tprint) THEN
|
||||
WRITE( stdout, * ) " PSEUDOPOTENTIAL REPORT; initial values "
|
||||
WRITE( stdout, * ) ' NQX :', nqx, ' NBETAM :', nbetam
|
||||
WRITE( stdout, * ) ' NDM :', ndm, ' KKBETA :', upf(1:ntyp)%kkbeta
|
||||
WRITE( stdout, * ) ' r(ndm):', (rgrid(nt)%r(upf(nt)%kkbeta),nt=1,ntyp), ' drcut :', drcut
|
||||
ENDIF
|
||||
WRITE(6,'(/5X,a)') 'Pseudopotential projectors are smoothed with filter( q / qmax, a, nn )'
|
||||
WRITE(6,'(5X,a,f6.2,a,i4,2(a,f11.8))') 'FILTER : a=',a,', nn=',nn,', filter(1.0/3)=', filter(1.d0/3,a,nn), &
|
||||
', filter(1.0)=', filter(1.d0,a,nn)
|
||||
WRITE(6,'(14X,a,1p,e12.4)') 'target error: eps =', eps
|
||||
!
|
||||
! ... fill the interpolation table tab
|
||||
!
|
||||
CALL divide( intra_bgrp_comm, nqx, startq, lastq )
|
||||
!
|
||||
!- loop over pseudopotentials
|
||||
DO nt = 1, ntyp
|
||||
! WRITE( stdout, * ) ' ityp = ', nt
|
||||
WRITE( stdout, '(/5X,a,i4)' ) 'Smoothing PSEUDO #', nt
|
||||
IF ( upf(nt)%is_gth ) CYCLE
|
||||
!
|
||||
tab0(:,:) = 0.d0 ; tab(:,:) = 0.d0
|
||||
|
@ -108,6 +107,17 @@ SUBROUTINE init_us_b0
|
|||
ENDDO
|
||||
CLOSE (4)
|
||||
ENDIF
|
||||
!- compute original beta normalization in real space where grid is sufficient by definition
|
||||
DO nb = 1, upf(nt)%nbeta
|
||||
l = upf(nt)%lll (nb)
|
||||
kktest = upf(nt)%kkbeta
|
||||
aux(1:kktest) = upf(nt)%beta(1:kktest,nb) * upf(nt)%beta(1:kktest,nb)
|
||||
CALL simpson( kktest, aux, rgrid(nt)%rab, power_r(nb) )
|
||||
ENDDO
|
||||
IF (tprint) THEN
|
||||
WRITE( stdout, '(5X,a)' ) "BETA FUNCTIONS NORMS:"
|
||||
WRITE( stdout, '(5X,a,1p,10e12.4)'), "norm2 R:", (power_r(nb), nb=1,upf(nt)%nbeta)
|
||||
END IF
|
||||
!-
|
||||
DO nb = 1, upf(nt)%nbeta
|
||||
l = upf(nt)%lll (nb)
|
||||
|
@ -120,27 +130,42 @@ SUBROUTINE init_us_b0
|
|||
CALL simpson( upf(nt)%kkbeta, aux, rgrid(nt)%rab, vqint )
|
||||
!
|
||||
tab0(iq, nb) = vqint
|
||||
! tab (iq, nb) = vqint * filter( 1.1d0 * qi / qmax, a, nn )
|
||||
tab (iq, nb) = vqint * filter( 1.2d0 * qi / qmax, a, nn )
|
||||
tab (iq, nb) = vqint * filter( qi / qmax, a, nn )
|
||||
!
|
||||
DO ir = 1, upf(nt)%kkbeta
|
||||
beta (ir,nb) = beta (ir,nb) + qi*qi * dq * tab0(iq,nb) * besr (ir) * rgrid(nt)%r(ir)
|
||||
betas(ir,nb) = betas(ir,nb) + qi*qi * dq * tab (iq,nb) * besr (ir) * rgrid(nt)%r(ir)
|
||||
ENDDO
|
||||
ENDDO
|
||||
!-
|
||||
ENDDO
|
||||
!
|
||||
CALL mp_sum( tab0, intra_bgrp_comm )
|
||||
CALL mp_sum( tab , intra_bgrp_comm )
|
||||
beta = beta * 8.d0/fpi ; CALL mp_sum( beta , intra_bgrp_comm )
|
||||
betas= betas* 8.d0/fpi ; CALL mp_sum( betas, intra_bgrp_comm )
|
||||
upf(nt)%beta(1:upf(nt)%kkbeta,1:upf(nt)%nbeta) = betas(1:upf(nt)%kkbeta,1:upf(nt)%nbeta)
|
||||
!-
|
||||
!- compute the beta normalization in q space to see how much of the original is lost
|
||||
DO nb = 1, upf(nt)%nbeta
|
||||
l = upf(nt)%lll (nb)
|
||||
power_q(nb) =0.d0
|
||||
DO iq = startq, lastq
|
||||
qi = (iq - 1) * dq
|
||||
power_q(nb) = power_q(nb) + qi*qi * dq * tab0(iq,nb) * tab0(iq,nb)
|
||||
ENDDO
|
||||
ENDDO
|
||||
power_q = power_q * 8.d0/fpi ; CALL mp_sum( power_q, intra_bgrp_comm )
|
||||
!
|
||||
IF (tprint) WRITE( stdout, '(5X,a,1p,10e12.4)'), "norm2 Q:", (power_q(nb), nb=1,upf(nt)%nbeta)
|
||||
IF (tprint) WRITE( stdout, '(5X,a,1p,10e12.4)'), "1-ratio:", (1.D0-power_q(nb)/power_r(nb), nb=1,upf(nt)%nbeta)
|
||||
!
|
||||
error_estimate= MAXVAL(1.d0-power_q(1:upf(nt)%nbeta)/power_r(1:upf(nt)%nbeta)) * filter(1.d0,a,nn)**2
|
||||
IF (error_estimate > eps ) CALL errore( 'init_us_b0','R and Q norms of beta too different',nt)
|
||||
WRITE ( stdout, '(5x,"Smoothing truncation error estimate in Q",1pe12.4," < eps")' ) error_estimate
|
||||
!
|
||||
!-Fix the core cutoff radius kkbeta so that no more than eps of the integral is lost.
|
||||
! Parseval's identity is used to monitor the real space completeness of the integral of the(squarred) betas
|
||||
! Parseval's identity is used to monitor the real space completeness of the integral of betas**2
|
||||
!
|
||||
! Compute the integral in Fourier space
|
||||
! Compute the integral in Fourier space of the smoothed projectors
|
||||
DO nb = 1, upf(nt)%nbeta
|
||||
l = upf(nt)%lll (nb)
|
||||
power_q(nb) =0.d0
|
||||
|
@ -151,20 +176,27 @@ SUBROUTINE init_us_b0
|
|||
ENDDO
|
||||
power_q = power_q * 8.d0/fpi ; CALL mp_sum( power_q, intra_bgrp_comm )
|
||||
!
|
||||
! Compute the integral in real space up to the current value of kkbeta and try to reduce it if close enough to power_q
|
||||
kktest = upf(nt)%kkbeta + 1 ; test = 0.d0 ! initialize mesh limit and test value
|
||||
! Compute norms in real space up to the current value of kkbeta and reduce it if close enough to power_q
|
||||
kktest = upf(nt)%kkbeta + 1 ; test = 0.d0 ; first =.TRUE. ! initialize mesh limit and test value
|
||||
DO WHILE ( test < eps )
|
||||
test_save = test ; kktest = kktest - 1
|
||||
safe_test = test ; kktest = kktest - 1
|
||||
DO nb = 1, upf(nt)%nbeta
|
||||
l = upf(nt)%lll (nb)
|
||||
aux(1:kktest) = betas(1:kktest,nb) * betas(1:kktest,nb)
|
||||
CALL simpson( kktest, aux, rgrid(nt)%rab, power_r(nb) )
|
||||
ENDDO
|
||||
IF ( first ) THEN
|
||||
first = .FALSE.
|
||||
test = MAXVAL( 1.d0 - power_r(1:upf(nt)%nbeta)/power_q(1:upf(nt)%nbeta) )
|
||||
if (test>eps) WRITE (stdout,'(5X,"WARNING: R and Q norms disagree by",6X,1pe12.4," > eps !")') test
|
||||
power_q = power_r
|
||||
END IF
|
||||
test = MAXVAL( 1.d0 - power_r(1:upf(nt)%nbeta)/power_q(1:upf(nt)%nbeta) )
|
||||
! write( stdout, * ) kktest, rgrid(nt)%r(kktest), test
|
||||
ENDDO
|
||||
if (kktest == upf(nt)%kkbeta ) CALL errore('init_us_b0',' initial kkbeta too small', nt)
|
||||
upf(nt)%kkbeta = kktest + 1 ! set the
|
||||
WRITE( stdout, * ) 'PSEUDO', nt, upf(nt)%kkbeta, rgrid(nt)%r(upf(nt)%kkbeta), test_save
|
||||
upf(nt)%kkbeta = kktest + 1 ! set the mesh limit to the last value fulfilling the test condition
|
||||
WRITE( stdout, '(5X,"kkbeta, r(kkbeta), error in R",i5,f7.3,1pe11.4," < eps")' ) &
|
||||
upf(nt)%kkbeta, rgrid(nt)%r(upf(nt)%kkbeta), safe_test
|
||||
!
|
||||
IF (tprint) THEN
|
||||
!-
|
||||
|
@ -209,6 +241,13 @@ SUBROUTINE init_us_b0
|
|||
ENDIF
|
||||
!
|
||||
ENDDO
|
||||
!- end of loop over pseudopotentials
|
||||
IF (tprint) THEN
|
||||
WRITE( stdout,'(/13X,"final rcut indices:",10i7)') upf(1:ntyp)%kkbeta
|
||||
WRITE( stdout,'(19X,"rcut (ityp) :",10f7.3)'), (rgrid(nt)%r(upf(nt)%kkbeta),nt=1,ntyp)
|
||||
END IF
|
||||
WRITE( stdout, '(/5X,"PSEUDOPOTENTIAL end Beta Smoothing: ==========================================")' )
|
||||
!
|
||||
IF (tprint) THEN
|
||||
filename(1:4) = 'ffqq' ! the filter in reciprocal space
|
||||
OPEN (4,FILE=filename, FORM='formatted', STATUS='unknown')
|
||||
|
@ -216,20 +255,42 @@ SUBROUTINE init_us_b0
|
|||
WRITE(4, *) '# fillter : a=',a,', nn=', nn
|
||||
DO iq = 1, nqx
|
||||
qi = (iq - 1) * dq
|
||||
WRITE(4,'(12f16.10)') qi, filter(1.1d0*qi/qmax, a, nn)
|
||||
WRITE(4,'(12f16.10)') qi, filter( qi/qmax, a, nn)
|
||||
ENDDO
|
||||
CLOSE (4)
|
||||
ENDIF
|
||||
IF (tprint) THEN
|
||||
WRITE( stdout, * ) " PSEUDOPOTENTIAL REPORT; updated values "
|
||||
WRITE( stdout, * ) ' NDM :', ndm, ' KKBETA :', upf(1:ntyp)%kkbeta
|
||||
WRITE( stdout, * ) ' r(ndm):', (rgrid(nt)%r(upf(nt)%kkbeta),nt=1,ntyp)
|
||||
ENDIF
|
||||
DO nt = 1, ntyp
|
||||
filename = 'ffr' ! the filter in real space
|
||||
WRITE(filename( 4:4 ),'(i1)') nt
|
||||
!- (spherically) fourier transform the filter in real space
|
||||
ffrr(:) = 0.d0
|
||||
DO iq = startq, lastq
|
||||
qi = (iq - 1) * dq
|
||||
CALL sph_bes (upf(nt)%kkbeta, rgrid(nt)%r, qi, 0, besr)
|
||||
DO ir = 1, upf(nt)%kkbeta
|
||||
ffrr (ir) = ffrr (ir) + qi*qi * dq * filter( qi / qmax, a, nn) * besr (ir)
|
||||
ENDDO
|
||||
END DO
|
||||
CALL mp_sum( ffrr, intra_bgrp_comm )
|
||||
ffrr(:) = 8.D0 / fpi * ffrr(:)
|
||||
aux(1:upf(nt)%kkbeta) = ffrr(1:upf(nt)%kkbeta) * rgrid(nt)%r(1:upf(nt)%kkbeta)**2
|
||||
CALL simpson( upf(nt)%kkbeta, aux, rgrid(nt)%rab, vqint )
|
||||
OPEN (4,FILE=filename, FORM='formatted', STATUS='unknown')
|
||||
WRITE(4, *) '# the filter in real space'
|
||||
WRITE(4, *) '# fillter : a=',a,', nn=', nn
|
||||
DO ir = 1, upf(nt)%kkbeta
|
||||
WRITE(4,'(12f16.10)') rgrid(nt)%r(ir), ffrr(ir)
|
||||
ENDDO
|
||||
CLOSE (4)
|
||||
!-
|
||||
END DO
|
||||
END IF
|
||||
!
|
||||
DEALLOCATE( power_r, power_q )
|
||||
DEALLOCATE( tab0, tab )
|
||||
DEALLOCATE( beta, betas )
|
||||
DEALLOCATE( besr, aux )
|
||||
DEALLOCATE( besr, aux, ffrr )
|
||||
!
|
||||
CALL stop_clock( 'init_us_b0' )
|
||||
!
|
||||
|
|
Loading…
Reference in New Issue