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

Fixed stupid error in "read_file" that was causing crashes in all 

phonon and postprocessing codes (Andrea please check the spin-orbit case)
Redundant variables in uspp_param (mostly) removed


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4313 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2007-10-06 08:23:39 +00:00
parent b8099605db
commit 6d97829a14
24 changed files with 123 additions and 363 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

15
CPV/cprsub.f90
View File

@ -27,7 +27,7 @@ subroutine formf( tfirst, eself )
use ions_base, ONLY : rcmax, zv, nsp, na use ions_base, ONLY : rcmax, zv, nsp, na
use local_pseudo, ONLY : vps, rhops, dvps, drhops use local_pseudo, ONLY : vps, rhops, dvps, drhops
use atom, ONLY : rgrid, numeric use atom, ONLY : rgrid, numeric
use uspp_param, ONLY : vloc_at, oldvan use uspp_param, ONLY : upf, oldvan
use pseudo_base, ONLY : compute_rhops, formfn, formfa, compute_eself use pseudo_base, ONLY : compute_rhops, formfn, formfa, compute_eself
use pseudopotential, ONLY : tpstab, vps_sp, dvps_sp use pseudopotential, ONLY : tpstab, vps_sp, dvps_sp
use cp_interfaces, ONLY : build_pstab use cp_interfaces, ONLY : build_pstab
@ -83,9 +83,9 @@ subroutine formf( tfirst, eself )
if ( numeric(is) ) then if ( numeric(is) ) then
call formfn( vps(:,is), dvps(:,is), rgrid(is)%r, rgrid(is)%rab, vloc_at(:,is), & call formfn( vps(:,is), dvps(:,is), rgrid(is)%r, rgrid(is)%rab, &
zv(is), rcmax(is), g, omega, tpiba2, rgrid(is)%mesh, & upf(is)%vloc(1), zv(is), rcmax(is), g, omega, tpiba2,&
ngs, oldvan(is), tpre ) rgrid(is)%mesh, ngs, oldvan(is), tpre )
else else
@ -332,8 +332,7 @@ subroutine nlinit
use constants, ONLY : pi, fpi use constants, ONLY : pi, fpi
use ions_base, ONLY : na, nsp use ions_base, ONLY : na, nsp
use uspp, ONLY : aainit, beta, qq, dvan, nhtol, nhtolm, indv use uspp, ONLY : aainit, beta, qq, dvan, nhtol, nhtolm, indv
use uspp_param, ONLY : kkbeta, qqq, nqlc, betar, lmaxq, dion,& use uspp_param, ONLY : upf, lmaxq, nbetam, lmaxkb, nhm, nh
nbeta, nbetam, lmaxkb, lll, nhm, nh, tvanp
use atom, ONLY : rgrid, nlcc, numeric use atom, ONLY : rgrid, nlcc, numeric
use qradb_mod, ONLY : qradb use qradb_mod, ONLY : qradb
use qgb_mod, ONLY : qgb use qgb_mod, ONLY : qgb
@ -410,8 +409,8 @@ subroutine nlinit
! !
WRITE( stdout,*) WRITE( stdout,*)
WRITE( stdout,'(20x,a)') ' dion ' WRITE( stdout,'(20x,a)') ' dion '
do iv = 1, nbeta(is) do iv = 1, upf(is)%nbeta
WRITE( stdout,'(8f9.4)') ( fac*dion(iv,jv,is), jv = 1, nbeta(is) ) WRITE( stdout,'(8f9.4)') ( fac*upf(is)%dion(iv,jv), jv = 1, upf(is)%nbeta )
end do end do
! !
end do end do

2
CPV/modules.f90
View File

@ -53,8 +53,6 @@ end module bhs
! !
! lqmax: maximum angular momentum of Q (Vanderbilt augmentation charges) ! lqmax: maximum angular momentum of Q (Vanderbilt augmentation charges)
! nqfx : maximum number of coefficients in Q smoothing
! nbrx : maximum number of distinct radial beta functions
! ndmx: maximum number of points in the radial grid ! ndmx: maximum number of points in the radial grid
! !

7
CPV/pseudopot_sub.f90
View File

@ -389,7 +389,7 @@
USE pseudo_base, ONLY : nlin_base USE pseudo_base, ONLY : nlin_base
USE pseudo_base, ONLY : nlin_stress_base USE pseudo_base, ONLY : nlin_stress_base
USE reciprocal_vectors, ONLY : g, gstart USE reciprocal_vectors, ONLY : g, gstart
USE uspp_param, ONLY : upf, nbeta, nbetam USE uspp_param, ONLY : upf, nbetam
USE read_pseudo_module_fpmd, ONLY : nspnl USE read_pseudo_module_fpmd, ONLY : nspnl
USE cp_interfaces, ONLY : compute_xgtab, chkpstab USE cp_interfaces, ONLY : compute_xgtab, chkpstab
USE pseudopotential, ONLY : wnl_sp, wnla_sp, xgtab USE pseudopotential, ONLY : wnl_sp, wnla_sp, xgtab
@ -561,7 +561,7 @@
USE cell_base, ONLY : tpiba USE cell_base, ONLY : tpiba
USE pseudo_base, ONLY : nlin_stress_base USE pseudo_base, ONLY : nlin_stress_base
USE splines, ONLY : spline USE splines, ONLY : spline
USE uspp_param, ONLY : upf, nbeta USE uspp_param, ONLY : upf
USE read_pseudo_module_fpmd, ONLY : nspnl USE read_pseudo_module_fpmd, ONLY : nspnl
USE reciprocal_vectors, ONLY : g, gstart USE reciprocal_vectors, ONLY : g, gstart
USE pseudopotential, ONLY : wnla_sp, tpstab USE pseudopotential, ONLY : wnla_sp, tpstab
@ -1624,8 +1624,7 @@
do l = lmin, lmax do l = lmin, lmax
do ir = 1, upf(is)%kkbeta do ir = 1, upf(is)%kkbeta
if ( rgrid(is)%r(ir) >= upf(is)%rinner(l) ) then if ( rgrid(is)%r(ir) >= upf(is)%rinner(l) ) then
! qrl(ir,ijv,l)=upf(is)%qfunc(ir,ijv) TEMP qrl(ir,ijv,l)=upf(is)%qfunc(ir,ijv)
qrl(ir,ijv,l)=upf(is)%qfunc(ir,iv,jv)
else else
qrl(ir,ijv,l)=upf(is)%qfcoef(1,l,iv,jv) qrl(ir,ijv,l)=upf(is)%qfcoef(1,l,iv,jv)
do i = 2, upf(is)%nqf do i = 2, upf(is)%nqf

26
CPV/qqberry.f90
View File

@ -17,7 +17,7 @@ subroutine qqberry2( gqq,gqqm, ipol)
use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, & use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, &
nr1bx, nr2bx, nr3bx, nnrb => nnrbx nr1bx, nr2bx, nr3bx, nnrb => nnrbx
use uspp_param, only: lmaxq, nqlc, kkbeta, nbeta, nbetam, nh, nhm, oldvan use uspp_param, only: upf, lmaxq, nbetam, nh, nhm, oldvan
use uspp, only: indv, lpx, lpl, ap,nhtolm use uspp, only: indv, lpx, lpl, ap,nhtolm
use atom, only: rgrid use atom, only: rgrid
use core use core
@ -50,7 +50,7 @@ subroutine qqberry2( gqq,gqqm, ipol)
integer :: ivs, jvs, ivl, jvl, lp, ijv integer :: ivs, jvs, ivl, jvl, lp, ijv
real(8), allocatable:: ylm(:,:) real(8), allocatable:: ylm(:,:)
ndm = MAXVAL (kkbeta(1:nsp)) ndm = MAXVAL (upf(1:nsp)%kkbeta)
allocate( fint( ndm), jl(ndm)) allocate( fint( ndm), jl(ndm))
allocate( qradb2(nbetam,nbetam,lmaxq,nsp)) allocate( qradb2(nbetam,nbetam,lmaxq,nsp))
allocate( ylm(ngw, lmaxq*lmaxq)) allocate( ylm(ngw, lmaxq*lmaxq))
@ -86,28 +86,30 @@ subroutine qqberry2( gqq,gqqm, ipol)
do is=1,nvb do is=1,nvb
c=fpi !/omegab c=fpi !/omegab
! !
ALLOCATE ( qrl(kkbeta(is), nbeta(is)*(nbeta(is)+1)/2, nqlc(is)) ) ALLOCATE ( qrl( upf(is)%kkbeta, upf(is)%nbeta*(upf(is)%nbeta+1)/2, &
upf(is)%nqlc ) )
! !
call fill_qrl ( is, qrl ) call fill_qrl ( is, qrl )
! now the radial part ! now the radial part
do l=1,nqlc(is) do l=1,upf(is)%nqlc
xg= gmes !only orthorombic cells xg= gmes !only orthorombic cells
!!!call bess(xg,l,kkbeta(is),r(1,is),jl) !!!call bess(xg,l,upf(is)%kkbeta,rgrid(is)%r,jl)
call sph_bes ( kkbeta(is), rgrid(is)%r, xg, l-1, jl ) call sph_bes ( upf(is)%kkbeta, rgrid(is)%r, xg, l-1, jl )
do iv= 1,nbeta(is) do iv= 1,upf(is)%nbeta
do jv=iv,nbeta(is) do jv=iv,upf(is)%nbeta
ijv = (jv-1)*jv/2 + iv ijv = (jv-1)*jv/2 + iv
! !
! note qrl(r)=r^2*q(r) ! note qrl(r)=r^2*q(r)
! !
do ir=1,kkbeta(is) do ir=1,upf(is)%kkbeta
fint(ir)=qrl(ir,ijv,l)*jl(ir) fint(ir)=qrl(ir,ijv,l)*jl(ir)
end do end do
if (oldvan(is)) then if (oldvan(is)) then
call herman_skillman_int & call herman_skillman_int ( upf(is)%kkbeta,fint,rgrid(is)%rab,&
& (kkbeta(is),fint,rgrid(is)%rab,qradb2(iv,jv,l,is)) qradb2(iv,jv,l,is) )
else else
call simpson (kkbeta(is),fint,rgrid(is)%rab,qradb2(iv,jv,l,is)) call simpson ( upf(is)%kkbeta,fint,rgrid(is)%rab,&
qradb2(iv,jv,l,is) )
endif endif
qradb2(iv,jv,l,is)= c*qradb2(iv,jv,l,is) qradb2(iv,jv,l,is)= c*qradb2(iv,jv,l,is)
if ( iv /= jv ) qradb2(jv,iv,l,is)= qradb2(iv,jv,l,is) if ( iv /= jv ) qradb2(jv,iv,l,is)= qradb2(iv,jv,l,is)

178
CPV/read_pseudo.f90
View File

@ -319,7 +319,7 @@ END FUNCTION calculate_dx
IF( program_name == 'FPMD' ) THEN IF( program_name == 'FPMD' ) THEN
CALL errore(' readpp ', ' file format not supported ', 1 ) CALL errore(' readpp ', ' file format not supported ', 1 )
ELSE ELSE
call readbhs(is,pseudounit) CALL errore(' readpp ', ' file format no longer supported ', 2 )
END IF END IF
END IF END IF
@ -533,179 +533,3 @@ END FUNCTION calculate_dx
END MODULE read_pseudo_module_fpmd END MODULE read_pseudo_module_fpmd
!=----------------------------------------------------------------------------=! !=----------------------------------------------------------------------------=!
! !
!
!
!
!---------------------------------------------------------------------
subroutine readbhs( is, iunps )
!---------------------------------------------------------------------
!
use atom, only: rgrid, nlcc, rho_atc, numeric
use uspp_param, only: zp, betar, dion, vloc_at, lll, nbeta, kkbeta
use bhs, only: rcl, rc2, bl, al, wrc1, lloc, wrc2, rc1
use funct, only: set_dft_from_name, dft_is_hybrid
use io_global, only: stdout
!
implicit none
!
integer is, iunps
!
integer meshp, ir, ib, il, i, j, jj
real(8), allocatable:: fint(:), vnl(:)
real(8) rdum, alpha, z, zval, cmesh, cmeshp, exfact
character(len=20) :: dft_name
!
! nlcc is unfortunately not read from file
!
numeric(is) = .false.
nlcc(is)=.false.
read(iunps,*) z,zp(is),nbeta(is),lloc(is),exfact
if (zp(is) < 1 .or. zp(is) > 100 ) then
call errore('readbhs','wrong potential read',15)
endif
call dftname_cp (nint(exfact), dft_name)
call set_dft_from_name( dft_name )
IF ( dft_is_hybrid() ) &
CALL errore( 'readbhs', 'HYBRID XC not implemented in CPV', 1 )
!
if(lloc(is).eq.2)then
lll(1,is)=0
lll(2,is)=1
else if(lloc(is).ne.2) then
call errore('readbhs','kb-ization for lloc=2 only',10)
endif
!
! see eqs. (2.21) and (2.22) of bhs, prb 26, 4199 (1982).
!
! wrc1 =c_core(1)
! wrc2 =c_core(2)
! rc1 =alpha_core(1)
! rc2 =alpha_core(2)
! al(i) =a(i) i=1,3
! bl(i) =a(i+3) i=1,3
! rcl(i)=alpha(i) i=1,3
!
! ------------------------------------------------------------------
! pp parameters are read from file iunps
! bhs 's coefficients have been turned into lengths
! ------------------------------------------------------------------
read(iunps,*) wrc1(is),rc1(is),wrc2(is),rc2(is)
rc1(is)=1.0d0/sqrt(rc1(is))
rc2(is)=1.0d0/sqrt(rc2(is))
do il=1,3
do ib=1,3
read(iunps,*) rcl(ib,is,il),al(ib,is,il),bl(ib,is,il)
rcl(ib,is,il)=1.0d0/sqrt(rcl(ib,is,il))
end do
end do
!
! ------------------------------------------------------------------
! wavefunctions are read from file iunps
! ------------------------------------------------------------------
do il=1,nbeta(is)
read(iunps,*) rgrid(is)%mesh,cmesh
!
! kkbeta is for compatibility with Vanderbilt PP
!
kkbeta(is)=rgrid(is)%mesh
do j=1,rgrid(is)%mesh
read(iunps,*) jj,rgrid(is)%r(j),betar(j,il,is)
end do
end do
!
! ------------------------------------------------------------------
! core charge is read from unit 15
! ------------------------------------------------------------------
!
if(nlcc(is)) then
read(15,*) meshp,cmeshp
if ( meshp.ne.rgrid(is)%mesh .or. cmeshp.ne.cmesh ) then
call errore('readbhs','core charge mesh mismatch',is)
endif
do ir=1,rgrid(is)%mesh
read(15,*) rdum, rho_atc(ir,is)
end do
endif
!
! rab(i) is the derivative of the radial mesh
!
do ir=1,rgrid(is)%mesh
rgrid(is)%rab(ir)=rgrid(is)%r(ir) * log(cmesh)
end do
!
! ------------------------------------------------------------------
! local potential
! ------------------------------------------------------------------
lloc(is)=lloc(is)+1
!
! NB: the following is NOT the local potential: the -ze^2/r term is missing
!
do ir=1,rgrid(is)%mesh
vloc_at(ir,is)=0.d0
do i=1,3
vloc_at(ir,is) = vloc_at(ir,is) &
& +(al(i,is,lloc(is))+bl(i,is,lloc(is))*rgrid(is)%r(ir)**2) &
& *exp(-(rgrid(is)%r(ir)/rcl(i,is,lloc(is)))**2)
end do
end do
!
! ------------------------------------------------------------------
! nonlocal potentials: kleinman-bylander form
! (1) definition of betar (2) calculation of dion
! ------------------------------------------------------------------
allocate(fint(rgrid(is)%mesh), vnl(rgrid(is)%mesh))
do il=1,nbeta(is)
do ir=1,rgrid(is)%mesh
vnl(ir)=0.d0
do i=1,3
vnl(ir) = vnl(ir) + (al(i,is,il)+bl(i,is,il)*rgrid(is)%r(ir)**2)&
& * exp(-(rgrid(is)%r(ir)/rcl(i,is,il))**2)
end do
vnl(ir) = vnl(ir) - vloc_at(ir,is)
fint(ir)= betar(ir,il,is)**2*vnl(ir)
betar(ir,il,is)=vnl(ir)*betar(ir,il,is)
end do
call simpson_cp90(rgrid(is)%mesh,fint,rgrid(is)%rab,dion(il,il,is))
dion(il,il,is) = 1.0d0/dion(il,il,is)
end do
deallocate(vnl, fint)
!
! ------------------------------------------------------------------
! output: pp info
! ------------------------------------------------------------------
WRITE( stdout,3000) z,zp(is)
3000 format(2x,'bhs pp for z=',f3.0,2x,'zv=',f3.0)
WRITE( stdout,'(2x,a20)') dft_name
WRITE( stdout,3002) lloc(is)-1
3002 format(2x,' local angular momentum: l=',i3)
WRITE( stdout,3005) nbeta(is)
3005 format(2x,'number of nl ang. mom. nbeta=',i3)
do il=1,nbeta(is)
WRITE( stdout,3010) lll(il,is)
3010 format(2x,'nonlocal angular momentum: l=',i3)
end do
WRITE( stdout,3030)
3030 format(2x,'pseudopotential parameters:')
WRITE( stdout,3035) wrc1(is),1.0d0/rc1(is)**2
3035 format(2x,'core:',2x,'c1_c=',f7.4,' alpha1_c=',f7.4)
WRITE( stdout,3036) wrc2(is),1.0d0/rc2(is)**2
3036 format(2x,' ',2x,'c2_c=',f7.4,' alpha2_c=',f7.4)
WRITE( stdout,3038)
3038 format(2x,'other table parameters:')
do il=1,3
WRITE( stdout,3040) il-1
3040 format(2x,'l=',i3)
do i =1,3
alpha=1.0d0/rcl(i,is,il)**2
WRITE( stdout,3050) i,alpha,i,al(i,is,il),i+3,bl(i,is,il)
end do
end do
3050 format(2x,'alpha',i1,'=',f6.2,' a',i1,'=',f16.7, &
& ' a',i1,'=',f16.7)
WRITE( stdout,*)
!
return
end subroutine readbhs

6
CPV/stress.f90
View File

@ -60,7 +60,7 @@
USE cell_base, ONLY: tpiba2, omega USE cell_base, ONLY: tpiba2, omega
USE reciprocal_vectors, ONLY: gstart, gzero, g, gx USE reciprocal_vectors, ONLY: gstart, gzero, g, gx
USE gvecw, ONLY: ngw USE gvecw, ONLY: ngw
USE uspp_param, only: nh, lmaxkb, nbeta, nhm USE uspp_param, only: nbetam, nh, lmaxkb, nhm
USE uspp, only: nhtol, nhtolm, indv, nkb USE uspp, only: nhtol, nhtolm, indv, nkb
USE electrons_base, only: nupdwn, iupdwn, nspin USE electrons_base, only: nupdwn, iupdwn, nspin
USE cdvan, only: dbec USE cdvan, only: dbec
@ -150,8 +150,8 @@
! ... initialize array wnla ! ... initialize array wnla
ALLOCATE( wsg ( nhm, nsp ) ) ALLOCATE( wsg ( nhm, nsp ) )
ALLOCATE( wnl ( ngw, MAXVAL( nbeta( 1:nsp ) ), nsp, 1 ) ) ALLOCATE( wnl ( ngw, nbetam, nsp, 1 ) )
ALLOCATE( wnla( ngw, MAXVAL( nbeta( 1:nsp ) ), nsp ) ) ALLOCATE( wnla( ngw, nbetam, nsp ) )
ALLOCATE( fnlb( nsanl, MAXVAL( nh( 1:nsp ) ), MAXVAL( nupdwn ), nspin ) ) ALLOCATE( fnlb( nsanl, MAXVAL( nh( 1:nsp ) ), MAXVAL( nupdwn ), nspin ) )
! !
fac = sqrt( omega ) / ( 2.0d0 * 4.0d0 * pi ) fac = sqrt( omega ) / ( 2.0d0 * 4.0d0 * pi )

2
Modules/grid_paw_variables.f90
View File

@ -8,7 +8,7 @@ module grid_paw_variables
! NO rinner > 0 ! NO rinner > 0
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
USE parameters, ONLY : lqmax, nbrx, npsx, nqfx USE parameters, ONLY : lqmax, nbrx, npsx
USE radial_grids, ONLY: ndmx USE radial_grids, ONLY: ndmx
! !
implicit none implicit none

4
Modules/parameters.f90
View File

@ -20,14 +20,12 @@ MODULE parameters
npsx = ntypx, &! max number of different PPs (obsolete) npsx = ntypx, &! max number of different PPs (obsolete)
npk = 40000, &! max number of k-points npk = 40000, &! max number of k-points
lmaxx = 3, &! max non local angular momentum (l=0 to lmaxx) lmaxx = 3, &! max non local angular momentum (l=0 to lmaxx)
cp_lmax = lmaxx+1,&! max number of channels !TMP FOR PAW: REMOVE
nchix = 6, & ! max number of atomic wavefunctions per atom nchix = 6, & ! max number of atomic wavefunctions per atom
nwfsx = 14 ! max number of beta functions per atom nwfsx = 14 ! max number of beta functions per atom
! !
INTEGER, PARAMETER :: & INTEGER, PARAMETER :: &
nbrx = 14, &! max number of beta functions nbrx = 14, &! max number of beta functions
lqmax= 2*lmaxx+1, &! max number of angular momenta of Q lqmax= 2*lmaxx+1 ! max number of angular momenta of Q
nqfx = 8 ! max number of coefficients in Q smoothing
! !
! ... More parameter for the CP codes ! ... More parameter for the CP codes

55
Modules/pseudo_types.f90
View File

@ -12,8 +12,7 @@
! together with their allocation/deallocation routines ! together with their allocation/deallocation routines
USE kinds, ONLY: DP USE kinds, ONLY: DP
USE parameters, ONLY: cp_lmax, lmaxx use radial_grids, ONLY: radial_grid_type
use radial_grids, ONLY: ndmx, radial_grid_type
IMPLICIT NONE IMPLICIT NONE
SAVE SAVE
@ -81,44 +80,46 @@ END TYPE paw_t
REAL(DP) :: rmax ! the maximum radius of the mesh REAL(DP) :: rmax ! the maximum radius of the mesh
REAL(DP) :: zmesh ! the nuclear charge used for mesh REAL(DP) :: zmesh ! the nuclear charge used for mesh
REAL(DP) :: dx ! the deltax of the linear mesh REAL(DP) :: dx ! the deltax of the linear mesh
INTEGER, POINTER :: nn(:) ! nn(nwfc) INTEGER, POINTER :: nn(:) ! nn(nwfc) quantum number of wfc
REAL(DP), POINTER :: rcut(:) ! cut-off radius(nbeta) REAL(DP), POINTER :: rcut(:) ! cut-off radius(nbeta)
REAL(DP), POINTER :: rcutus(:)! cut-off ultrasoft radius (nbeta) REAL(DP), POINTER :: rcutus(:)! ultrasoft cut-off radius (nbeta)
REAL(DP), POINTER :: epseu(:) ! energy (nwfc) REAL(DP), POINTER :: epseu(:) ! energy (nwfc)
REAL(DP), POINTER :: jchi(:) ! jchi(nwfc) REAL(DP), POINTER :: jchi(:) ! jchi(nwfc) j=l+1/2 or l-1/2 of wfc
REAL(DP), POINTER :: jjj(:) ! jjj(nbeta) REAL(DP), POINTER :: jjj(:) ! jjj(nbeta) j=l+1/2 or l-1/2 of beta
INTEGER :: nv ! UPF file version number INTEGER :: nv ! UPF file version number
INTEGER :: lmax ! maximum angular momentum component INTEGER :: lmax ! maximum l component in beta
INTEGER :: mesh ! number of point in the radial mesh INTEGER :: mesh ! number of point in the radial mesh
INTEGER :: nwfc ! number of wavefunctions INTEGER :: nwfc ! number of atomic wavefunctions
INTEGER :: nbeta ! number of projectors INTEGER :: nbeta ! number of projectors
INTEGER :: kkbeta ! kkbeta=max(kbeta(:)) INTEGER :: kkbeta ! kkbeta=max(kbeta(:))
! kbeta<=mesh is the number of grid points for each beta function ! kbeta<=mesh is the number of grid points for each beta function
! beta(r,nb) = 0 for r > r(kbeta(nb)) ! beta(r,nb) = 0 for r > r(kbeta(nb))
! kkbeta<=mesh is the largest of such number so that for all beta ! kkbeta<=mesh is the largest of such number so that for all beta
! beta(r,nb) = 0 for r > r(kkbeta) ! beta(r,nb) = 0 for r > r(kkbeta)
CHARACTER(LEN=2), POINTER :: els(:) ! els(nwfc) CHARACTER(LEN=2), POINTER :: els(:) ! els(nwfc) label of wfc
CHARACTER(LEN=2), POINTER :: els_beta(:) ! els(nbeta) CHARACTER(LEN=2), POINTER :: els_beta(:) ! els(nbeta) label of beta
INTEGER, POINTER :: lchi(:) ! lchi(nwfc) INTEGER, POINTER :: lchi(:) ! lchi(nwfc) value of l for wavefcts
REAL(DP), POINTER :: oc(:) ! oc(nwfc) REAL(DP), POINTER :: oc(:) ! oc(nwfc) occupancies for wavefcts
REAL(DP), POINTER :: r(:) ! r(mesh) REAL(DP), POINTER :: r(:) ! r(mesh) radial grid
REAL(DP), POINTER :: rab(:) ! rab(mesh) REAL(DP), POINTER :: rab(:) ! rab(mesh) dr(x)/dx (x=linear grid)
REAL(DP), POINTER :: rho_atc(:) ! rho_atc(mesh) REAL(DP), POINTER :: rho_atc(:) ! rho_atc(mesh) atomic core charge
REAL(DP), POINTER :: vloc(:) ! vloc(mesh) REAL(DP), POINTER :: vloc(:) ! vloc(mesh) local atomic potential
INTEGER, POINTER :: lll(:) ! lll(nbeta) INTEGER, POINTER :: lll(:) ! lll(nbeta) l of each projector
INTEGER, POINTER :: kbeta(:) ! kbeta(nbeta) INTEGER, POINTER :: kbeta(:) ! kbeta(nbeta) see above kkbeta
REAL(DP), POINTER :: beta(:,:) ! beta(mesh,nbeta) REAL(DP), POINTER :: beta(:,:) ! beta(mesh,nbeta) projectors
INTEGER :: nd INTEGER :: nd
REAL(DP), POINTER :: dion(:,:) ! dion(nbeta,nbeta) REAL(DP), POINTER :: dion(:,:) ! dion(nbeta,nbeta) atomic D_{mu,nu}
INTEGER :: nqf INTEGER :: nqf ! number of Q coefficients
INTEGER :: nqlc INTEGER :: nqlc ! number of angular momenta in Q
REAL(DP), POINTER :: rinner(:) ! rinner(0:2*lmax) REAL(DP), POINTER :: rinner(:) ! rinner(0:2*lmax) r_L
REAL(DP), POINTER :: qqq(:,:) ! qqq(nbeta,nbeta) REAL(DP), POINTER :: qqq(:,:) ! qqq(nbeta,nbeta) q_{mu,nu}
REAL(DP), POINTER :: qfunc(:,:,:) ! qfunc(mesh,nbeta,nbeta) REAL(DP), POINTER :: qfunc(:,:) ! qfunc(mesh,nbeta*(nbeta+1)/2)
! Q_{mu,nu}(|r|) function for |r|> r_L
REAL(DP), POINTER :: qfcoef(:,:,:,:) ! qfcoef(nqf,0:2*lmax,nbeta,nbeta) REAL(DP), POINTER :: qfcoef(:,:,:,:) ! qfcoef(nqf,0:2*lmax,nbeta,nbeta)
REAL(DP), POINTER :: chi(:,:) ! chi(mesh,nwfc) ! coefficients for Q for |r|<r_L
REAL(DP), POINTER :: rho_at(:) ! rho_at(mesh) REAL(DP), POINTER :: chi(:,:) ! chi(mesh,nwfc) atomic wavefcts
REAL(DP), POINTER :: rho_at(:) ! rho_at(mesh) atomic charge
LOGICAL :: has_paw ! Whether PAW data is included LOGICAL :: has_paw ! Whether PAW data is included
REAL(DP) :: paw_data_format ! The version of the format REAL(DP) :: paw_data_format ! The version of the format

5
Modules/read_paw.f90
View File

@ -459,12 +459,11 @@ subroutine set_pseudo_paw (is, pawset)
end do end do
! triangularize matrix of qfunc's ! triangularize matrix of qfunc's
allocate ( upf(is)%qfunc(1:pawset%grid%mesh,pawset%nwfc,pawset%nwfc) ) allocate ( upf(is)%qfunc(1:pawset%grid%mesh,pawset%nwfc*(pawset%nwfc+1)/2) )
do nb = 1, pawset%nwfc do nb = 1, pawset%nwfc
do mb = nb, pawset%nwfc do mb = nb, pawset%nwfc
ijv = mb * (mb-1) / 2 + nb ijv = mb * (mb-1) / 2 + nb
!!! qfunc (1:pawset%grid%mesh, ijv, is) = pawset%augfun(1:pawset%grid%mesh,nb,mb,0) upf(is)%qfunc (1:pawset%grid%mesh,ijv) = &
upf(is)%qfunc (1:pawset%grid%mesh, nb,mb) = &
pawset%augfun(1:pawset%grid%mesh,nb,mb,0) pawset%augfun(1:pawset%grid%mesh,nb,mb,0)
enddo enddo
enddo enddo

16
Modules/read_upf.f90
View File

@ -333,7 +333,7 @@ subroutine read_pseudo_nl (upf, iunps)
integer :: iunps integer :: iunps
TYPE (pseudo_upf), INTENT(INOUT) :: upf TYPE (pseudo_upf), INTENT(INOUT) :: upf
! !
integer :: nb, mb, n, ir, ios, idum, ldum, icon, lp, i, ikk integer :: nb, mb, ijv, n, ir, ios, idum, ldum, icon, lp, i, ikk
! counters ! counters
character (len=75) :: dummy character (len=75) :: dummy
! !
@ -346,7 +346,7 @@ subroutine read_pseudo_nl (upf, iunps)
ALLOCATE( upf%dion( 1, 1 ) ) ALLOCATE( upf%dion( 1, 1 ) )
ALLOCATE( upf%rinner( 1 ) ) ALLOCATE( upf%rinner( 1 ) )
ALLOCATE( upf%qqq ( 1, 1 ) ) ALLOCATE( upf%qqq ( 1, 1 ) )
ALLOCATE( upf%qfunc ( 0:upf%mesh, 1, 1 ) ) ALLOCATE( upf%qfunc ( 0:upf%mesh, 1 ) )
ALLOCATE( upf%qfcoef( 1, 1, 1, 1 ) ) ALLOCATE( upf%qfcoef( 1, 1, 1, 1 ) )
ALLOCATE( upf%rcut( 1 ) ) ALLOCATE( upf%rcut( 1 ) )
ALLOCATE( upf%rcutus( 1 ) ) ALLOCATE( upf%rcutus( 1 ) )
@ -399,7 +399,7 @@ subroutine read_pseudo_nl (upf, iunps)
upf%nqlc = 2 * upf%lmax + 1 upf%nqlc = 2 * upf%lmax + 1
ALLOCATE( upf%rinner( upf%nqlc ) ) ALLOCATE( upf%rinner( upf%nqlc ) )
ALLOCATE( upf%qqq ( upf%nbeta, upf%nbeta ) ) ALLOCATE( upf%qqq ( upf%nbeta, upf%nbeta ) )
ALLOCATE( upf%qfunc ( 0:upf%mesh, upf%nbeta, upf%nbeta ) ) ALLOCATE( upf%qfunc ( 0:upf%mesh, upf%nbeta*(upf%nbeta+1)/2 ) )
ALLOCATE( upf%qfcoef( MAX( upf%nqf,1 ), upf%nqlc, upf%nbeta, upf%nbeta ) ) ALLOCATE( upf%qfcoef( MAX( upf%nqf,1 ), upf%nqlc, upf%nbeta, upf%nbeta ) )
upf%rinner = 0.0_DP upf%rinner = 0.0_DP
upf%qqq = 0.0_DP upf%qqq = 0.0_DP
@ -422,11 +422,9 @@ subroutine read_pseudo_nl (upf, iunps)
read (iunps,*,err=104,end=104) upf%qqq(nb,mb), dummy read (iunps,*,err=104,end=104) upf%qqq(nb,mb), dummy
! "Q_int" ! "Q_int"
upf%qqq(mb,nb) = upf%qqq(nb,mb) upf%qqq(mb,nb) = upf%qqq(nb,mb)
! ijv is the combined (nb,mb) index
read (iunps, *, err=105, end=105) (upf%qfunc(n,nb,mb), n=1,upf%mesh) ijv = mb * (mb-1) / 2 + nb
do n = 0, upf%mesh read (iunps, *, err=105, end=105) (upf%qfunc(n,ijv), n=1,upf%mesh)
upf%qfunc(n,mb,nb) = upf%qfunc(n,nb,mb)
enddo
if ( upf%nqf > 0 ) then if ( upf%nqf > 0 ) then
call scan_begin (iunps, "QFCOEF", .false.) call scan_begin (iunps, "QFCOEF", .false.)
@ -448,7 +446,7 @@ subroutine read_pseudo_nl (upf, iunps)
upf%nqlc = 2 * upf%lmax + 1 upf%nqlc = 2 * upf%lmax + 1
ALLOCATE( upf%rinner( upf%nqlc ) ) ALLOCATE( upf%rinner( upf%nqlc ) )
ALLOCATE( upf%qqq ( upf%nbeta, upf%nbeta ) ) ALLOCATE( upf%qqq ( upf%nbeta, upf%nbeta ) )
ALLOCATE( upf%qfunc ( 0:upf%mesh, upf%nbeta, upf%nbeta ) ) ALLOCATE( upf%qfunc ( 0:upf%mesh, upf%nbeta*(upf%nbeta+1)/2 ) )
ALLOCATE( upf%qfcoef( upf%nqf, upf%nqlc, upf%nbeta, upf%nbeta ) ) ALLOCATE( upf%qfcoef( upf%nqf, upf%nqlc, upf%nbeta, upf%nbeta ) )
upf%rinner = 0.0_DP upf%rinner = 0.0_DP
upf%qqq = 0.0_DP upf%qqq = 0.0_DP

48
Modules/read_uspp.f90
View File

@ -13,7 +13,7 @@ MODULE read_uspp_module
! Vanderbilt's code and Andrea's RRKJ3 format ! Vanderbilt's code and Andrea's RRKJ3 format
! !
USE kinds, ONLY: DP USE kinds, ONLY: DP
USE parameters, ONLY: nchix, lmaxx, nbrx, nsx, lqmax, nqfx USE parameters, ONLY: nchix, lmaxx, nsx, lqmax
USE io_global, ONLY: stdout, meta_ionode USE io_global, ONLY: stdout, meta_ionode
USE funct, ONLY: set_dft_from_name, dft_is_hybrid, dft_is_meta, & USE funct, ONLY: set_dft_from_name, dft_is_hybrid, dft_is_meta, &
set_dft_from_indices set_dft_from_indices
@ -101,13 +101,16 @@ CONTAINS
& exfact, &! index of the exchange and correlation used & exfact, &! index of the exchange and correlation used
& etotpseu, &! total pseudopotential energy & etotpseu, &! total pseudopotential energy
& ee(nchix), &! the energy of the valence states & ee(nchix), &! the energy of the valence states
& rc(nchix), &! the cut-off radii of the pseudopotential
& eloc, &! energy of the local potential & eloc, &! energy of the local potential
& dummy, &! dummy real variable & dummy, &! dummy real variable
& rc(nchix), &! the cut-off radii of the pseudopotential
& eee(nbrx), &! energies of the beta function
& ddd(nbrx,nbrx),&! the screened D_{\mu,\nu} parameters
& rinner1, &! rinner if only one is present & rinner1, &! rinner if only one is present
& rcloc ! the cut-off radius of the local potential & rcloc ! the cut-off radius of the local potential
real(DP), allocatable:: &
& eee(:), &! energies of the beta function
& ddd(:,:) ! the screened D_{\mu,\nu} parameters
integer, allocatable :: &
& iptype(:) ! more recent parameters
integer & integer &
& idmy(3), &! contains the date of creation of the pseudo & idmy(3), &! contains the date of creation of the pseudo
& ifpcor, &! for core correction, 0 otherwise & ifpcor, &! for core correction, 0 otherwise
@ -115,9 +118,8 @@ CONTAINS
& i, &! dummy counter & i, &! dummy counter
& nnlz(nchix), &! The nlm values of the valence states & nnlz(nchix), &! The nlm values of the valence states
& keyps, &! the type of pseudopotential. Only US allowed & keyps, &! the type of pseudopotential. Only US allowed
& irel, &! it says if the pseudopotential is relativistic & irel, &! says if the pseudopotential is relativistic
& ifqopt, &! level of Q optimization & ifqopt, &! level of Q optimization
& iptype(nbrx), &! more recent parameters
& npf, &! as above & npf, &! as above
& nang, &! number of angular momenta in pseudopotentials & nang, &! number of angular momenta in pseudopotentials
& lloc, &! angular momentum of the local part of PPs & lloc, &! angular momentum of the local part of PPs
@ -213,7 +215,7 @@ CONTAINS
if ( lloc == -1 ) lloc = nang+1 if ( lloc == -1 ) lloc = nang+1
if ( lloc > nang+1 .or. lloc < 0 ) & if ( lloc > nang+1 .or. lloc < 0 ) &
call errore( 'readvan', 'wrong lloc read', is ) call errore( 'readvan', 'wrong lloc read', is )
if ( upf%nqf > nqfx .or. upf%nqf < 0 ) & if ( upf%nqf < 0 ) &
call errore(' readvan', 'Wrong nqf read', upf%nqf) call errore(' readvan', 'Wrong nqf read', upf%nqf)
if ( ifqopt < 0 ) & if ( ifqopt < 0 ) &
call errore( 'readvan', 'wrong ifqopt read', is ) call errore( 'readvan', 'wrong ifqopt read', is )
@ -272,8 +274,8 @@ CONTAINS
ALLOCATE ( upf%kbeta(upf%nbeta) ) ALLOCATE ( upf%kbeta(upf%nbeta) )
upf%kbeta(:) = upf%kkbeta upf%kbeta(:) = upf%kkbeta
! !
if( upf%nbeta > nbrx .or. upf%nbeta <0 ) & if( upf%nbeta < 0 ) &
call errore( 'readvan','nbeta wrong or too large', is ) call errore( 'readvan','nbeta wrong', is )
if( upf%kkbeta > upf%mesh .or. upf%kkbeta < 0 ) & if( upf%kkbeta > upf%mesh .or. upf%kkbeta < 0 ) &
call errore( 'readvan','kkbeta wrong or too large', is ) call errore( 'readvan','kkbeta wrong or too large', is )
! !
@ -282,8 +284,9 @@ CONTAINS
ALLOCATE ( upf%lll(upf%nbeta) ) ALLOCATE ( upf%lll(upf%nbeta) )
ALLOCATE ( upf%beta(upf%mesh,upf%nbeta) ) ALLOCATE ( upf%beta(upf%mesh,upf%nbeta) )
ALLOCATE ( upf%dion(upf%nbeta,upf%nbeta), upf%qqq(upf%nbeta,upf%nbeta) ) ALLOCATE ( upf%dion(upf%nbeta,upf%nbeta), upf%qqq(upf%nbeta,upf%nbeta) )
ALLOCATE ( upf%qfunc(upf%mesh,upf%nbeta,upf%nbeta) ) ALLOCATE ( upf%qfunc(upf%mesh,upf%nbeta*(upf%nbeta+1)/2) )
ALLOCATE ( upf%qfcoef(upf%nqf, upf%nqlc, upf%nbeta, upf%nbeta) ) ALLOCATE ( upf%qfcoef(upf%nqf, upf%nqlc, upf%nbeta, upf%nbeta) )
ALLOCATE ( eee(upf%nbeta), ddd(upf%nbeta,upf%nbeta) )
do iv=1,upf%nbeta do iv=1,upf%nbeta
read( iunps, '(i5)',err=100, iostat=ios ) upf%lll(iv) read( iunps, '(i5)',err=100, iostat=ios ) upf%lll(iv)
read( iunps, '(1p4e19.11)',err=100, iostat=ios ) & read( iunps, '(1p4e19.11)',err=100, iostat=ios ) &
@ -297,15 +300,14 @@ CONTAINS
! !
! the symmetric matric Q_{nb,mb} is stored in packed form ! the symmetric matric Q_{nb,mb} is stored in packed form
! Q(iv,jv) => qfunc(ijv) as defined below (for jv >= iv) ! Q(iv,jv) => qfunc(ijv) as defined below (for jv >= iv)
! FIXME: no longer
! !
ijv = jv * (jv-1) / 2 + iv ijv = jv * (jv-1) / 2 + iv
read( iunps, '(1p4e19.11)', err=100, iostat=ios ) & read( iunps, '(1p4e19.11)', err=100, iostat=ios ) &
upf%dion(iv,jv), ddd(iv,jv), upf%qqq(iv,jv), & upf%dion(iv,jv), ddd(iv,jv), upf%qqq(iv,jv), &
(upf%qfunc(ir,iv,jv),ir=1,upf%kkbeta), & (upf%qfunc(ir,ijv),ir=1,upf%kkbeta), &
((upf%qfcoef(i,lp,iv,jv),i=1,upf%nqf),lp=1,upf%nqlc) ((upf%qfcoef(i,lp,iv,jv),i=1,upf%nqf),lp=1,upf%nqlc)
do ir=upf%kkbeta+1,upf%mesh do ir=upf%kkbeta+1,upf%mesh
upf%qfunc(ir,iv,jv)=0.0_DP upf%qfunc(ir,ijv)=0.0_DP
enddo enddo
! !
! Use the symmetry of the coefficients ! Use the symmetry of the coefficients
@ -313,19 +315,21 @@ CONTAINS
if ( iv /= jv ) then if ( iv /= jv ) then
upf%dion(jv,iv)=upf%dion(iv,jv) upf%dion(jv,iv)=upf%dion(iv,jv)
upf%qqq(jv,iv) =upf%qqq(iv,jv) upf%qqq(jv,iv) =upf%qqq(iv,jv)
upf%qfunc(:,jv,iv) = upf%qfunc(:,iv,jv)
upf%qfcoef(:,:,jv,iv)=upf%qfcoef(:,:,iv,jv) upf%qfcoef(:,:,jv,iv)=upf%qfcoef(:,:,iv,jv)
end if end if
enddo enddo
enddo enddo
DEALLOCATE (ddd)
! !
! for versions later than 7.2 ! for versions later than 7.2
! !
if (10*iver(1,is)+iver(2,is) >= 72) then if (10*iver(1,is)+iver(2,is) >= 72) then
ALLOCATE (iptype(upf%nbeta))
read( iunps, '(6i5)',err=100, iostat=ios ) & read( iunps, '(6i5)',err=100, iostat=ios ) &
(iptype(iv), iv=1,upf%nbeta) (iptype(iv), iv=1,upf%nbeta)
read( iunps, '(i5,f15.9)',err=100, iostat=ios ) & read( iunps, '(i5,f15.9)',err=100, iostat=ios ) &
npf, dummy npf, dummy
DEALLOCATE (iptype)
end if end if
! !
@ -449,6 +453,7 @@ CONTAINS
WRITE( stdout,1300) WRITE( stdout,1300)
1300 format (4x,60('=')) 1300 format (4x,60('='))
! !
DEALLOCATE (eee)
return return
100 call errore('readvan','error reading pseudo file', abs(ios) ) 100 call errore('readvan','error reading pseudo file', abs(ios) )
! !
@ -490,14 +495,14 @@ CONTAINS
read(iunps,*, err=100) & read(iunps,*, err=100) &
( (qrl(ir,l),ir=1,upf%kkbeta), l=lmin,lmax) ( (qrl(ir,l),ir=1,upf%kkbeta), l=lmin,lmax)
! !
!!! ijv = jv * (jv-1) / 2 + iv ijv = jv * (jv-1) / 2 + iv
! !
do l=lmin,lmax do l=lmin,lmax
! !
! reconstruct rinner ! reconstruct rinner
! !
do ir=upf%kkbeta,1,-1 do ir=upf%kkbeta,1,-1
if ( abs(qrl(ir,l)-upf%qfunc(ir,iv,jv)) > 1.0d-6) go to 10 if ( abs(qrl(ir,l)-upf%qfunc(ir,ijv)) > 1.0d-6) go to 10
end do end do
10 irinner = ir+1 10 irinner = ir+1
upf%rinner(l) = upf%r(irinner) upf%rinner(l) = upf%r(irinner)
@ -669,7 +674,7 @@ CONTAINS
read( iunps, '(2i5)', err=100, iostat=ios ) & read( iunps, '(2i5)', err=100, iostat=ios ) &
upf%nwfc, upf%nbeta upf%nwfc, upf%nbeta
! !
if ( upf%nbeta > nbrx .or. upf%nbeta < 0) & if ( upf%nbeta < 0) &
call errore('readrrkj', 'wrong nbeta', is) call errore('readrrkj', 'wrong nbeta', is)
if ( upf%nwfc > nchix .or. upf%nwfc < 0) & if ( upf%nwfc > nchix .or. upf%nwfc < 0) &
call errore('readrrkj', 'wrong nchi', is) call errore('readrrkj', 'wrong nchi', is)
@ -694,7 +699,7 @@ CONTAINS
ALLOCATE ( upf%kbeta(upf%nbeta) ) ALLOCATE ( upf%kbeta(upf%nbeta) )
ALLOCATE ( upf%dion(upf%nbeta,upf%nbeta), upf%qqq(upf%nbeta,upf%nbeta) ) ALLOCATE ( upf%dion(upf%nbeta,upf%nbeta), upf%qqq(upf%nbeta,upf%nbeta) )
ALLOCATE ( upf%beta(upf%mesh,upf%nbeta) ) ALLOCATE ( upf%beta(upf%mesh,upf%nbeta) )
ALLOCATE ( upf%qfunc(upf%mesh,upf%nbeta,upf%nbeta) ) ALLOCATE ( upf%qfunc(upf%mesh,upf%nbeta*(upf%nbeta+1)/2) )
upf%kkbeta = 0 upf%kkbeta = 0
do nb=1,upf%nbeta do nb=1,upf%nbeta
read ( iunps, '(i6)',err=100, iostat=ios ) upf%kbeta(nb) read ( iunps, '(i6)',err=100, iostat=ios ) upf%kbeta(nb)
@ -708,7 +713,7 @@ CONTAINS
! !
! the symmetric matric Q_{nb,mb} is stored in packed form ! the symmetric matric Q_{nb,mb} is stored in packed form
! Q(nb,mb) => qfunc(ijv) as defined below (for mb <= nb) ! Q(nb,mb) => qfunc(ijv) as defined below (for mb <= nb)
! FIXME: no longer !
ijv = nb * (nb - 1) / 2 + mb ijv = nb * (nb - 1) / 2 + mb
read( iunps, '(1p4e19.11)', err=100, iostat=ios ) & read( iunps, '(1p4e19.11)', err=100, iostat=ios ) &
upf%dion(nb,mb) upf%dion(nb,mb)
@ -716,15 +721,14 @@ CONTAINS
read(iunps,'(1p4e19.11)',err=100,iostat=ios) & read(iunps,'(1p4e19.11)',err=100,iostat=ios) &
upf%qqq(nb,mb) upf%qqq(nb,mb)
read(iunps,'(1p4e19.11)',err=100,iostat=ios) & read(iunps,'(1p4e19.11)',err=100,iostat=ios) &
(upf%qfunc(n,nb,mb),n=1,upf%mesh) (upf%qfunc(n,ijv),n=1,upf%mesh)
else else
upf%qqq(nb,mb)=0.0_DP upf%qqq(nb,mb)=0.0_DP
upf%qfunc(:,nb,mb)=0.0_DP upf%qfunc(:,ijv)=0.0_DP
endif endif
if ( mb /= nb ) then if ( mb /= nb ) then
upf%dion(mb,nb)=upf%dion(nb,mb) upf%dion(mb,nb)=upf%dion(nb,mb)
upf%qqq(mb,nb)=upf%qqq(nb,mb) upf%qqq(mb,nb)=upf%qqq(nb,mb)
upf%qfunc(:,mb,nb)=upf%qfunc(:,nb,mb)
end if end if
enddo enddo
enddo enddo

32
Modules/upf_to_internal.f90
View File

@ -30,8 +30,7 @@ subroutine set_pseudo_upf (is, upf)
! !
USE atom, ONLY: rgrid, chi, oc, nchi, lchi, jchi, rho_at, & USE atom, ONLY: rgrid, chi, oc, nchi, lchi, jchi, rho_at, &
rho_atc, nlcc rho_atc, nlcc
USE uspp_param, ONLY: zp, vloc_at, dion, betar, qqq, qfcoef, qfunc, nqf, & USE uspp_param, ONLY: tvanp
nqlc, rinner, nbeta, kkbeta, lll, jjj, psd, tvanp
USE funct, ONLY: set_dft_from_name, set_dft_from_indices, dft_is_meta USE funct, ONLY: set_dft_from_name, set_dft_from_indices, dft_is_meta
! !
USE pseudo_types USE pseudo_types
@ -47,8 +46,6 @@ subroutine set_pseudo_upf (is, upf)
TYPE (pseudo_upf) :: upf TYPE (pseudo_upf) :: upf
! !
! !
zp(is) = upf%zp
psd (is)= upf%psd
tvanp(is)=upf%tvanp tvanp(is)=upf%tvanp
nlcc(is) = upf%nlcc nlcc(is) = upf%nlcc
! workaround for rrkj format - it contains the indices, not the name ! workaround for rrkj format - it contains the indices, not the name
@ -64,30 +61,6 @@ subroutine set_pseudo_upf (is, upf)
oc(1:upf%nwfc, is) = upf%oc(1:upf%nwfc) oc(1:upf%nwfc, is) = upf%oc(1:upf%nwfc)
chi(1:upf%mesh, 1:upf%nwfc, is) = upf%chi(1:upf%mesh, 1:upf%nwfc) chi(1:upf%mesh, 1:upf%nwfc, is) = upf%chi(1:upf%mesh, 1:upf%nwfc)
! !
nbeta(is)= upf%nbeta
kkbeta(is)=upf%kkbeta
betar(1:upf%mesh, 1:upf%nbeta, is) = upf%beta(1:upf%mesh, 1:upf%nbeta)
dion(1:upf%nbeta, 1:upf%nbeta, is) = upf%dion(1:upf%nbeta, 1:upf%nbeta)
!
nqlc(is) = upf%nqlc
nqf (is) = upf%nqf
lll(1:upf%nbeta,is) = upf%lll(1:upf%nbeta)
rinner(1:upf%nqlc,is) = upf%rinner(1:upf%nqlc)
lll(1:upf%nbeta,is) = upf%lll(1:upf%nbeta)
if ( upf%tvanp ) then
qqq(1:upf%nbeta,1:upf%nbeta,is) = upf%qqq(1:upf%nbeta,1:upf%nbeta)
do nb = 1, upf%nbeta
do mb = nb, upf%nbeta
ijv = mb * (mb-1) / 2 + nb
qfunc (1:upf%mesh, ijv, is) = upf%qfunc(1:upf%mesh, nb, mb)
end do
end do
if ( upf%nqf > 0 ) then
qfcoef(1:upf%nqf, 1:upf%nqlc, 1:upf%nbeta, 1:upf%nbeta, is ) = &
upf%qfcoef( 1:upf%nqf, 1:upf%nqlc, 1:upf%nbeta, 1:upf%nbeta )
end if
end if
!
rgrid(is)%dx = upf%dx rgrid(is)%dx = upf%dx
rgrid(is)%xmin = upf%xmin rgrid(is)%xmin = upf%xmin
rgrid(is)%zmesh= upf%zmesh rgrid(is)%zmesh= upf%zmesh
@ -100,10 +73,8 @@ subroutine set_pseudo_upf (is, upf)
! !
if (upf%has_so) then if (upf%has_so) then
jchi(1:upf%nwfc, is) = upf%jchi(1:upf%nwfc) jchi(1:upf%nwfc, is) = upf%jchi(1:upf%nwfc)
jjj(1:upf%nbeta, is) = upf%jjj(1:upf%nbeta)
else else
jchi(1:upf%nwfc, is) = 0.0_DP jchi(1:upf%nwfc, is) = 0.0_DP
jjj(1:upf%nbeta, is) = 0.0_DP
endif endif
! !
if ( upf%nlcc) then if ( upf%nlcc) then
@ -112,7 +83,6 @@ subroutine set_pseudo_upf (is, upf)
rho_atc(:,is) = 0.0_DP rho_atc(:,is) = 0.0_DP
end if end if
rho_at (1:upf%mesh, is) = upf%rho_at (1:upf%mesh) rho_at (1:upf%mesh, is) = upf%rho_at (1:upf%mesh)
vloc_at(1:upf%mesh,is) = upf%vloc(1:upf%mesh)
end subroutine set_pseudo_upf end subroutine set_pseudo_upf

21
Modules/uspp.f90
View File

@ -10,7 +10,7 @@ MODULE uspp_param
! ... Ultrasoft and Norm-Conserving pseudopotential parameters ! ... Ultrasoft and Norm-Conserving pseudopotential parameters
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
USE parameters, ONLY : lqmax, nbrx, npsx, nqfx USE parameters, ONLY : lqmax, npsx
USE radial_grids, ONLY: ndmx USE radial_grids, ONLY: ndmx
USE pseudo_types, ONLY: pseudo_upf USE pseudo_types, ONLY: pseudo_upf
! !
@ -18,29 +18,10 @@ MODULE uspp_param
! !
TYPE (pseudo_upf), ALLOCATABLE, TARGET :: upf(:) TYPE (pseudo_upf), ALLOCATABLE, TARGET :: upf(:)
CHARACTER(LEN=2 ) :: psd(npsx) ! name of the pseudopotential
REAL(DP) :: &
zp(npsx), &! the charge of the pseudopotential
dion(nbrx,nbrx,npsx), &! D_{mu,nu} parameters (in the atomic case)
betar(ndmx,nbrx,npsx), &! radial beta_{mu} functions
jjj(nbrx,npsx), &! total angular momentum of the beta function
qqq(nbrx,nbrx,npsx), &! q_{mu,nu} parameters (in the atomic case)
qfunc(ndmx,nbrx*(nbrx+1)/2,npsx), &! Q_{mu,nu}(|r|) function for |r|> r_L
qfcoef(nqfx,lqmax,nbrx,nbrx,npsx), &! coefficients for Q for |r|<r_L
! augfun(ndmx,nbrx,nbrx,0:lqmax,npsx), &! moved to grid_paw_variables.f90
vloc_at(ndmx,npsx), &! local potential
rinner(lqmax,npsx) ! values of r_L
INTEGER :: & INTEGER :: &
nbeta(npsx), &! number of beta functions
nh(npsx), &! number of beta functions per atomic type nh(npsx), &! number of beta functions per atomic type
nhm, &! max number of different beta functions per atom nhm, &! max number of different beta functions per atom
nbetam, &! max number of beta functions nbetam, &! max number of beta functions
kkbeta(npsx), &! point where the beta are zero
nqf(npsx), &! number of coefficients for Q
nqlc(npsx), &! number of angular momenta in Q
ifqopt(npsx), &! level of q optimization
lll(nbrx,npsx), &! angular momentum of the beta function
iver(3,npsx) ! version of the atomic code iver(3,npsx) ! version of the atomic code
INTEGER :: & INTEGER :: &
lmaxkb, &! max angular momentum lmaxkb, &! max angular momentum

3
PW/bp_calc_btq.f90
View File

@ -44,8 +44,7 @@ SUBROUTINE calc_btq(ql,qr_k,idbes)
aux(:) = 0.0_DP aux(:) = 0.0_DP
DO i = upf(np)%kkbeta,2,-1 DO i = upf(np)%kkbeta,2,-1
IF (rgrid(np)%r(i) .LT. upf(np)%rinner(l+1)) GOTO 100 IF (rgrid(np)%r(i) .LT. upf(np)%rinner(l+1)) GOTO 100
!!! aux(i) = qfunc(i,ijv,np) TEMP aux(i) = upf(np)%qfunc(i,ijv)
aux(i) = upf(np)%qfunc(i,iv,jv)
ENDDO ENDDO
100 CALL setqf ( upf(np)%qfcoef(1,l+1,iv,jv), aux(1), & 100 CALL setqf ( upf(np)%qfcoef(1,l+1,iv,jv), aux(1), &
rgrid(np)%r, upf(np)%nqf, l, i ) rgrid(np)%r, upf(np)%nqf, l, i )

9
PW/compute_qdipol.f90
View File

@ -47,8 +47,7 @@ SUBROUTINE compute_qdipol(dpqq)
(mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then (mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then
do ir = 1, upf(nt)%kkbeta do ir = 1, upf(nt)%kkbeta
if (rgrid(nt)%r(ir) >= upf(nt)%rinner(l+1)) then if (rgrid(nt)%r(ir) >= upf(nt)%rinner(l+1)) then
! qtot(ir, nb, mb)=qfunc(ir,ijv,nt) TEMP qtot(ir, nb, mb)=upf(nt)%qfunc(ir,ijv)
qtot(ir, nb, mb)=upf(nt)%qfunc(ir,nb,mb)
else else
ilast = ir ilast = ir
endif endif
@ -97,8 +96,10 @@ SUBROUTINE compute_qdipol(dpqq)
nb=indv(jh,nt) nb=indv(jh,nt)
if (ivl > nlx) call errore('compute_qdipol',' ivl > nlx', ivl) if (ivl > nlx) call errore('compute_qdipol',' ivl > nlx', ivl)
if (jvl > nlx) call errore('compute_qdipol',' jvl > nlx', jvl) if (jvl > nlx) call errore('compute_qdipol',' jvl > nlx', jvl)
if (nb > nbetam) call errore('compute_qdipol',' nb > nbrx', nb) if (nb > nbetam) &
if (mb > nbetam) call errore('compute_qdipol',' mb > nbrx', mb) call errore('compute_qdipol',' nb out of bounds', nb)
if (mb > nbetam) &
call errore('compute_qdipol',' mb out of bounds', mb)
if (mb > nb) call errore('compute_qdipol',' mb > nb', 1) if (mb > nb) call errore('compute_qdipol',' mb > nb', 1)
dpqq(ih,jh,ipol,nt)=fact*ap(lp,ivl,jvl)*qrad2(mb,nb,nt) dpqq(ih,jh,ipol,nt)=fact*ap(lp,ivl,jvl)*qrad2(mb,nb,nt)
dpqq(jh,ih,ipol,nt)=dpqq(ih,jh,ipol,nt) dpqq(jh,ih,ipol,nt)=dpqq(ih,jh,ipol,nt)

4
PW/gen_us_dj.f90
View File

@ -24,7 +24,7 @@ subroutine gen_us_dj (ik, dvkb)
USE uspp, ONLY : nkb, indv, nhtol, nhtolm USE uspp, ONLY : nkb, indv, nhtol, nhtolm
USE us, ONLY : nqx, tab, tab_d2y, dq, spline_ps USE us, ONLY : nqx, tab, tab_d2y, dq, spline_ps
USE splinelib USE splinelib
USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh USE uspp_param, ONLY : upf, lmaxkb, nbetam, nh
! !
implicit none implicit none
! !
@ -88,7 +88,7 @@ subroutine gen_us_dj (ik, dvkb)
endif endif
do nt = 1, ntyp do nt = 1, ntyp
do nb = 1, nbeta (nt) do nb = 1, upf(nt)%nbeta
do ig = 1, npw do ig = 1, npw
qt = sqrt(q (ig)) * tpiba qt = sqrt(q (ig)) * tpiba
if (spline_ps) then if (spline_ps) then

4
PW/gen_us_dy.f90
View File

@ -25,7 +25,7 @@ subroutine gen_us_dy (ik, u, dvkb)
USE uspp, ONLY : nkb, indv, nhtol, nhtolm USE uspp, ONLY : nkb, indv, nhtol, nhtolm
USE us, ONLY : nqx, tab, tab_d2y, dq, spline_ps USE us, ONLY : nqx, tab, tab_d2y, dq, spline_ps
USE splinelib USE splinelib
USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh USE uspp_param, ONLY : upf, lmaxkb, nbetam, nh
! !
implicit none implicit none
! !
@ -82,7 +82,7 @@ subroutine gen_us_dy (ik, u, dvkb)
do nt = 1, ntyp do nt = 1, ntyp
! calculate beta in G-space using an interpolation table ! calculate beta in G-space using an interpolation table
do nb = 1, nbeta (nt) do nb = 1, upf(nt)%nbeta
do ig = 1, npw do ig = 1, npw
if (spline_ps) then if (spline_ps) then
vkb0(ig,nb,nt) = splint(xdata, tab(:,nb,nt), & vkb0(ig,nb,nt) = splint(xdata, tab(:,nb,nt), &

11
PW/grid_paw_routines.f90
View File

@ -28,7 +28,6 @@ CONTAINS
SUBROUTINE allocate_paw_internals SUBROUTINE allocate_paw_internals
USE gvect, ONLY : nrxx USE gvect, ONLY : nrxx
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE parameters, ONLY : nbrx
USE radial_grids, ONLY : ndmx USE radial_grids, ONLY : ndmx
USE ions_base, ONLY : nsp, nat, ntyp => nsp USE ions_base, ONLY : nsp, nat, ntyp => nsp
USE us, ONLY : nqxq USE us, ONLY : nqxq
@ -89,7 +88,6 @@ CONTAINS
SUBROUTINE deallocate_paw_internals SUBROUTINE deallocate_paw_internals
USE gvect, ONLY : nrxx USE gvect, ONLY : nrxx
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE parameters, ONLY : nbrx
USE radial_grids, ONLY : ndmx USE radial_grids, ONLY : ndmx
USE ions_base, ONLY : nsp, nat, ntyp => nsp USE ions_base, ONLY : nsp, nat, ntyp => nsp
USE us, ONLY : nqxq USE us, ONLY : nqxq
@ -192,7 +190,7 @@ CONTAINS
qi ! q-point grid for interpolation qi ! q-point grid for interpolation
REAL(DP), ALLOCATABLE :: ylmk0 (:) ! the spherical harmonics REAL(DP), ALLOCATABLE :: ylmk0 (:) ! the spherical harmonics
INTEGER :: n1, m0, m1, n, li, mi, vi, vj, ijs, is1, is2, & INTEGER :: n1, m0, m1, n, li, mi, vi, vj, ijs, is1, is2, &
lk, mk, vk, kh, lh, sph_ind, nnbrx, ll lk, mk, vk, kh, lh, sph_ind, ll
COMPLEX(DP) :: coeff, qgm(1) COMPLEX(DP) :: coeff, qgm(1)
REAL(DP) :: spinor, ji, jk REAL(DP) :: spinor, ji, jk
@ -745,7 +743,8 @@ CONTAINS
!#include "f_defs.h" !#include "f_defs.h"
USE kinds, ONLY: DP USE kinds, ONLY: DP
USE us, ONLY: dq!, qrad USE us, ONLY: dq!, qrad
USE uspp_param, ONLY: lmaxq, nbrx USE parameters, ONLY: nbrx
USE uspp_param, ONLY: lmaxq
USE uspp, ONLY: nlx, lpl, lpx, ap, indv, nhtolm USE uspp, ONLY: nlx, lpl, lpx, ap, indv, nhtolm
IMPLICIT NONE IMPLICIT NONE
@ -1466,7 +1465,7 @@ END SUBROUTINE atomic_becsum
USE cell_base, ONLY : omega, tpiba2 USE cell_base, ONLY : omega, tpiba2
USE gvect, ONLY : ngl, gl USE gvect, ONLY : ngl, gl
USE pseud, ONLY : lloc, lmax, cc, nlc, nnl, alpc, alps, aps USE pseud, ONLY : lloc, lmax, cc, nlc, nnl, alpc, alps, aps
USE uspp_param, ONLY : zp USE uspp_param, ONLY : upf
USE grid_paw_variables, ONLY: aevloc_at, psvloc_at, aevloc, psvloc USE grid_paw_variables, ONLY: aevloc_at, psvloc_at, aevloc, psvloc
! !
USE radial_grids, ONLY : ndmx USE radial_grids, ONLY : ndmx
@ -1498,7 +1497,7 @@ END SUBROUTINE atomic_becsum
! !
CALL vloc_of_g_noerf (lloc (nt), lmax (nt), numeric (nt), rgrid(nt)%mesh, & CALL vloc_of_g_noerf (lloc (nt), lmax (nt), numeric (nt), rgrid(nt)%mesh, &
msh (nt), rgrid(nt)%rab(1), rgrid(nt)%r(1), vloc_at_ (:, nt), cc (1, & msh (nt), rgrid(nt)%rab(1), rgrid(nt)%r(1), vloc_at_ (:, nt), cc (1, &
nt), alpc (1, nt), nlc (nt), nnl (nt), zp (nt), aps (1, 0, nt), & nt), alpc (1, nt), nlc (nt), nnl (nt), upf(nt)%zp, aps (1, 0, nt), &
alps (1, 0, nt), tpiba2, ngl, gl, omega, vloc_ (:, nt) ) alps (1, 0, nt), tpiba2, ngl, gl, omega, vloc_ (:, nt) )
END DO END DO
END DO whattodo END DO whattodo

3
PW/init_us_1.f90
View File

@ -252,8 +252,7 @@ subroutine init_us_1
(mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then (mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then
do ir = 1, upf(nt)%kkbeta do ir = 1, upf(nt)%kkbeta
if (rgrid(nt)%r(ir) >=upf(nt)%rinner (l+1) ) then if (rgrid(nt)%r(ir) >=upf(nt)%rinner (l+1) ) then
! qtot (ir, ijv) = qfunc (ir, ijv, nt) TEMP qtot (ir, ijv) = upf(nt)%qfunc(ir,ijv)
qtot (ir, ijv) = upf(nt)%qfunc(ir,nb,mb)
else else
ilast = ir ilast = ir
endif endif

6
PW/pwscf.f90
View File

@ -12,7 +12,7 @@ PROGRAM pwscf
! ... Plane Wave Self-Consistent Field code ! ... Plane Wave Self-Consistent Field code
! !
USE io_global, ONLY : stdout, ionode USE io_global, ONLY : stdout, ionode
USE parameters, ONLY : ntypx, npk, lmaxx, nchix, nqfx, nbrx USE parameters, ONLY : ntypx, npk, lmaxx, nchix, nbrx
use radial_grids, ONLY : ndmx use radial_grids, ONLY : ndmx
USE global_version, ONLY : version_number USE global_version, ONLY : version_number
USE wvfct, ONLY : gamma_only USE wvfct, ONLY : gamma_only
@ -59,7 +59,7 @@ PROGRAM pwscf
FMT = '(/5X,"Ultrasoft (Vanderbilt) Pseudopotentials")') FMT = '(/5X,"Ultrasoft (Vanderbilt) Pseudopotentials")')
! !
WRITE( unit = stdout, FMT = 9010 ) & WRITE( unit = stdout, FMT = 9010 ) &
ntypx, npk, lmaxx, nchix, ndmx, nbrx, nqfx ntypx, npk, lmaxx, nchix, ndmx, nbrx
! !
END IF END IF
! !
@ -133,6 +133,6 @@ PROGRAM pwscf
! !
9010 FORMAT( /,5X,'Current dimensions of program pwscf are:', /, & 9010 FORMAT( /,5X,'Current dimensions of program pwscf are:', /, &
& /,5X,'ntypx = ',I2,' npk = ',I5,' lmax = ',I2 & & /,5X,'ntypx = ',I2,' npk = ',I5,' lmax = ',I2 &
& /,5X,'nchix = ',I2,' ndmx = ',I5,' nbrx = ',I2,' nqfx = ',I2 ) & /,5X,'nchix = ',I2,' ndmx = ',I5,' nbrx = ',I2 )
! !
END PROGRAM pwscf END PROGRAM pwscf

10
PW/read_file.f90
View File

@ -36,7 +36,7 @@ SUBROUTINE read_file()
USE vlocal, ONLY : strf USE vlocal, ONLY : strf
USE io_files, ONLY : tmp_dir, prefix, iunpun, nwordwfc, iunwfc USE io_files, ONLY : tmp_dir, prefix, iunpun, nwordwfc, iunwfc
USE buffers, ONLY : open_buffer, close_buffer USE buffers, ONLY : open_buffer, close_buffer
USE uspp_param, ONLY : upf, nbeta, jjj USE uspp_param, ONLY : upf
USE noncollin_module, ONLY : noncolin, npol USE noncollin_module, ONLY : noncolin, npol
USE mp_global, ONLY : kunit USE mp_global, ONLY : kunit
USE pw_restart, ONLY : pw_readfile USE pw_restart, ONLY : pw_readfile
@ -155,13 +155,7 @@ SUBROUTINE read_file()
! !
DO nt = 1, nsp DO nt = 1, nsp
! !
so(nt) = ( upf(nt)%nbeta > 0 ) so(nt) = upf(nt)%has_so
!
DO nb = 1, upf(nt)%nbeta
!
so(nt) = so(nt) .AND. ( ABS( upf(nt)%jjj(nb) ) > 1.D-7 )
!
END DO
! !
END DO END DO
! !

6
PW/read_pseudo.f90
View File

@ -82,7 +82,7 @@ subroutine readpp
call read_pseudo_upf(iunps, upf(nt), isupf) call read_pseudo_upf(iunps, upf(nt), isupf)
! !
if (isupf == 0) then if (isupf == 0) then
call set_pseudo_upf (nt, upf(nt)) ! TEMP call set_pseudo_upf (nt, upf(nt))
call set_paw_upf (nt, upf(nt)) call set_paw_upf (nt, upf(nt))
! for compatibility with old formats ! for compatibility with old formats
newpseudo (nt) = .true. newpseudo (nt) = .true.
@ -109,7 +109,7 @@ subroutine readpp
ELSE ELSE
CALL readvan (iunps, nt, upf(nt)) CALL readvan (iunps, nt, upf(nt))
ENDIF ENDIF
CALL set_pseudo_upf (nt, upf(nt)) ! TEMP CALL set_pseudo_upf (nt, upf(nt))
! !
else if (pseudo_type (psfile (nt) ) ==3) then else if (pseudo_type (psfile (nt) ) ==3) then
! !
@ -130,7 +130,7 @@ subroutine readpp
! !
call read_ncpp (iunps, nt, upf(nt)) call read_ncpp (iunps, nt, upf(nt))
! !
CALL set_pseudo_upf (nt, upf(nt)) ! TEMP CALL set_pseudo_upf (nt, upf(nt))
! !
endif endif
! for compatibility with old formats - maybe obsolete? ! for compatibility with old formats - maybe obsolete?

13
PW/realus.f90
View File

@ -72,7 +72,7 @@ MODULE realus
IMPLICIT NONE IMPLICIT NONE
! !
INTEGER :: qsdim, ia, mbia, iqs, iqsia INTEGER :: qsdim, ia, mbia, iqs, iqsia
INTEGER :: indm, inbrx, idimension, & INTEGER :: indm, idimension, &
ih, jh, ijh, lllnbnt, lllmbnt ih, jh, ijh, lllnbnt, lllmbnt
INTEGER :: roughestimate, goodestimate, lamx2, l, nt INTEGER :: roughestimate, goodestimate, lamx2, l, nt
INTEGER, ALLOCATABLE :: buffpoints(:,:) INTEGER, ALLOCATABLE :: buffpoints(:,:)
@ -106,13 +106,10 @@ MODULE realus
boxrad(:) = 0.D0 boxrad(:) = 0.D0
! !
DO nt = 1, nsp DO nt = 1, nsp
!DO inbrx = 1, upf(nt)%nbeta*(upf(nt)%nbeta+1)/2 TEMP DO ijv = 1, upf(nt)%nbeta*(upf(nt)%nbeta+1)/2
DO nb = 1, upf(nt)%nbeta
DO mb = nb, upf(nt)%nbeta
DO indm = upf(nt)%kkbeta, 1, -1 DO indm = upf(nt)%kkbeta, 1, -1
! !
! IF ( ABS( qfunc(indm,inbrx,nt) ) > eps16 ) THEN TEMP IF ( ABS( upf(nt)%qfunc(indm,ijv) ) > eps16 ) THEN
IF ( ABS( upf(nt)%qfunc(indm,nb,mb) ) > eps16 ) THEN
! !
boxrad(nt) = MAX( rgrid(nt)%r(indm), boxrad(nt) ) boxrad(nt) = MAX( rgrid(nt)%r(indm), boxrad(nt) )
! !
@ -122,7 +119,6 @@ MODULE realus
! !
END DO END DO
END DO END DO
END DO
END DO END DO
! !
boxrad(:) = boxrad(:) / alat boxrad(:) = boxrad(:) / alat
@ -368,8 +364,7 @@ MODULE realus
! !
DO ir = 1, upf(nt)%kkbeta DO ir = 1, upf(nt)%kkbeta
IF ( rgrid(nt)%r(ir) >= upf(nt)%rinner(l+1) ) THEN IF ( rgrid(nt)%r(ir) >= upf(nt)%rinner(l+1) ) THEN
!qtot(ir,nb,mb) = qfunc(ir,ijv,nt) / rgrid(nt)%r(ir)**2 qtot(ir,nb,mb) = upf(nt)%qfunc(ir,ijv) / &
qtot(ir,nb,mb) = upf(nt)%qfunc(ir,nb,mb) / &
rgrid(nt)%r(ir)**2 rgrid(nt)%r(ir)**2
ELSE ELSE
ilast = ir ilast = ir