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

final (?) installment of the beta function smoothing. for each pseudo the chosen filter is the mildest that deliver the required accuracy.

This commit is contained in:
Stefano de Gironcoli 2020-04-18 16:40:49 +02:00
parent e257858304
commit f8b044582e
1 changed files with 175 additions and 117 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

292
PW/src/init_us_b0.f90
View File

@ -23,11 +23,12 @@ 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=14 ! 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.
! FILTER PARAMETERS: see REAL(DP) FUNCTION filter( x, a, n ) below for full definition and meaning.
INTEGER :: nf ! Smoothing parameter, order of the polynomial in the inverse gaussian approximant.
REAL(DP):: af ! Smoothing parameter, exponent of the gaussian decaying factor, to be chosen so that
! the filter has a flex at x=2/3 => af = 1.125d0 * ( 2 * nf + 1) .
!
LOGICAL, PARAMETER :: tprint=.FALSE. ! Whether the beta_l(r) and its relatives are printed or not.
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.
@ -45,11 +46,9 @@ SUBROUTINE init_us_b0
! q-point grid for interpolation
! interpolated value
INTEGER :: kktest
REAL(DP) :: test, safe_test, error_estimate
REAL(DP) :: test, safe_test, error_estimate, missing_norm
LOGICAL :: first
!
CHARACTER(LEN=4) :: filename
!
CALL start_clock( 'init_us_b0' )
!
! ... Initialization of variables
@ -80,9 +79,7 @@ SUBROUTINE init_us_b0
ALLOCATE( aux(ndm), besr(ndm), ffrr(ndm) )
!
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
WRITE(6,'(5X,a,1p,e12.4)') 'target error: eps =', eps
!
! ... fill the interpolation table tab
!
@ -93,20 +90,6 @@ SUBROUTINE init_us_b0
WRITE( stdout, '(/5X,a,i4)' ) 'Smoothing PSEUDO #', nt
IF ( upf(nt)%is_gth ) CYCLE
!
tab0(:,:) = 0.d0 ; tab(:,:) = 0.d0
beta(:,:) = 0.d0 ; betas(:,:) = 0.d0
!
IF (tprint) THEN
filename = 'br_' ! the radial beta_l(r) as defined in the pseudopotential
WRITE(filename( 4:4 ),'(i1)') nt
OPEN(4, FILE=filename, FORM='formatted', STATUS='unknown')
WRITE(4,*) '# the radial beta_l(r) as defined in the pseudopotential'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' kkbeta :',upf(nt)%kkbeta
DO ir = 1, upf(nt)%kkbeta
WRITE(4,'(12f16.10)') rgrid(nt)%r(ir), (upf(nt)%beta(ir,nb), nb=1,upf(nt)%nbeta)
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)
@ -117,8 +100,12 @@ SUBROUTINE init_us_b0
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)
!- write the radial beta_l(r) as defined in the pseudopotential
CALL write_beta_r( 'br_', nt, upf(nt)%beta, upf(nt)%kkbeta, upf(nt)%nbeta )
END IF
!-
!
!- compute the fourier transform of the beta and their (truncated) back f.t. in real space
tab0(:,:) = 0.d0 ; beta(:,:) = 0.d0
DO nb = 1, upf(nt)%nbeta
l = upf(nt)%lll (nb)
DO iq = startq, lastq
@ -130,20 +117,15 @@ SUBROUTINE init_us_b0
CALL simpson( upf(nt)%kkbeta, aux, rgrid(nt)%rab, vqint )
!
tab0(iq, nb) = vqint
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)
@ -154,18 +136,47 @@ SUBROUTINE init_us_b0
ENDDO
ENDDO
power_q = power_q * 8.d0/fpi ; CALL mp_sum( power_q, intra_bgrp_comm )
IF (tprint) THEN
WRITE( stdout, '(5X,a,1p,10e12.4)'), "norm2 Q:", (power_q(nb), nb=1,upf(nt)%nbeta)
WRITE( stdout, '(5X,a,1p,10e12.4)'), "1-ratio:", (1.D0-power_q(nb)/power_r(nb), nb=1,upf(nt)%nbeta)
END IF
!- fraction of beta normalizations that is lost due to Q space truncation
missing_norm = MAXVAL(1.d0-power_q(1:upf(nt)%nbeta)/power_r(1:upf(nt)%nbeta))
!
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)
!- determine the mildest filter that deliver the required accuracy in reciprocal space
nf = 20 ; af = 1.125D0 * ( 2 * nf + 1) ; error_estimate= missing_norm * filter(1.d0,af,nf)**2
do while ( error_estimate < eps .and. nf > 0)
nf = nf - 1 ; af = 1.125D0 * ( 2 * nf + 1) ; error_estimate = missing_norm * filter(1.d0,af,nf)**2
end do
nf = min(nf+1,20) ; af = 1.125D0 * ( 2 * nf + 1); error_estimate = missing_norm * filter(1.d0,af,nf)**2
WRITE ( stdout, '(5x,"Smoothing truncation error estimate in Q",1pe12.4," < eps")' ) error_estimate
IF (error_estimate > eps ) CALL errore( 'init_us_b0','R and Q norms of beta are too different',nt)
WRITE(6,'(5X,a,a,f6.2,a,i4,2(a,f11.8))') 'FILTER parameters :', &
' a=',af,', nn=',nf,', filter(1.0/3)=', filter(1.d0/3,af,nf), ', filter(1.0)=', filter(1.d0,af,nf)
!
!- compute the fourier transform of the smoothed beta and their (truncated) back f.t. in real space
tab(:,:) = 0.d0 ; betas(:,:) = 0.d0
DO nb = 1, upf(nt)%nbeta
l = upf(nt)%lll (nb)
DO iq = startq, lastq
qi = (iq - 1) * dq
CALL sph_bes (upf(nt)%kkbeta, rgrid(nt)%r, qi, l, besr)
!
tab (iq, nb) = tab0(iq, nb) * filter( qi / qmax, af, nf )
!
DO ir = 1, upf(nt)%kkbeta
betas(ir,nb) = betas(ir,nb) + qi*qi * dq * tab (iq,nb) * besr (ir) * rgrid(nt)%r(ir)
ENDDO
ENDDO
ENDDO
CALL mp_sum( tab , 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)
!
!-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 betas**2
!
! Compute the integral in Fourier space of the smoothed projectors
!- 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
@ -176,7 +187,7 @@ SUBROUTINE init_us_b0
ENDDO
power_q = power_q * 8.d0/fpi ; CALL mp_sum( power_q, intra_bgrp_comm )
!
! Compute norms in real space up to the current value of kkbeta and reduce it if close enough to power_q
!- 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 )
safe_test = test ; kktest = kktest - 1
@ -199,45 +210,18 @@ SUBROUTINE init_us_b0
upf(nt)%kkbeta, rgrid(nt)%r(upf(nt)%kkbeta), safe_test
!
IF (tprint) THEN
!-
filename(1:3) = 'bq_' ! the radial fourier transform of beta_l in reciprcal space up to qmax
OPEN( 4, FILE=filename, FORM='formatted', STATUS='unknown' )
WRITE(4,*) '# the radial fourier transform of beta_l in reciprcal space up to qmax'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' nqx :',nqx, dq*nqx
DO iq = 1, nqx
qi = (iq - 1) * dq
WRITE(4,'(12f16.10)') qi,(tab0(iq,nb), nb=1,upf(nt)%nbeta)
ENDDO
CLOSE(4)
!-
filename(1:3) = 'bqs' ! the smoothed radial fourier transform of beta_l in reciprcal space up to qmax
OPEN (4,FILE=filename,FORM='formatted', STATUS='unknown')
WRITE(4,*) '# the smoothed radial fourier transform of beta_l in reciprcal space up to qmax'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' nqx :',nqx, dq*nqx
do iq=1,nqx
qi = (iq - 1) * dq
WRITE(4,'(12f16.10)') qi,(tab(iq,nb), nb=1,upf(nt)%nbeta)
ENDDO
CLOSE(4)
!-
filename(1:3) = 'brq' ! the back radial fourier transform of beta_l in real space
OPEN (4,FILE=filename, FORM='formatted', STATUS='unknown')
WRITE(4,*) '# the back radial fourier transform of beta_l in real space'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' kkbeta :',upf(nt)%kkbeta
do ir =1,upf(nt)%kkbeta
WRITE(4,'(12f16.10)') rgrid(nt)%r(ir), (beta(ir,nb), nb=1,upf(nt)%nbeta)
ENDDO
CLOSE(4)
!-
filename(1:3) = 'brs' ! the back radial fourier transform of the smoothed beta_l in real space
OPEN(4,FILE=filename, FORM='formatted', STATUS='unknown')
WRITE(4,*) '# the back radial fourier transform of the smoothed beta_l in real space'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' kkbeta :',upf(nt)%kkbeta
do ir = 1, upf(nt)%kkbeta
WRITE(4,'(12f16.10)') rgrid(nt)%r(ir), (betas(ir,nb), nb=1,upf(nt)%nbeta)
ENDDO
CLOSE(4)
!-
!- write the radial fourier transform of beta_l in reciprcal space up to qmax
CALL write_beta_q( 'bq_', nt, tab0 )
!- write the smoothed radial fourier transform of beta_l in reciprcal space up to qmax
CALL write_beta_q( 'bqs', nt, tab )
!- write the back radial fourier transform of beta_l in real space
CALL write_beta_r( 'brq', nt, beta, ndm, nbetam )
!- write the back radial fourier transform of the smoothed beta_l in real space
CALL write_beta_r( 'brs', nt, betas, ndm, nbetam )
!- write the filter in reciprocal space
CALL write_filter_q( 'ffq', nt, af, nf )
!- write the filter in real space
CALL write_filter_r( 'ffr', nt, af, nf )
ENDIF
!
ENDDO
@ -248,45 +232,6 @@ SUBROUTINE init_us_b0
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')
WRITE(4, *) '# the filter in reciprocal space'
WRITE(4, *) '# fillter : a=',a,', nn=', nn
DO iq = 1, nqx
qi = (iq - 1) * dq
WRITE(4,'(12f16.10)') qi, filter( qi/qmax, a, nn)
ENDDO
CLOSE (4)
ENDIF
IF (tprint) THEN
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 )
@ -319,5 +264,118 @@ SUBROUTINE init_us_b0
!
END FUNCTION filter
!
! I/O routines used when tprint = .TRUE.
!
SUBROUTINE write_beta_r( filename_, nt, beta_, ndm, nbetam )
!
IMPLICIT NONE
!
CHARACTER(LEN=3), INTENT(IN) :: filename_
INTEGER, INTENT(IN) :: nt, ndm, nbetam
REAL(DP), INTENT(IN) :: beta_(ndm, nbetam)
!
CHARACTER(LEN=4) :: filename
INTEGER :: ir, nb
filename(1:3) = filename_; WRITE(filename( 4:4 ),'(i1)') nt
OPEN(4, FILE=filename, FORM='formatted', STATUS='unknown')
if (filename(1:3)=='br_') WRITE(4,*) '# the radial beta_l(r) as defined in the pseudopotential'
if (filename(1:3)=='brq') WRITE(4,*) '# the back radial fourier transform of beta_l in real space'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' kkbeta :',upf(nt)%kkbeta
if ( upf(nt)%nbeta > nbetam ) CALL errore('init_us_b0','wrong nbetam in write_beta_r', nbetam )
if ( upf(nt)%kkbeta > ndm ) CALL errore('init_us_b0','wrong ndm in write_beta_r', ndm )
DO ir = 1, upf(nt)%kkbeta
WRITE(4,'(12f16.10)') rgrid(nt)%r(ir), (beta_(ir,nb), nb=1,upf(nt)%nbeta)
ENDDO
CLOSE(4)
!
END SUBROUTINE write_beta_r
SUBROUTINE write_beta_q( filename_, nt, tab_ )
!
IMPLICIT NONE
!
CHARACTER(LEN=3), INTENT(IN) :: filename_
INTEGER, INTENT(IN) :: nt
REAL(DP), INTENT(IN) :: tab_(nqx, nbetam)
!
CHARACTER(LEN=4) :: filename
INTEGER :: iq, nb
REAL(DP) :: qi
!-
filename(1:3) = filename_; WRITE(filename( 4:4 ),'(i1)') nt
OPEN(4, FILE=filename, FORM='formatted', STATUS='unknown')
if (filename(1:3)=='bq_') WRITE(4,*) &
'# the radial fourier transform of beta_l in reciprcal space up to qmax'
if (filename(1:3)=='bqs') WRITE(4,*) &
'# the smoothed radial fourier transform of beta_l in reciprcal space up to qmax'
WRITE(4,*) '# nbeta :', upf(nt)%nbeta,' nqx :',nqx, dq*nqx
DO iq = 1, nqx
qi = (iq - 1) * dq
WRITE(4,'(12f16.10)') qi,(tab_(iq,nb), nb=1,upf(nt)%nbeta)
ENDDO
CLOSE(4)
!
END SUBROUTINE write_beta_q
!
SUBROUTINE write_filter_q( filename_, nt, af, nf )
!
IMPLICIT NONE
!
CHARACTER(LEN=3), INTENT(IN) :: filename_
INTEGER, INTENT(IN) :: nt, nf
REAL(DP), INTENT(IN) :: af
!
CHARACTER(LEN=4) :: filename
INTEGER :: iq
REAL(DP) :: qi
filename(1:3) = filename_; WRITE(filename( 4:4 ),'(i1)') nt
OPEN(4,FILE=filename, FORM='formatted', STATUS='unknown')
WRITE(4, *) '# the filter in reciprocal space'
WRITE(4, *) '# fillter : a=',af,', nn=', nf
DO iq = 1, nqx
qi = (iq - 1) * dq
WRITE(4,'(12f16.10)') qi, filter( qi/qmax, af, nf)
ENDDO
CLOSE(4)
!
END SUBROUTINE write_filter_q
SUBROUTINE write_filter_r( filename_, nt, af, nf )
!
IMPLICIT NONE
!
CHARACTER(LEN=3), INTENT(IN) :: filename_
INTEGER, INTENT(IN) :: nt, nf
REAL(DP), INTENT(IN) :: af
!
CHARACTER(LEN=4) :: filename
INTEGER :: iq
REAL(DP) :: qi
!- compute the (spherical) fourier transform of 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, af, nf) * besr (ir)
ENDDO
END DO
CALL mp_sum( ffrr, intra_bgrp_comm )
ffrr(:) = 8.D0 / fpi * ffrr(:)
!- write it on file
filename(1:3) = filename_; WRITE(filename( 4:4 ),'(i1)') nt
OPEN(4,FILE=filename, FORM='formatted', STATUS='unknown')
WRITE(4, *) '# the filter in real space'
WRITE(4, *) '# fillter : a=',af,', nn=', nf
DO ir = 1, upf(nt)%kkbeta
WRITE(4,'(12f16.10)') rgrid(nt)%r(ir), ffrr(ir)
ENDDO
CLOSE(4)
!
END SUBROUTINE write_filter_r
END SUBROUTINE init_us_b0