quantum-espresso/upflib/init_us_0.f90

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