mirror of https://gitlab.com/QEF/q-e.git
579 lines
22 KiB
Fortran
579 lines
22 KiB
Fortran
|
|
! Copyright (C) 2001-2007 Quantum ESPRESSO group
|
|
! This file is distributed under the terms of the
|
|
! GNU General Public License. See the file `License'
|
|
! in the root directory of the present distribution,
|
|
! or http://www.gnu.org/copyleft/gpl.txt .
|
|
!
|
|
!
|
|
!-----------------------------------------------------------------------------------
|
|
SUBROUTINE init_us_0(ecutrho,intra_bgrp_comm)
|
|
!---------------------------------------------------------------------------------
|
|
!! This routine performs the following task: for each uspp or paw pseudopotential
|
|
!! the l-dependent aumentation charge \(\text{ q_nb_mb_l}\)(r), stored in
|
|
!! \(\text{qfuncl}\)(ir,nmb,l), is:
|
|
!
|
|
!! * transformed in reciprocal space by bessel transform up to qmax = sqrt(ecutrho);
|
|
!! * smoothed by multiplying with a filter function \(\textrm{filter}\)(q/qmax,a,nn);
|
|
!! * brought back in real space,
|
|
!
|
|
!! where it overwrites the original array. The filter function is:
|
|
!! \[ \text{filter}(x,a,\text{nn}) = e^{-\text{axx}} \sum_{k=0,\text{nn}}
|
|
!! \frac{\text{axx}^k}{k!}\ . \]
|
|
!
|
|
USE upf_kinds, ONLY: DP
|
|
USE upf_io, ONLY: stdout
|
|
USE upf_const, ONLY: fpi, sqrt2, eps8, eps6
|
|
USE atom, ONLY: rgrid
|
|
USE uspp_param, ONLY: nsp, upf, lmaxq, nbetam
|
|
USE mp, ONLY: mp_sum
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
REAL(DP), INTENT(IN) :: ecutrho
|
|
INTEGER, INTENT(IN) :: intra_bgrp_comm
|
|
!
|
|
! ... local variables
|
|
!
|
|
! sdg
|
|
! LOGICAL, PARAMETER :: tprint=.true. ! whether the q_l(r) and its relatives are printed or not
|
|
LOGICAL, PARAMETER :: tprint=.FALSE. ! whether the q_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.0_DP ! smoothing parameter, exponent of the gaussian decaying factor
|
|
REAL(DP), PARAMETER :: dq = 0.01_dp !nterpolation table
|
|
! a=0.0 ; nn=0.0 would be no smoothing.
|
|
!
|
|
INTEGER :: nt, ih, jh, nb, mb, ijv, l, m, ir, ir0, iq, is, startq, lastq, ilast, ndm, ia
|
|
! various counters
|
|
INTEGER :: nqxq
|
|
REAL(DP), ALLOCATABLE :: qrad_q(:,:,:), qrad_r(:,:,:), qrad_rs(:,:,:), ffrr(:)
|
|
REAL(DP), ALLOCATABLE :: power_0(:,:), power_r(:,:), power_q(:,:), power_rs(:,:), power_qs(:,:)
|
|
REAL(DP), ALLOCATABLE :: aux(:), aux1(:)
|
|
! various work space
|
|
REAL(DP) :: q, qi
|
|
! the modulus of g for each shell
|
|
! q-point grid for interpolation
|
|
REAL(DP), ALLOCATABLE :: ylmk0(:)
|
|
! the spherical harmonics
|
|
INTEGER :: lnb, lmb
|
|
REAL(DP) :: qmax, rcut, drcut
|
|
REAL(DP) :: target_ratio, ratio, ratio_s, fac
|
|
!
|
|
CHARACTER(LEN=6) :: filename
|
|
!
|
|
CALL start_clock( 'init_us_0' )
|
|
!
|
|
! ... Initialization of the variables
|
|
!
|
|
drcut = ABS(LOG(eps8))/SQRT(ecutrho)
|
|
!
|
|
DO nt = 1, nsp
|
|
!
|
|
rcut = rgrid(nt)%r(upf(nt)%kkbeta)
|
|
WRITE (stdout,*) 'RCUT:', rcut, drcut, rcut+drcut
|
|
DO ir = upf(nt)%kkbeta, upf(nt)%mesh
|
|
IF ( rgrid(nt)%r(ir) < rcut+drcut ) upf(nt)%kkbeta=ir
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
ndm = MAXVAL( upf(:)%kkbeta )
|
|
nqxq = INT( SQRT(ecutrho) / dq + 4 )
|
|
!
|
|
IF (tprint) THEN
|
|
WRITE (stdout,*) " PSEUDOPOTENTIAL REPORT "
|
|
WRITE (stdout,*) ' NDM :', ndm, ' ', upf(1:nsp)%kkbeta
|
|
WRITE (stdout,*) ' LMAXQ :', lmaxq, ' NBETAM :', nbetam
|
|
ENDIF
|
|
!
|
|
ALLOCATE( aux(ndm), aux1(ndm) )
|
|
ALLOCATE( qrad_q(nqxq, nbetam*(nbetam+1)/2,lmaxq) )
|
|
ALLOCATE( qrad_r(ndm , nbetam*(nbetam+1)/2,lmaxq) )
|
|
ALLOCATE( qrad_rs(ndm, nbetam*(nbetam+1)/2,lmaxq) )
|
|
ALLOCATE( ffrr(ndm) )
|
|
ALLOCATE( power_0(nbetam*(nbetam+1)/2 ,0:lmaxq) )
|
|
ALLOCATE( power_r(nbetam*(nbetam+1)/2 ,0:lmaxq), power_q(nbetam*(nbetam+1)/2 ,0:lmaxq) )
|
|
ALLOCATE( power_rs(nbetam*(nbetam+1)/2,0:lmaxq), power_qs(nbetam*(nbetam+1)/2,0:lmaxq) )
|
|
ALLOCATE( ylmk0(lmaxq*lmaxq) )
|
|
!
|
|
! ... the following prevents an out-of-bound error: upf(nt)%nqlc=2*lmax+1
|
|
! but in some versions of the PP files lmax is not set to the maximum
|
|
! l of the beta functions but includes the l of the local potential
|
|
!
|
|
DO nt = 1, nsp
|
|
!
|
|
upf(nt)%nqlc = MIN( upf(nt)%nqlc, lmaxq )
|
|
IF (upf(nt)%nqlc < 0) upf(nt)%nqlc = 0
|
|
!
|
|
ENDDO
|
|
!
|
|
! ... Here for the US types we compute the Fourier transform of the
|
|
! Q functions.
|
|
!
|
|
CALL divide( intra_bgrp_comm, nqxq, startq, lastq )
|
|
!
|
|
qmax = SQRT(ecutrho)
|
|
WRITE (stdout, *) ' qmax : sqrt(ecutrho) =', SQRT(ecutrho), dq*nqxq
|
|
WRITE (stdout,'(a,f6.2,a,i4,4(a,f11.8))') 'FILTER : a=',a,', nn=',nn, &
|
|
', filter(1.1d0)=', filter(1.1d0,a,nn), &
|
|
', filter(1.0d0)=', filter(1.0d0,a,nn), &
|
|
', filter(0.9d0)=', filter(0.9d0,a,nn), &
|
|
', filter(0.8d0)=', filter(0.8d0,a,nn)
|
|
!
|
|
DO nt = 1, nsp
|
|
WRITE (stdout,*) ' NT = ', nt
|
|
!
|
|
IF ( upf(nt)%tvanp ) THEN
|
|
!-
|
|
IF (tprint) THEN
|
|
filename = 'qr_' ! the radial q_l(r) as defined in the pseudopotential
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb - 1) / 2 + nb
|
|
lnb = upf(nt)%lll(nb) ; lmb = upf(nt)%lll(mb)
|
|
WRITE (filename(4:4),'(i1)') nb
|
|
WRITE (filename(5:5),'(i1)') mb
|
|
WRITE (filename(6:6),'(i1)') nt
|
|
OPEN (4,file=filename,form='formatted', status='unknown')
|
|
WRITE (4, *) '# the radial q_l(r) as defined in the pseudopotential'
|
|
WRITE (4, *) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, &
|
|
' kkbeta :',upf(nt)%kkbeta
|
|
!
|
|
DO ir =1,upf(nt)%kkbeta
|
|
write (4,'(12f16.10)') rgrid(nt)%r(ir), (upf(nt)%qfuncl(ir,ijv,l), l=0, lnb+lmb)
|
|
ENDDO
|
|
!
|
|
CLOSE (4)
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDIF
|
|
!-
|
|
qrad_q(:,:,:) = 0.0_DP ; qrad_r(:,:,:) = 0.0_DP
|
|
power_0(:,:) = 0.0_DP ; power_r(:,:) = 0.0_DP ; power_q(:,:) = 0.0_DP
|
|
!
|
|
DO l = 0, upf(nt)%nqlc-1
|
|
!
|
|
! 1) first of all compute the integrated power spectrum of the Qs in real space.
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb - 1) / 2 + nb
|
|
aux1(1) = 0.0_DP
|
|
!
|
|
ir0 = 1
|
|
IF (rgrid(nt)%r(ir0) < eps8) ir0 = 2
|
|
!
|
|
DO ir = ir0, upf(nt)%kkbeta
|
|
aux1(ir) = (upf(nt)%qfuncl(ir,ijv,l)/rgrid(nt)%r(ir))**2
|
|
ENDDO
|
|
!
|
|
CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, power_0(ijv,l+1) )
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! 2) compute the fourier transform of the Qs and their integrated power spectum in
|
|
! reciprocal space - FIXME: use routine init_tab_qrad
|
|
!
|
|
DO iq = startq, lastq
|
|
!
|
|
q = (iq - 1) * dq
|
|
!
|
|
! ... here we compute the spherical bessel function for the given |q|
|
|
!
|
|
CALL sph_bes( upf(nt)%kkbeta, rgrid(nt)%r, q, l, aux )
|
|
!
|
|
! .. and here we integrate with all the Q functions
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb - 1) / 2 + nb
|
|
lnb = upf(nt)%lll(nb)
|
|
lmb = upf(nt)%lll(mb)
|
|
!
|
|
IF ( (l >= ABS(lnb-lmb)) .AND. (l <= lnb+lmb) .AND. (MOD(l+lnb+lmb,2)==0) ) THEN
|
|
!
|
|
DO ir = 1, upf(nt)%kkbeta
|
|
aux1 (ir) = aux (ir) * upf(nt)%qfuncl(ir,ijv,l)
|
|
ENDDO
|
|
!
|
|
! ... computes the Fourier transform of q_nb_mb_l(r) up to qmax
|
|
!
|
|
CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, qrad_q(iq,ijv,l+1) )
|
|
!
|
|
! ... and update the integrated power spectrum in reciprocal space
|
|
!
|
|
power_q(ijv,l+1) = power_q(ijv,l+1) + q*q * dq * qrad_q(iq,ijv,l+1)**2
|
|
!
|
|
ENDIF
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! 3) back-fourier transform of the Qs to real space.
|
|
!
|
|
DO ir =1, upf(nt)%kkbeta
|
|
!
|
|
! ... q_nb_mb_l(r) from the back fourier transform up to qmax of q_nb_mb_l(q)
|
|
!
|
|
qrad_r(ir,1:nbetam*(nbetam+1)/2,l+1) = qrad_r(ir,1:nbetam*(nbetam+1)/2,l+1) + &
|
|
q*q * dq * aux(ir) * rgrid(nt)%r(ir)**2 * &
|
|
qrad_q(iq,1:nbetam*(nbetam+1)/2,l+1)
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
CALL mp_sum( qrad_q , intra_bgrp_comm )
|
|
CALL mp_sum( power_q, intra_bgrp_comm ) ; power_q (:,:) = power_q (:,:) * 8.0_DP/fpi
|
|
CALL mp_sum( qrad_r , intra_bgrp_comm ) ; qrad_r(:,:,:) = qrad_r(:,:,:) * 8.0_DP/fpi
|
|
!
|
|
! 4) compute integrated power spectrum of the Qs in real space (completeness check).
|
|
!
|
|
DO l = 0, upf(nt)%nqlc-1
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb - 1)/2 + nb
|
|
aux1(1) = 0.0_DP
|
|
!
|
|
ir0 = 1
|
|
IF (rgrid(nt)%r(ir0) < eps8) ir0 = 2
|
|
!
|
|
DO ir = ir0, upf(nt)%kkbeta
|
|
aux1(ir) = (qrad_r(ir,ijv,l+1)/rgrid(nt)%r(ir))**2
|
|
ENDDO
|
|
!
|
|
CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, power_r(ijv,l+1) )
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
! ... computes the ratio of power_0 and power_q to see how complete is the reciprocal
|
|
! space expansion.
|
|
!
|
|
power_0(:,0) = 0.0_DP
|
|
power_q(:,0) = 0.0_DP
|
|
power_r(:,0) = 0.0_DP
|
|
target_ratio = 1.0_DP
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb-1) / 2 + nb
|
|
lnb = upf(nt)%lll(nb)
|
|
lmb = upf(nt)%lll(mb)
|
|
!
|
|
DO l = 0, lnb+lmb
|
|
power_0(ijv,0) = power_0(ijv,0) + power_0(ijv,l+1)
|
|
power_q(ijv,0) = power_q(ijv,0) + power_q(ijv,l+1)
|
|
power_r(ijv,0) = power_r(ijv,0) + power_r(ijv,l+1)
|
|
ENDDO
|
|
!
|
|
!write (stdout, *) ' nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb
|
|
!write (stdout,'(a,12f16.10)') 'power_0 ',(power_0 (ijv,l+1), l=0,lnb+lmb)
|
|
!write (stdout,'(a,12f16.10)') 'power_r ',(power_r (ijv,l+1), l=0,lnb+lmb)
|
|
!write (stdout,'(a,12f16.10)') 'power_q ',(power_q (ijv,l+1), l=0,lnb+lmb)
|
|
!write (stdout,*) 'ratio ',1.d0-(power_r (ijv,0)/power_q (ijv,0)), &
|
|
! 1.d0-(power_r (ijv,0)/power_0 (ijv,0))
|
|
IF (power_0(ijv,0)>eps8) target_ratio = MIN(target_ratio, power_r(ijv,0)/power_0(ijv,0))
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
WRITE (stdout,*) ' TARGET Qs SPILLOVER : 1.d0-target_ratio, target_ratio ', &
|
|
1.d0-target_ratio, target_ratio
|
|
!
|
|
fac = 1.2_DP
|
|
!
|
|
99 CONTINUE
|
|
!
|
|
! ... smooth the Fourier transform of the Qs and compute its integrated power spectrum
|
|
!
|
|
power_qs(:,:) = 0.0_DP
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb - 1) / 2 + nb
|
|
lnb = upf(nt)%lll(nb)
|
|
lmb = upf(nt)%lll(mb)
|
|
!
|
|
DO l = 0, upf(nt)%nqlc-1
|
|
!
|
|
IF ( (l >= ABS(lnb-lmb)) .AND. (l <= lnb+lmb) .AND. (MOD(l+lnb+lmb,2)==0) ) THEN
|
|
DO iq = startq, lastq
|
|
q = (iq - 1) * dq
|
|
power_qs(ijv,l+1) = power_qs(ijv,l+1) + q*q * dq * &
|
|
(qrad_q(iq,ijv,l+1)*filter(fac*q/qmax,a,nn))**2
|
|
ENDDO
|
|
ENDIF
|
|
!
|
|
power_qs(ijv,0) = power_qs(ijv,0) + power_qs(ijv,l+1)
|
|
!
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
power_qs(:,:) = power_qs(:,:) * 8.0_DP/fpi
|
|
!
|
|
CALL mp_sum( power_qs, intra_bgrp_comm )
|
|
!
|
|
ratio = 1.0_DP ; ratio_s = 1.0_DP
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
ijv = mb * (mb - 1) / 2 + nb
|
|
IF (power_q(ijv,0)>eps8) ratio = MIN(ratio, power_qs(ijv,0)/power_q(ijv,0))
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
! WRITE (stdout,*) ' filter factor and power_qs/power_q ratio:', fac, ratio
|
|
!
|
|
! ... with a given fac the smoothed Qs qre built
|
|
!
|
|
qrad_rs(:,:,:) = 0.0_DP ; ffrr(:) = 0.0_DP
|
|
!
|
|
DO l = 0, upf(nt)%nqlc-1
|
|
!
|
|
DO iq = startq, lastq
|
|
!
|
|
q = (iq - 1) * dq
|
|
!
|
|
! ... here we compute the spherical bessel function for the given |q| ...
|
|
!
|
|
CALL sph_bes( upf(nt)%kkbeta, rgrid(nt)%r, q, l, aux )
|
|
!
|
|
! ... and here we integrate with all the Q functions
|
|
!
|
|
! 5) back-fourier transform of the smoothed Qs to real space.
|
|
!
|
|
DO ir = 1, upf(nt)%kkbeta
|
|
!
|
|
! ... q_nb_mb_l(r) from the back fourier transform up to qmax of q_nb_mb_l(q)
|
|
!
|
|
qrad_rs(ir,1:nbetam*(nbetam+1)/2,l+1) = qrad_rs(ir,1:nbetam*(nbetam+1)/2,l+1) &
|
|
+ aux(ir) * q*q * dq * rgrid(nt)%r(ir)**2 &
|
|
* qrad_q(iq,1:nbetam*(nbetam+1)/2,l+1) * filter(fac*q/qmax,a,nn)
|
|
!
|
|
! ... build the filter function in real space from the back fourier transform up
|
|
! to qmax
|
|
!
|
|
IF (l==0) ffrr(ir) = ffrr(ir) + q*q * dq * aux(ir) * rgrid(nt)%r(ir)**2 &
|
|
* filter(fac*q/qmax,a,nn)
|
|
!
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
CALL mp_sum( qrad_rs, intra_bgrp_comm ) ; qrad_rs(:,:,:) = qrad_rs(:,:,:) * 8.0_DP/fpi
|
|
CALL mp_sum( ffrr, intra_bgrp_comm ) ; ffrr(:) = ffrr(:) * 8.0_DP/fpi
|
|
!
|
|
! 6) compute intergrated power spectrum of the Qs in real space (completeness check).
|
|
!
|
|
DO l = 0, upf(nt)%nqlc-1
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb-1)/2 + nb
|
|
aux1(1) = 0.0_DP
|
|
!
|
|
ir0 = 1
|
|
IF (rgrid(nt)%r(ir0) < eps8) ir0 = 2
|
|
!
|
|
DO ir = ir0, upf(nt)%kkbeta
|
|
aux1(ir) = (qrad_rs(ir,ijv,l+1)/rgrid(nt)%r(ir))**2
|
|
ENDDO
|
|
!
|
|
CALL simpson( upf(nt)%kkbeta, aux1, rgrid(nt)%rab, power_rs(ijv,l+1) )
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
!
|
|
power_rs(:,0) = 0.0_DP
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb-1)/2 + nb
|
|
lnb = upf(nt)%lll(nb)
|
|
lmb = upf(nt)%lll(mb)
|
|
!
|
|
DO l = 0, lnb+lmb
|
|
power_rs(ijv,0) = power_rs(ijv,0) + power_rs(ijv,l+1)
|
|
ENDDO
|
|
!
|
|
!write (stdout, *) ' nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb
|
|
!write (stdout,'(a,12f16.10)') 'power_rs',(power_rs(ijv,l+1), l=0,lnb+lmb)
|
|
!write (stdout,'(a,12f16.10)') 'power_qs',(power_qs(ijv,l+1), l=0,lnb+lmb)
|
|
!write (stdout,*) 'ratio ',1.d0-(power_rs(ijv,0)/power_qs(ijv,0)), &
|
|
! 1.d0-(power_qs(ijv,0)/power_q (ijv,0))
|
|
IF (power_qs(ijv,0) > eps8) ratio_s = MIN(ratio_s, power_rs(ijv,0)/power_qs(ijv,0))
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
WRITE (stdout,*) ' filter factor, 1.d0-power_rs/power_qs and power_qs/power_q ratio :', &
|
|
fac, 1.d0-ratio_s, ratio
|
|
fac = fac - 0.05_DP
|
|
!sdg
|
|
IF (ratio < target_ratio .AND. (1.d0-ratio_s) < eps8) GOTO 99
|
|
!
|
|
! IF (ratio < target_ratio .AND. (1.d0-ratio_s) < eps6) GOTO 99
|
|
!
|
|
fac = fac + 0.05_DP ! reset the last successful value of fac
|
|
!
|
|
!- save the smoothed real space qs in qfuncl
|
|
!
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb-1)/2 + nb
|
|
lnb = upf(nt)%lll(nb)
|
|
lmb = upf(nt)%lll(mb)
|
|
!
|
|
DO l = 0, upf(nt)%nqlc-1
|
|
!
|
|
IF ( (l>=ABS(lnb-lmb)) .AND. (l<=lnb+lmb) .AND. (MOD(l+lnb+lmb,2)==0) ) THEN
|
|
!
|
|
upf(nt)%qfuncl(1:upf(nt)%kkbeta,ijv,l) = qrad_rs(1:upf(nt)%kkbeta,ijv,l+1)
|
|
!
|
|
ENDIF
|
|
!
|
|
ENDDO
|
|
!
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDIF
|
|
!
|
|
IF (tprint) THEN
|
|
!-
|
|
filename = 'ffff'
|
|
OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown')
|
|
WRITE (4, *) '# filter function : a=',a,', nn=',nn,', fac=', fac
|
|
WRITE (4, *) '# nqxq :', nqxq,' dq :',dq, ' qmax :',qmax
|
|
DO iq = 1, nqxq
|
|
q = (iq-1)*dq
|
|
WRITE (4,'(2f16.10)') q, filter( fac*q/qmax, a, nn )
|
|
ENDDO
|
|
CLOSE (4)
|
|
!-
|
|
filename = 'ffrr'
|
|
OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown')
|
|
WRITE (4, *) '# filter function : a=',a,', nn=',nn,', fac=', fac
|
|
WRITE (4, *) '# kkbeta :', upf(nt)%kkbeta
|
|
DO ir = 1, upf(nt)%kkbeta
|
|
write (4,'(2f16.10)') rgrid(nt)%r(ir), ffrr(ir)
|
|
ENDDO
|
|
CLOSE (4)
|
|
!-
|
|
DO nb = 1, upf(nt)%nbeta
|
|
DO mb = nb, upf(nt)%nbeta
|
|
!
|
|
ijv = mb * (mb - 1) / 2 + nb
|
|
lnb = upf(nt)%lll(nb) ; lmb = upf(nt)%lll(mb)
|
|
WRITE (filename(4:4),'(i1)') nb
|
|
WRITE (filename(5:5),'(i1)') mb
|
|
WRITE (filename(6:6),'(i1)') nt
|
|
!-
|
|
filename(1:3) = 'qq_' ! the radial fourier transform of q_l in reciprocal space
|
|
OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown')
|
|
WRITE (4,*) '# the radial fourier transform of q_l in reciprocal space'
|
|
WRITE (4,*) '# nb :', nb, lnb,' mb :', mb, lmb,' lmax :', lnb+lmb, ' nqxq :', nqxq
|
|
DO iq=1,nqxq
|
|
q = (iq-1)*dq
|
|
WRITE (4,'(12f16.10)') q, (qrad_q(iq,ijv,l+1), l=0,lnb+lmb )
|
|
ENDDO
|
|
CLOSE (4)
|
|
!-
|
|
filename(1:3) = 'qqs' ! the smoothed radial fourier transform of q_l in reciprocal space
|
|
OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown')
|
|
WRITE (4,*) '# the smoothed radial fourier transform of q_l in reciprocal space'
|
|
WRITE (4,*) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, ' nqxq :',nqxq
|
|
DO iq = 1, nqxq
|
|
q = (iq-1)*dq
|
|
WRITE (4,'(12f16.10)') q,(qrad_q(iq,ijv,l+1)*filter(fac*q/qmax,a,nn), l=0,lnb+lmb )
|
|
ENDDO
|
|
CLOSE (4)
|
|
!-
|
|
filename(1:3) = 'qrq' ! the radial q_l(r) as obtained back-transforming the q_l(q)
|
|
OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown')
|
|
WRITE (4,*) '# the radial q_l(r) as obtained back-transforming the q_l(q)'
|
|
WRITE (4,*) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, ' kkbeta :',upf(nt)%kkbeta
|
|
DO ir = 1, upf(nt)%kkbeta
|
|
WRITE (4,'(12f16.10)') rgrid(nt)%r(ir),(qrad_r(ir,ijv,l+1), l=0,lnb+lmb )
|
|
ENDDO
|
|
CLOSE (4)
|
|
!-
|
|
filename(1:3) = 'qrs' ! the radial q_l(r) as obtained back-transforming the smoothed q_l(q)
|
|
OPEN (4, FILE=filename, FORM='formatted', STATUS='unknown')
|
|
WRITE (4,*) '# the radial q_l(r) as obtained back-transforming the smoothed q_l(q)'
|
|
WRITE (4,*) '# nb :', nb,lnb,' mb :',mb,lmb,' lmax :',lnb+lmb, ' kkbeta :',upf(nt)%kkbeta
|
|
DO ir = 1, upf(nt)%kkbeta
|
|
WRITE (4,'(12f16.10)') rgrid(nt)%r(ir),(qrad_rs(ir,ijv,l+1), l=0,lnb+lmb )
|
|
ENDDO
|
|
CLOSE (4)
|
|
!-
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ENDIF
|
|
! nsp
|
|
ENDDO
|
|
!
|
|
DEALLOCATE( ylmk0 )
|
|
DEALLOCATE( power_0, power_r, power_q, power_rs, power_qs )
|
|
DEALLOCATE( qrad_q, qrad_r, qrad_rs, ffrr )
|
|
DEALLOCATE( aux1, aux )
|
|
!
|
|
CALL stop_clock( 'init_us_0' )
|
|
!
|
|
RETURN
|
|
!
|
|
!
|
|
!
|
|
CONTAINS
|
|
!
|
|
REAL(DP) FUNCTION filter( x, a, n )
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
REAL(DP), INTENT(IN) :: x, a
|
|
INTEGER, INTENT(IN) :: n
|
|
REAL(DP) :: axx, ff
|
|
INTEGER :: k
|
|
!
|
|
axx = a * x * x
|
|
ff = 1.0_DP
|
|
!
|
|
DO k = n, 1, -1
|
|
ff = (1._DP + axx/DBLE(k)*ff)
|
|
ENDDO
|
|
!
|
|
filter = ff * EXP(-axx)
|
|
!
|
|
RETURN
|
|
!
|
|
END FUNCTION filter
|
|
!
|
|
!
|
|
END SUBROUTINE init_us_0
|