From 6d97829a147e9a66046101a00006caf4f99119fe Mon Sep 17 00:00:00 2001 From: giannozz Date: Sat, 6 Oct 2007 08:23:39 +0000 Subject: [PATCH] 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 --- CPV/cprsub.f90 | 15 ++- CPV/modules.f90 | 2 - CPV/pseudopot_sub.f90 | 7 +- CPV/qqberry.f90 | 26 ++--- CPV/read_pseudo.f90 | 178 +-------------------------------- CPV/stress.f90 | 6 +- Modules/grid_paw_variables.f90 | 2 +- Modules/parameters.f90 | 4 +- Modules/pseudo_types.f90 | 55 +++++----- Modules/read_paw.f90 | 5 +- Modules/read_upf.f90 | 16 ++- Modules/read_uspp.f90 | 48 +++++---- Modules/upf_to_internal.f90 | 32 +----- Modules/uspp.f90 | 21 +--- PW/bp_calc_btq.f90 | 3 +- PW/compute_qdipol.f90 | 9 +- PW/gen_us_dj.f90 | 4 +- PW/gen_us_dy.f90 | 4 +- PW/grid_paw_routines.f90 | 11 +- PW/init_us_1.f90 | 3 +- PW/pwscf.f90 | 6 +- PW/read_file.f90 | 10 +- PW/read_pseudo.f90 | 6 +- PW/realus.f90 | 13 +-- 24 files changed, 123 insertions(+), 363 deletions(-) diff --git a/CPV/cprsub.f90 b/CPV/cprsub.f90 index 9d9240045..2891f9705 100644 --- a/CPV/cprsub.f90 +++ b/CPV/cprsub.f90 @@ -27,7 +27,7 @@ subroutine formf( tfirst, eself ) use ions_base, ONLY : rcmax, zv, nsp, na use local_pseudo, ONLY : vps, rhops, dvps, drhops 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 pseudopotential, ONLY : tpstab, vps_sp, dvps_sp use cp_interfaces, ONLY : build_pstab @@ -83,9 +83,9 @@ subroutine formf( tfirst, eself ) if ( numeric(is) ) then - call formfn( vps(:,is), dvps(:,is), rgrid(is)%r, rgrid(is)%rab, vloc_at(:,is), & - zv(is), rcmax(is), g, omega, tpiba2, rgrid(is)%mesh, & - ngs, oldvan(is), tpre ) + call formfn( vps(:,is), dvps(:,is), rgrid(is)%r, rgrid(is)%rab, & + upf(is)%vloc(1), zv(is), rcmax(is), g, omega, tpiba2,& + rgrid(is)%mesh, ngs, oldvan(is), tpre ) else @@ -332,8 +332,7 @@ subroutine nlinit use constants, ONLY : pi, fpi use ions_base, ONLY : na, nsp use uspp, ONLY : aainit, beta, qq, dvan, nhtol, nhtolm, indv - use uspp_param, ONLY : kkbeta, qqq, nqlc, betar, lmaxq, dion,& - nbeta, nbetam, lmaxkb, lll, nhm, nh, tvanp + use uspp_param, ONLY : upf, lmaxq, nbetam, lmaxkb, nhm, nh use atom, ONLY : rgrid, nlcc, numeric use qradb_mod, ONLY : qradb use qgb_mod, ONLY : qgb @@ -410,8 +409,8 @@ subroutine nlinit ! WRITE( stdout,*) WRITE( stdout,'(20x,a)') ' dion ' - do iv = 1, nbeta(is) - WRITE( stdout,'(8f9.4)') ( fac*dion(iv,jv,is), jv = 1, nbeta(is) ) + do iv = 1, upf(is)%nbeta + WRITE( stdout,'(8f9.4)') ( fac*upf(is)%dion(iv,jv), jv = 1, upf(is)%nbeta ) end do ! end do diff --git a/CPV/modules.f90 b/CPV/modules.f90 index f3b0ee5d1..bce547ad5 100644 --- a/CPV/modules.f90 +++ b/CPV/modules.f90 @@ -53,8 +53,6 @@ end module bhs ! ! 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 ! diff --git a/CPV/pseudopot_sub.f90 b/CPV/pseudopot_sub.f90 index f6d1af297..4ae0a33e6 100644 --- a/CPV/pseudopot_sub.f90 +++ b/CPV/pseudopot_sub.f90 @@ -389,7 +389,7 @@ USE pseudo_base, ONLY : nlin_base USE pseudo_base, ONLY : nlin_stress_base 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 cp_interfaces, ONLY : compute_xgtab, chkpstab USE pseudopotential, ONLY : wnl_sp, wnla_sp, xgtab @@ -561,7 +561,7 @@ USE cell_base, ONLY : tpiba USE pseudo_base, ONLY : nlin_stress_base USE splines, ONLY : spline - USE uspp_param, ONLY : upf, nbeta + USE uspp_param, ONLY : upf USE read_pseudo_module_fpmd, ONLY : nspnl USE reciprocal_vectors, ONLY : g, gstart USE pseudopotential, ONLY : wnla_sp, tpstab @@ -1624,8 +1624,7 @@ do l = lmin, lmax do ir = 1, upf(is)%kkbeta 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,iv,jv) + qrl(ir,ijv,l)=upf(is)%qfunc(ir,ijv) else qrl(ir,ijv,l)=upf(is)%qfcoef(1,l,iv,jv) do i = 2, upf(is)%nqf diff --git a/CPV/qqberry.f90 b/CPV/qqberry.f90 index 328b02808..32d15b88c 100644 --- a/CPV/qqberry.f90 +++ b/CPV/qqberry.f90 @@ -17,7 +17,7 @@ subroutine qqberry2( gqq,gqqm, ipol) use smallbox_grid_dimensions, only: nr1b, nr2b, nr3b, & 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 atom, only: rgrid use core @@ -50,7 +50,7 @@ subroutine qqberry2( gqq,gqqm, ipol) integer :: ivs, jvs, ivl, jvl, lp, ijv real(8), allocatable:: ylm(:,:) - ndm = MAXVAL (kkbeta(1:nsp)) + ndm = MAXVAL (upf(1:nsp)%kkbeta) allocate( fint( ndm), jl(ndm)) allocate( qradb2(nbetam,nbetam,lmaxq,nsp)) allocate( ylm(ngw, lmaxq*lmaxq)) @@ -86,28 +86,30 @@ subroutine qqberry2( gqq,gqqm, ipol) do is=1,nvb 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 ) ! now the radial part - do l=1,nqlc(is) + do l=1,upf(is)%nqlc xg= gmes !only orthorombic cells - !!!call bess(xg,l,kkbeta(is),r(1,is),jl) - call sph_bes ( kkbeta(is), rgrid(is)%r, xg, l-1, jl ) - do iv= 1,nbeta(is) - do jv=iv,nbeta(is) + !!!call bess(xg,l,upf(is)%kkbeta,rgrid(is)%r,jl) + call sph_bes ( upf(is)%kkbeta, rgrid(is)%r, xg, l-1, jl ) + do iv= 1,upf(is)%nbeta + do jv=iv,upf(is)%nbeta ijv = (jv-1)*jv/2 + iv ! ! 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) end do if (oldvan(is)) then - call herman_skillman_int & - & (kkbeta(is),fint,rgrid(is)%rab,qradb2(iv,jv,l,is)) + call herman_skillman_int ( upf(is)%kkbeta,fint,rgrid(is)%rab,& + qradb2(iv,jv,l,is) ) 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 qradb2(iv,jv,l,is)= c*qradb2(iv,jv,l,is) if ( iv /= jv ) qradb2(jv,iv,l,is)= qradb2(iv,jv,l,is) diff --git a/CPV/read_pseudo.f90 b/CPV/read_pseudo.f90 index 8753651d1..8406ce3b6 100644 --- a/CPV/read_pseudo.f90 +++ b/CPV/read_pseudo.f90 @@ -319,7 +319,7 @@ END FUNCTION calculate_dx IF( program_name == 'FPMD' ) THEN CALL errore(' readpp ', ' file format not supported ', 1 ) ELSE - call readbhs(is,pseudounit) + CALL errore(' readpp ', ' file format no longer supported ', 2 ) END IF END IF @@ -533,179 +533,3 @@ END FUNCTION calculate_dx 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 diff --git a/CPV/stress.f90 b/CPV/stress.f90 index bcab97d70..63b12aa32 100644 --- a/CPV/stress.f90 +++ b/CPV/stress.f90 @@ -60,7 +60,7 @@ USE cell_base, ONLY: tpiba2, omega USE reciprocal_vectors, ONLY: gstart, gzero, g, gx 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 electrons_base, only: nupdwn, iupdwn, nspin USE cdvan, only: dbec @@ -150,8 +150,8 @@ ! ... initialize array wnla ALLOCATE( wsg ( nhm, nsp ) ) - ALLOCATE( wnl ( ngw, MAXVAL( nbeta( 1:nsp ) ), nsp, 1 ) ) - ALLOCATE( wnla( ngw, MAXVAL( nbeta( 1:nsp ) ), nsp ) ) + ALLOCATE( wnl ( ngw, nbetam, nsp, 1 ) ) + ALLOCATE( wnla( ngw, nbetam, nsp ) ) ALLOCATE( fnlb( nsanl, MAXVAL( nh( 1:nsp ) ), MAXVAL( nupdwn ), nspin ) ) ! fac = sqrt( omega ) / ( 2.0d0 * 4.0d0 * pi ) diff --git a/Modules/grid_paw_variables.f90 b/Modules/grid_paw_variables.f90 index 84a070ffb..2adcdd5bc 100644 --- a/Modules/grid_paw_variables.f90 +++ b/Modules/grid_paw_variables.f90 @@ -8,7 +8,7 @@ module grid_paw_variables ! NO rinner > 0 ! USE kinds, ONLY : DP - USE parameters, ONLY : lqmax, nbrx, npsx, nqfx + USE parameters, ONLY : lqmax, nbrx, npsx USE radial_grids, ONLY: ndmx ! implicit none diff --git a/Modules/parameters.f90 b/Modules/parameters.f90 index 26cffe325..85d1bb902 100644 --- a/Modules/parameters.f90 +++ b/Modules/parameters.f90 @@ -20,14 +20,12 @@ MODULE parameters npsx = ntypx, &! max number of different PPs (obsolete) npk = 40000, &! max number of k-points 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 nwfsx = 14 ! max number of beta functions per atom ! INTEGER, PARAMETER :: & nbrx = 14, &! max number of beta functions - lqmax= 2*lmaxx+1, &! max number of angular momenta of Q - nqfx = 8 ! max number of coefficients in Q smoothing + lqmax= 2*lmaxx+1 ! max number of angular momenta of Q ! ! ... More parameter for the CP codes diff --git a/Modules/pseudo_types.f90 b/Modules/pseudo_types.f90 index f7cfd7d2e..ec477bb6d 100644 --- a/Modules/pseudo_types.f90 +++ b/Modules/pseudo_types.f90 @@ -12,8 +12,7 @@ ! together with their allocation/deallocation routines USE kinds, ONLY: DP - USE parameters, ONLY: cp_lmax, lmaxx - use radial_grids, ONLY: ndmx, radial_grid_type + use radial_grids, ONLY: radial_grid_type IMPLICIT NONE SAVE @@ -81,44 +80,46 @@ END TYPE paw_t REAL(DP) :: rmax ! the maximum radius of the mesh REAL(DP) :: zmesh ! the nuclear charge used for 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 :: rcutus(:)! cut-off ultrasoft radius (nbeta) + REAL(DP), POINTER :: rcutus(:)! ultrasoft cut-off radius (nbeta) REAL(DP), POINTER :: epseu(:) ! energy (nwfc) - REAL(DP), POINTER :: jchi(:) ! jchi(nwfc) - REAL(DP), POINTER :: jjj(:) ! jjj(nbeta) + REAL(DP), POINTER :: jchi(:) ! jchi(nwfc) j=l+1/2 or l-1/2 of wfc + REAL(DP), POINTER :: jjj(:) ! jjj(nbeta) j=l+1/2 or l-1/2 of beta 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 :: nwfc ! number of wavefunctions + INTEGER :: nwfc ! number of atomic wavefunctions INTEGER :: nbeta ! number of projectors INTEGER :: kkbeta ! kkbeta=max(kbeta(:)) ! kbeta<=mesh is the number of grid points for each beta function ! beta(r,nb) = 0 for r > r(kbeta(nb)) ! kkbeta<=mesh is the largest of such number so that for all beta ! beta(r,nb) = 0 for r > r(kkbeta) - CHARACTER(LEN=2), POINTER :: els(:) ! els(nwfc) - CHARACTER(LEN=2), POINTER :: els_beta(:) ! els(nbeta) - INTEGER, POINTER :: lchi(:) ! lchi(nwfc) - REAL(DP), POINTER :: oc(:) ! oc(nwfc) - REAL(DP), POINTER :: r(:) ! r(mesh) - REAL(DP), POINTER :: rab(:) ! rab(mesh) - REAL(DP), POINTER :: rho_atc(:) ! rho_atc(mesh) - REAL(DP), POINTER :: vloc(:) ! vloc(mesh) - INTEGER, POINTER :: lll(:) ! lll(nbeta) - INTEGER, POINTER :: kbeta(:) ! kbeta(nbeta) - REAL(DP), POINTER :: beta(:,:) ! beta(mesh,nbeta) + CHARACTER(LEN=2), POINTER :: els(:) ! els(nwfc) label of wfc + CHARACTER(LEN=2), POINTER :: els_beta(:) ! els(nbeta) label of beta + INTEGER, POINTER :: lchi(:) ! lchi(nwfc) value of l for wavefcts + REAL(DP), POINTER :: oc(:) ! oc(nwfc) occupancies for wavefcts + REAL(DP), POINTER :: r(:) ! r(mesh) radial grid + REAL(DP), POINTER :: rab(:) ! rab(mesh) dr(x)/dx (x=linear grid) + REAL(DP), POINTER :: rho_atc(:) ! rho_atc(mesh) atomic core charge + REAL(DP), POINTER :: vloc(:) ! vloc(mesh) local atomic potential + INTEGER, POINTER :: lll(:) ! lll(nbeta) l of each projector + INTEGER, POINTER :: kbeta(:) ! kbeta(nbeta) see above kkbeta + REAL(DP), POINTER :: beta(:,:) ! beta(mesh,nbeta) projectors INTEGER :: nd - REAL(DP), POINTER :: dion(:,:) ! dion(nbeta,nbeta) - INTEGER :: nqf - INTEGER :: nqlc - REAL(DP), POINTER :: rinner(:) ! rinner(0:2*lmax) - REAL(DP), POINTER :: qqq(:,:) ! qqq(nbeta,nbeta) - REAL(DP), POINTER :: qfunc(:,:,:) ! qfunc(mesh,nbeta,nbeta) + REAL(DP), POINTER :: dion(:,:) ! dion(nbeta,nbeta) atomic D_{mu,nu} + INTEGER :: nqf ! number of Q coefficients + INTEGER :: nqlc ! number of angular momenta in Q + REAL(DP), POINTER :: rinner(:) ! rinner(0:2*lmax) r_L + REAL(DP), POINTER :: qqq(:,:) ! qqq(nbeta,nbeta) q_{mu,nu} + 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 :: chi(:,:) ! chi(mesh,nwfc) - REAL(DP), POINTER :: rho_at(:) ! rho_at(mesh) + ! coefficients for Q for |r| 0 ) then call scan_begin (iunps, "QFCOEF", .false.) @@ -448,7 +446,7 @@ subroutine read_pseudo_nl (upf, iunps) upf%nqlc = 2 * upf%lmax + 1 ALLOCATE( upf%rinner( upf%nqlc ) ) 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 ) ) upf%rinner = 0.0_DP upf%qqq = 0.0_DP diff --git a/Modules/read_uspp.f90 b/Modules/read_uspp.f90 index 8e46a7e5d..b6f903a70 100644 --- a/Modules/read_uspp.f90 +++ b/Modules/read_uspp.f90 @@ -13,7 +13,7 @@ MODULE read_uspp_module ! Vanderbilt's code and Andrea's RRKJ3 format ! 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 funct, ONLY: set_dft_from_name, dft_is_hybrid, dft_is_meta, & set_dft_from_indices @@ -101,13 +101,16 @@ CONTAINS & exfact, &! index of the exchange and correlation used & etotpseu, &! total pseudopotential energy & ee(nchix), &! the energy of the valence states + & rc(nchix), &! the cut-off radii of the pseudopotential & eloc, &! energy of the local potential & 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 & 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 & & idmy(3), &! contains the date of creation of the pseudo & ifpcor, &! for core correction, 0 otherwise @@ -115,9 +118,8 @@ CONTAINS & i, &! dummy counter & nnlz(nchix), &! The nlm values of the valence states & 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 - & iptype(nbrx), &! more recent parameters & npf, &! as above & nang, &! number of angular momenta in pseudopotentials & lloc, &! angular momentum of the local part of PPs @@ -213,7 +215,7 @@ CONTAINS if ( lloc == -1 ) lloc = nang+1 if ( lloc > nang+1 .or. lloc < 0 ) & 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) if ( ifqopt < 0 ) & call errore( 'readvan', 'wrong ifqopt read', is ) @@ -272,8 +274,8 @@ CONTAINS ALLOCATE ( upf%kbeta(upf%nbeta) ) upf%kbeta(:) = upf%kkbeta ! - if( upf%nbeta > nbrx .or. upf%nbeta <0 ) & - call errore( 'readvan','nbeta wrong or too large', is ) + if( upf%nbeta < 0 ) & + call errore( 'readvan','nbeta wrong', is ) if( upf%kkbeta > upf%mesh .or. upf%kkbeta < 0 ) & call errore( 'readvan','kkbeta wrong or too large', is ) ! @@ -282,8 +284,9 @@ CONTAINS ALLOCATE ( upf%lll(upf%nbeta) ) ALLOCATE ( upf%beta(upf%mesh,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 ( eee(upf%nbeta), ddd(upf%nbeta,upf%nbeta) ) do iv=1,upf%nbeta read( iunps, '(i5)',err=100, iostat=ios ) upf%lll(iv) read( iunps, '(1p4e19.11)',err=100, iostat=ios ) & @@ -297,15 +300,14 @@ CONTAINS ! ! the symmetric matric Q_{nb,mb} is stored in packed form ! Q(iv,jv) => qfunc(ijv) as defined below (for jv >= iv) - ! FIXME: no longer ! ijv = jv * (jv-1) / 2 + iv read( iunps, '(1p4e19.11)', err=100, iostat=ios ) & 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) do ir=upf%kkbeta+1,upf%mesh - upf%qfunc(ir,iv,jv)=0.0_DP + upf%qfunc(ir,ijv)=0.0_DP enddo ! ! Use the symmetry of the coefficients @@ -313,19 +315,21 @@ CONTAINS if ( iv /= jv ) then upf%dion(jv,iv)=upf%dion(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) end if enddo enddo + DEALLOCATE (ddd) ! ! for versions later than 7.2 ! if (10*iver(1,is)+iver(2,is) >= 72) then + ALLOCATE (iptype(upf%nbeta)) read( iunps, '(6i5)',err=100, iostat=ios ) & (iptype(iv), iv=1,upf%nbeta) read( iunps, '(i5,f15.9)',err=100, iostat=ios ) & npf, dummy + DEALLOCATE (iptype) end if ! @@ -449,6 +453,7 @@ CONTAINS WRITE( stdout,1300) 1300 format (4x,60('=')) ! + DEALLOCATE (eee) return 100 call errore('readvan','error reading pseudo file', abs(ios) ) ! @@ -490,14 +495,14 @@ CONTAINS read(iunps,*, err=100) & ( (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 ! ! reconstruct rinner ! 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 10 irinner = ir+1 upf%rinner(l) = upf%r(irinner) @@ -669,7 +674,7 @@ CONTAINS read( iunps, '(2i5)', err=100, iostat=ios ) & upf%nwfc, upf%nbeta ! - if ( upf%nbeta > nbrx .or. upf%nbeta < 0) & + if ( upf%nbeta < 0) & call errore('readrrkj', 'wrong nbeta', is) if ( upf%nwfc > nchix .or. upf%nwfc < 0) & call errore('readrrkj', 'wrong nchi', is) @@ -694,7 +699,7 @@ CONTAINS ALLOCATE ( upf%kbeta(upf%nbeta) ) ALLOCATE ( upf%dion(upf%nbeta,upf%nbeta), upf%qqq(upf%nbeta,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 do nb=1,upf%nbeta 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 ! Q(nb,mb) => qfunc(ijv) as defined below (for mb <= nb) - ! FIXME: no longer + ! ijv = nb * (nb - 1) / 2 + mb read( iunps, '(1p4e19.11)', err=100, iostat=ios ) & upf%dion(nb,mb) @@ -716,15 +721,14 @@ CONTAINS read(iunps,'(1p4e19.11)',err=100,iostat=ios) & upf%qqq(nb,mb) 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 upf%qqq(nb,mb)=0.0_DP - upf%qfunc(:,nb,mb)=0.0_DP + upf%qfunc(:,ijv)=0.0_DP endif if ( mb /= nb ) then upf%dion(mb,nb)=upf%dion(nb,mb) upf%qqq(mb,nb)=upf%qqq(nb,mb) - upf%qfunc(:,mb,nb)=upf%qfunc(:,nb,mb) end if enddo enddo diff --git a/Modules/upf_to_internal.f90 b/Modules/upf_to_internal.f90 index e11a51259..0a6bfb490 100644 --- a/Modules/upf_to_internal.f90 +++ b/Modules/upf_to_internal.f90 @@ -30,8 +30,7 @@ subroutine set_pseudo_upf (is, upf) ! USE atom, ONLY: rgrid, chi, oc, nchi, lchi, jchi, rho_at, & rho_atc, nlcc - USE uspp_param, ONLY: zp, vloc_at, dion, betar, qqq, qfcoef, qfunc, nqf, & - nqlc, rinner, nbeta, kkbeta, lll, jjj, psd, tvanp + USE uspp_param, ONLY: tvanp USE funct, ONLY: set_dft_from_name, set_dft_from_indices, dft_is_meta ! USE pseudo_types @@ -47,8 +46,6 @@ subroutine set_pseudo_upf (is, upf) TYPE (pseudo_upf) :: upf ! ! - zp(is) = upf%zp - psd (is)= upf%psd tvanp(is)=upf%tvanp nlcc(is) = upf%nlcc ! 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) 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)%xmin = upf%xmin rgrid(is)%zmesh= upf%zmesh @@ -100,10 +73,8 @@ subroutine set_pseudo_upf (is, upf) ! if (upf%has_so) then jchi(1:upf%nwfc, is) = upf%jchi(1:upf%nwfc) - jjj(1:upf%nbeta, is) = upf%jjj(1:upf%nbeta) else jchi(1:upf%nwfc, is) = 0.0_DP - jjj(1:upf%nbeta, is) = 0.0_DP endif ! if ( upf%nlcc) then @@ -112,7 +83,6 @@ subroutine set_pseudo_upf (is, upf) rho_atc(:,is) = 0.0_DP end if 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 diff --git a/Modules/uspp.f90 b/Modules/uspp.f90 index 9cde92121..59472fdca 100644 --- a/Modules/uspp.f90 +++ b/Modules/uspp.f90 @@ -10,7 +10,7 @@ MODULE uspp_param ! ... Ultrasoft and Norm-Conserving pseudopotential parameters ! USE kinds, ONLY : DP - USE parameters, ONLY : lqmax, nbrx, npsx, nqfx + USE parameters, ONLY : lqmax, npsx USE radial_grids, ONLY: ndmx USE pseudo_types, ONLY: pseudo_upf ! @@ -18,29 +18,10 @@ MODULE uspp_param ! 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|= upf(nt)%rinner(l+1)) then - ! qtot(ir, nb, mb)=qfunc(ir,ijv,nt) TEMP - qtot(ir, nb, mb)=upf(nt)%qfunc(ir,nb,mb) + qtot(ir, nb, mb)=upf(nt)%qfunc(ir,ijv) else ilast = ir endif @@ -97,8 +96,10 @@ SUBROUTINE compute_qdipol(dpqq) nb=indv(jh,nt) if (ivl > nlx) call errore('compute_qdipol',' ivl > nlx', ivl) if (jvl > nlx) call errore('compute_qdipol',' jvl > nlx', jvl) - if (nb > nbetam) call errore('compute_qdipol',' nb > nbrx', nb) - if (mb > nbetam) call errore('compute_qdipol',' mb > nbrx', mb) + if (nb > nbetam) & + 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) dpqq(ih,jh,ipol,nt)=fact*ap(lp,ivl,jvl)*qrad2(mb,nb,nt) dpqq(jh,ih,ipol,nt)=dpqq(ih,jh,ipol,nt) diff --git a/PW/gen_us_dj.f90 b/PW/gen_us_dj.f90 index 735dfe770..072a8c7bb 100644 --- a/PW/gen_us_dj.f90 +++ b/PW/gen_us_dj.f90 @@ -24,7 +24,7 @@ subroutine gen_us_dj (ik, dvkb) USE uspp, ONLY : nkb, indv, nhtol, nhtolm USE us, ONLY : nqx, tab, tab_d2y, dq, spline_ps USE splinelib - USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh + USE uspp_param, ONLY : upf, lmaxkb, nbetam, nh ! implicit none ! @@ -88,7 +88,7 @@ subroutine gen_us_dj (ik, dvkb) endif do nt = 1, ntyp - do nb = 1, nbeta (nt) + do nb = 1, upf(nt)%nbeta do ig = 1, npw qt = sqrt(q (ig)) * tpiba if (spline_ps) then diff --git a/PW/gen_us_dy.f90 b/PW/gen_us_dy.f90 index 4db481563..c02b4c94c 100644 --- a/PW/gen_us_dy.f90 +++ b/PW/gen_us_dy.f90 @@ -25,7 +25,7 @@ subroutine gen_us_dy (ik, u, dvkb) USE uspp, ONLY : nkb, indv, nhtol, nhtolm USE us, ONLY : nqx, tab, tab_d2y, dq, spline_ps USE splinelib - USE uspp_param, ONLY : lmaxkb, nbeta, nbetam, nh + USE uspp_param, ONLY : upf, lmaxkb, nbetam, nh ! implicit none ! @@ -82,7 +82,7 @@ subroutine gen_us_dy (ik, u, dvkb) do nt = 1, ntyp ! calculate beta in G-space using an interpolation table - do nb = 1, nbeta (nt) + do nb = 1, upf(nt)%nbeta do ig = 1, npw if (spline_ps) then vkb0(ig,nb,nt) = splint(xdata, tab(:,nb,nt), & diff --git a/PW/grid_paw_routines.f90 b/PW/grid_paw_routines.f90 index 4f6976663..0e5e51ea5 100644 --- a/PW/grid_paw_routines.f90 +++ b/PW/grid_paw_routines.f90 @@ -28,7 +28,6 @@ CONTAINS SUBROUTINE allocate_paw_internals USE gvect, ONLY : nrxx USE lsda_mod, ONLY : nspin - USE parameters, ONLY : nbrx USE radial_grids, ONLY : ndmx USE ions_base, ONLY : nsp, nat, ntyp => nsp USE us, ONLY : nqxq @@ -89,7 +88,6 @@ CONTAINS SUBROUTINE deallocate_paw_internals USE gvect, ONLY : nrxx USE lsda_mod, ONLY : nspin - USE parameters, ONLY : nbrx USE radial_grids, ONLY : ndmx USE ions_base, ONLY : nsp, nat, ntyp => nsp USE us, ONLY : nqxq @@ -192,7 +190,7 @@ CONTAINS qi ! q-point grid for interpolation REAL(DP), ALLOCATABLE :: ylmk0 (:) ! the spherical harmonics 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) REAL(DP) :: spinor, ji, jk @@ -745,7 +743,8 @@ CONTAINS !#include "f_defs.h" USE kinds, ONLY: DP 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 IMPLICIT NONE @@ -1466,7 +1465,7 @@ END SUBROUTINE atomic_becsum USE cell_base, ONLY : omega, tpiba2 USE gvect, ONLY : ngl, gl 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 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, & 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) ) END DO END DO whattodo diff --git a/PW/init_us_1.f90 b/PW/init_us_1.f90 index 6b53ba182..b5f9409a0 100644 --- a/PW/init_us_1.f90 +++ b/PW/init_us_1.f90 @@ -252,8 +252,7 @@ subroutine init_us_1 (mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) ) then do ir = 1, upf(nt)%kkbeta 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,nb,mb) + qtot (ir, ijv) = upf(nt)%qfunc(ir,ijv) else ilast = ir endif diff --git a/PW/pwscf.f90 b/PW/pwscf.f90 index 53ca94450..141d49ccc 100644 --- a/PW/pwscf.f90 +++ b/PW/pwscf.f90 @@ -12,7 +12,7 @@ PROGRAM pwscf ! ... Plane Wave Self-Consistent Field code ! 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 global_version, ONLY : version_number USE wvfct, ONLY : gamma_only @@ -59,7 +59,7 @@ PROGRAM pwscf FMT = '(/5X,"Ultrasoft (Vanderbilt) Pseudopotentials")') ! WRITE( unit = stdout, FMT = 9010 ) & - ntypx, npk, lmaxx, nchix, ndmx, nbrx, nqfx + ntypx, npk, lmaxx, nchix, ndmx, nbrx ! END IF ! @@ -133,6 +133,6 @@ PROGRAM pwscf ! 9010 FORMAT( /,5X,'Current dimensions of program pwscf are:', /, & & /,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 diff --git a/PW/read_file.f90 b/PW/read_file.f90 index 777f19420..501d8d8c4 100644 --- a/PW/read_file.f90 +++ b/PW/read_file.f90 @@ -36,7 +36,7 @@ SUBROUTINE read_file() USE vlocal, ONLY : strf USE io_files, ONLY : tmp_dir, prefix, iunpun, nwordwfc, iunwfc 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 mp_global, ONLY : kunit USE pw_restart, ONLY : pw_readfile @@ -155,13 +155,7 @@ SUBROUTINE read_file() ! DO nt = 1, nsp ! - so(nt) = ( upf(nt)%nbeta > 0 ) - ! - DO nb = 1, upf(nt)%nbeta - ! - so(nt) = so(nt) .AND. ( ABS( upf(nt)%jjj(nb) ) > 1.D-7 ) - ! - END DO + so(nt) = upf(nt)%has_so ! END DO ! diff --git a/PW/read_pseudo.f90 b/PW/read_pseudo.f90 index 4b2ed8f87..e8755fa35 100644 --- a/PW/read_pseudo.f90 +++ b/PW/read_pseudo.f90 @@ -82,7 +82,7 @@ subroutine readpp call read_pseudo_upf(iunps, upf(nt), isupf) ! 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)) ! for compatibility with old formats newpseudo (nt) = .true. @@ -109,7 +109,7 @@ subroutine readpp ELSE CALL readvan (iunps, nt, upf(nt)) ENDIF - CALL set_pseudo_upf (nt, upf(nt)) ! TEMP + CALL set_pseudo_upf (nt, upf(nt)) ! else if (pseudo_type (psfile (nt) ) ==3) then ! @@ -130,7 +130,7 @@ subroutine readpp ! call read_ncpp (iunps, nt, upf(nt)) ! - CALL set_pseudo_upf (nt, upf(nt)) ! TEMP + CALL set_pseudo_upf (nt, upf(nt)) ! endif ! for compatibility with old formats - maybe obsolete? diff --git a/PW/realus.f90 b/PW/realus.f90 index 7fadfd05b..3f754b206 100644 --- a/PW/realus.f90 +++ b/PW/realus.f90 @@ -72,7 +72,7 @@ MODULE realus IMPLICIT NONE ! INTEGER :: qsdim, ia, mbia, iqs, iqsia - INTEGER :: indm, inbrx, idimension, & + INTEGER :: indm, idimension, & ih, jh, ijh, lllnbnt, lllmbnt INTEGER :: roughestimate, goodestimate, lamx2, l, nt INTEGER, ALLOCATABLE :: buffpoints(:,:) @@ -106,13 +106,10 @@ MODULE realus boxrad(:) = 0.D0 ! DO nt = 1, nsp - !DO inbrx = 1, upf(nt)%nbeta*(upf(nt)%nbeta+1)/2 TEMP - DO nb = 1, upf(nt)%nbeta - DO mb = nb, upf(nt)%nbeta + DO ijv = 1, upf(nt)%nbeta*(upf(nt)%nbeta+1)/2 DO indm = upf(nt)%kkbeta, 1, -1 ! - ! IF ( ABS( qfunc(indm,inbrx,nt) ) > eps16 ) THEN TEMP - IF ( ABS( upf(nt)%qfunc(indm,nb,mb) ) > eps16 ) THEN + IF ( ABS( upf(nt)%qfunc(indm,ijv) ) > eps16 ) THEN ! boxrad(nt) = MAX( rgrid(nt)%r(indm), boxrad(nt) ) ! @@ -122,7 +119,6 @@ MODULE realus ! END DO END DO - END DO END DO ! boxrad(:) = boxrad(:) / alat @@ -368,8 +364,7 @@ MODULE realus ! DO ir = 1, upf(nt)%kkbeta 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,nb,mb) / & + qtot(ir,nb,mb) = upf(nt)%qfunc(ir,ijv) / & rgrid(nt)%r(ir)**2 ELSE ilast = ir