diff --git a/CPV/bessel.f90 b/CPV/bessel.f90 index 09f1c53ff..9f1c840cf 100644 --- a/CPV/bessel.f90 +++ b/CPV/bessel.f90 @@ -116,7 +116,7 @@ ARG_S(1:mmax) = RW(1:mmax) * XG xrgm1 = 1.0d0 / arg_s - LMAX = MAXVAL( INDL ) + LMAX = MAXVAL( INDL ) + 1 IF ( LMAX > 0 ) THEN ! ... Calculate J0(|G||r|) = SIN(|G||r|) / (|G||r|) @@ -140,16 +140,16 @@ DO L = 1,LNL LL = INDL(L) - IF(LL == 1) THEN + IF(LL == 0) THEN ! ... FINT = FUNT * J0 FINT(1:mmax,L) = J0(1:mmax) - ELSE IF (LL == 2) THEN + ELSE IF (LL == 1) THEN ! ... FINT = FUNT * J1 FINT(1:mmax,L) = J1(1:mmax) - ELSE IF (LL == 3) THEN + ELSE IF (LL == 2) THEN ! ... FINT = FUNT * J2 FINT(1:mmax,L) = J2(1:mmax) - ELSE IF (LL == 4) THEN + ELSE IF (LL == 3) THEN ! ... FINT = FUNT * J3 FINT(1:mmax,L) = J3(1:mmax) ELSE @@ -195,7 +195,7 @@ ! ... Subroutine Body - LMAX = MAXVAL( INDL ) + LMAX = MAXVAL( INDL ) + 1 IF( ABS(XG) < small ) THEN CALL errore( ' bessel3 ',' xg too small ', 1) @@ -241,13 +241,13 @@ DO L = 1,LNL LL = INDL(L) - IF(LL.EQ.1) THEN + IF(LL.EQ.0) THEN FINT(1:mmax, L) = F0(1:mmax) - ELSE IF (LL.EQ.2) THEN + ELSE IF (LL.EQ.1) THEN FINT(1:mmax, L) = F1(1:mmax) - ELSE IF (LL.EQ.3) THEN + ELSE IF (LL.EQ.2) THEN FINT(1:mmax, L) = F2(1:mmax) - ELSE IF (LL.EQ.4) THEN + ELSE IF (LL.EQ.3) THEN FINT(1:mmax, L) = F3(1:mmax) ELSE CALL errore(" bessel3 "," ll value not programmed ", MAX( 1, ABS(ll) ) ) diff --git a/CPV/cg_sub.f90 b/CPV/cg_sub.f90 index 8d706a2f1..61aaef0b6 100644 --- a/CPV/cg_sub.f90 +++ b/CPV/cg_sub.f90 @@ -51,15 +51,13 @@ use cell_base, only: h, hold, deth, wmass, tpiba2 use smooth_grid_dimensions, only: nnrsx, nr1s, nr2s, nr3s use smallbox_grid_dimensions, only: nnrb => nnrbx, nr1b, nr2b, nr3b - use pseu, only: vps, rhops - use pseu, only: deallocate_pseu + use local_pseudo, only: vps, rhops use io_global, ONLY: io_global_start, stdout, ionode use mp_global, ONLY: mp_global_start use mp, ONLY: mp_end use para_mod use dener use derho - use dpseu use cdvan use stre use parameters, only: nacx, natx, nsx, nbndxx diff --git a/CPV/cp_fpmd.f90 b/CPV/cp_fpmd.f90 index cf210ca78..30c46e467 100644 --- a/CPV/cp_fpmd.f90 +++ b/CPV/cp_fpmd.f90 @@ -1419,7 +1419,7 @@ END SUBROUTINE gshcount tv0rd, taurdr, nv0rd, nbeg, tcp, tcap, & program_name USE ions_base, ONLY: tau_srt, tau_units, if_pos, ind_srt, nsp, na, & - pmass, nat, fricp, greasp + pmass, nat, fricp, greasp, rcmax USE ions_nose, ONLY: tempw, ndega USE constants, ONLY: scmass @@ -1460,7 +1460,7 @@ END SUBROUTINE gshcount WRITE(stdout,660) isa = 0 DO IS = 1, nsp - WRITE(stdout,1000) is, na(is), pmass(is), pmass(is) / scmass + WRITE(stdout,1000) is, na(is), pmass(is), pmass(is) / scmass, rcmax(is) DO IA = 1, na(is) isa = isa + 1 WRITE(stdout,1010) ( tau_srt(k,isa), K = 1,3 ) @@ -1549,7 +1549,8 @@ END SUBROUTINE gshcount 850 FORMAT( 3X,'Initial ion velocities read from unit : ',I4) - 1000 FORMAT(3X,'Species ',I3,' atoms = ',I4,' mass = ',F12.2, ' (a.u.), ', F12.2, ' (uma)' ) + 1000 FORMAT(3X,'Species ',I3,' atoms = ',I4,' mass = ',F12.2, ' (a.u.), ', & + & F12.2, ' (uma)', ' rcmax = ', F6.2, ' (a.u.)' ) 1010 FORMAT(3X,3(1X,F12.6)) 1020 FORMAT(/,3X,'NOT all atoms are allowed to move ') 1021 FORMAT(/,3X,'All atoms are allowed to move') diff --git a/CPV/cplib.f90 b/CPV/cplib.f90 index 0123782e5..b60cff13e 100644 --- a/CPV/cplib.f90 +++ b/CPV/cplib.f90 @@ -502,8 +502,7 @@ use reciprocal_vectors, only: gstart, gx, g use cell_base, only: omega use cell_base, only: ainv, tpiba2 - use pseu - use dpseu + use local_pseudo, only: rhops, drhops use mp, only: mp_sum implicit none @@ -634,8 +633,7 @@ use reciprocal_vectors, only: gstart, gx use cell_base, only: omega use cell_base, only: ainv, tpiba2 - use pseu, only: vps - use dpseu, only: dvps + use local_pseudo, only: vps, dvps use mp, only: mp_sum implicit none @@ -1340,7 +1338,7 @@ use ions_base, only: nsp, na, nas => nax, nat use grid_dimensions, only: nr1, nr2, nr3 use parameters, only: nsx, natx - use pseu + use local_pseudo, only: vps, rhops ! implicit none ! input @@ -3013,7 +3011,7 @@ nr1sx, nr2sx, nr3sx, nnrsx use electrons_base, only: nx => nbspx, n => nbsp, f, ispin => fspin, nspin use constants, only: pi, fpi - use pseu + ! use local_pseudo ! use cdvan use dener @@ -4108,7 +4106,6 @@ use cell_base, only: a1, a2, a3, tpiba2 use reciprocal_vectors, only: gstart, g use recvecs_indexes, only: np, nm - !use parm use grid_dimensions, only: nr1, nr2, nr3, & nr1x, nr2x, nr3x, nnr => nnrx use smooth_grid_dimensions, only: nr1s, nr2s, nr3s, & @@ -4116,13 +4113,11 @@ use electrons_base, only: nspin use constants, only: pi, fpi use energies, only: etot, eself, enl, ekin, epseu, esr, eht, exc - use pseu + use local_pseudo, only: vps, rhops use core use gvecb -! use dener use derho - use dpseu use mp, only: mp_sum ! implicit none @@ -4514,7 +4509,7 @@ use electrons_base, only: nspin, qbac use constants, only: pi, fpi use energies, only: etot, eself, enl, ekin, epseu, esr, eht, exc - use pseu + use local_pseudo, only: rhops, vps use core use gvecb use atom, only: nlcc @@ -4525,7 +4520,6 @@ ! use dener use derho - use dpseu ! implicit none ! diff --git a/CPV/cpr.f90 b/CPV/cpr.f90 index 79fc28b5e..4cdaf7e9e 100644 --- a/CPV/cpr.f90 +++ b/CPV/cpr.f90 @@ -103,13 +103,12 @@ SUBROUTINE cprmain( tau, fion_out, etot_out ) USE grid_dimensions, ONLY : nnrx, nr1, nr2, nr3 USE smooth_grid_dimensions, ONLY : nnrsx, nr1s, nr2s, nr3s USE smallbox_grid_dimensions, ONLY : nr1b, nr2b, nr3b - USE pseu, ONLY : vps, rhops + USE local_pseudo, ONLY : allocate_local_pseudo USE io_global, ONLY : io_global_start, stdout, ionode USE mp_global, ONLY : mp_global_start USE mp, ONLY : mp_sum, mp_barrier USE dener, ONLY : detot USE derho, ONLY : drhor, drhog - USE dpseu, ONLY : dvps, drhops USE cdvan, ONLY : dbec, drhovan USE stre, ONLY : stress USE gvecw, ONLY : ggp @@ -279,14 +278,13 @@ SUBROUTINE cprmain( tau, fion_out, etot_out ) ALLOCATE( c0( ngw, nbspx, 1, 1 ) ) ALLOCATE( cm( ngw, nbspx, 1, 1 ) ) ALLOCATE( phi( ngw, nbspx, 1, 1 ) ) - ALLOCATE( rhops( ngs, nsp ) ) - ALLOCATE( vps( ngs, nsp ) ) + ! + CALL allocate_local_pseudo( ngs, nsp ) + ! ALLOCATE( vkb( ngw, nkb ) ) ALLOCATE( deeq( nhm, nhm, nat, nspin ) ) ALLOCATE( becsum( nhm*(nhm+1)/2, nat, nspin ) ) ALLOCATE( dbec( nkb, nbsp, 3, 3 ) ) - ALLOCATE( dvps( ngs, nsp ) ) - ALLOCATE( drhops( ngs, nsp ) ) ALLOCATE( drhog( ngm, nspin, 3, 3 ) ) ALLOCATE( drhor( nnrx, nspin, 3, 3 ) ) ALLOCATE( drhovan( nhm*(nhm+1)/2, nat, nspin, 3, 3 ) ) diff --git a/CPV/cprsub.f90 b/CPV/cprsub.f90 index e3de4b2e9..8366b6f43 100644 --- a/CPV/cprsub.f90 +++ b/CPV/cprsub.f90 @@ -126,8 +126,7 @@ use cell_base, only: omega, tpiba2 use ions_base, only: rcmax, zv, nsp, na use cvan, only: oldvan - use pseu, only: vps, rhops - use dpseu, only: dvps, drhops + use local_pseudo, only: vps, rhops, dvps, drhops use atom, only: r, rab, mesh, numeric use uspp_param, only: vloc_at use qrl_mod, only: cmesh @@ -563,7 +562,7 @@ use io_global, only: stdout use gvecw, only: ngw use cvan, only: ish, nvb, oldvan - use core + use core, only: rhocb, nlcc_any use constants, only: pi, fpi use ions_base, only: na, nsp use uspp, only: aainit, beta, qq, dvan, nhtol, nhtolm, indv, & @@ -572,13 +571,16 @@ nbeta, lmaxkb, lll, nhm, nh, tvanp use qrl_mod, only: qrl, cmesh use atom, only: mesh, r, rab, nlcc, numeric - use qradb_mod - use qgb_mod - use gvecb - use cdvan - use dqrad_mod - use dqgb_mod - use betax + use qradb_mod, only: qradb + use qgb_mod, only: qgb + use gvecb, only: ngb + use cdvan, only: dbeta + use dqrad_mod, only: dqrad + use dqgb_mod, only: dqgb + use betax, only: qradx, dqradx, refg, betagx, mmx, dbetagx + use pseudopotential, only: pseudopotential_indexes, compute_dvan, & + compute_betagx, compute_qradx + ! implicit none ! @@ -586,308 +588,70 @@ real(kind=8), allocatable:: fint(:), jl(:), jltmp(:), djl(:), & & dfint(:) real(kind=8) xg, xrg, fac -! ------------------------------------------------------------------ -! find number of beta functions per species, max dimensions, -! total number of beta functions (all and Vanderbilt only) -! ------------------------------------------------------------------ - lmaxkb=-1 - nhm=0 - nhsa=0 - nhsavb=0 - nlcc_any=.false. - do is=1,nsp - ind=0 - do iv=1,nbeta(is) - lmaxkb = max(lmaxkb,lll(iv,is)) - ind=ind+2*lll(iv,is)+1 - end do - nh(is)=ind - nhm=max(nhm,nh(is)) - ish(is)=nhsa - nhsa=nhsa+na(is)*nh(is) - if(tvanp(is)) nhsavb=nhsavb+na(is)*nh(is) - nlcc_any = nlcc_any .OR. nlcc(is) - end do - if (lmaxkb > lmaxx) call errore('nlinit ',' l > lmax ',lmaxkb) - lmaxq = 2*lmaxkb + 1 + ! - ! the following prevents an out-of-bound error: nqlc(is)=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 + ! initialize indexes ! - do is=1,nsp - nqlc(is) = MIN ( nqlc(is), lmaxq ) - end do - if (nhsa <= 0) call errore('nlinit ','not implemented ?',nhsa) -! -! initialize array ap -! - call aainit(lmaxkb+1) + CALL pseudopotential_indexes() + ! + ! initialize array ap + ! + call aainit( lmaxkb + 1 ) ! allocate(beta(ngw,nhm,nsp)) allocate(qradb(ngb,nbrx,nbrx,lmaxq,nsp)) allocate(qgb(ngb,nhm*(nhm+1)/2,nsp)) allocate(qq(nhm,nhm,nsp)) - allocate(dvan(nhm,nhm,nsp)) if (nlcc_any) allocate(rhocb(ngb,nsp)) - allocate(nhtol(nhm,nsp)) - allocate(indv (nhm,nsp)) - allocate(nhtolm(nhm,nsp)) ! allocate(dqrad(ngb,nbrx,nbrx,lmaxq,nsp,3,3)) allocate(dqgb(ngb,nhm*(nhm+1)/2,nsp,3,3)) allocate(dbeta(ngw,nhm,nsp,3,3)) - allocate(betagx(mmx,nhm,nsp)) - allocate(dbetagx(mmx,nhm,nsp)) - allocate(qradx(mmx,nbrx,nbrx,lmaxq,nsp)) - allocate(dqradx(mmx,nbrx,nbrx,lmaxq,nsp)) ! qradb(:,:,:,:,:) = 0.d0 qq (:,:,:) =0.d0 - dvan(:,:,:) =0.d0 if(tpre) dqrad(:,:,:,:,:,:,:) = 0.d0 -! -! ------------------------------------------------------------------ -! definition of indices nhtol, indv, nhtolm -! ------------------------------------------------------------------ - do is=1,nsp - ind=0 - do iv=1,nbeta(is) - lm = lll(iv,is)**2 - do il=1,2*lll(iv,is)+1 - lm=lm+1 - ind=ind+1 - nhtolm(ind,is)=lm - nhtol(ind,is)=lll(iv,is) - indv(ind,is)=iv - end do - end do - end do -! -! =============================================================== -! initialization for vanderbilt species -! =============================================================== - do is=1,nvb - if (tpre) then - allocate(dfint(kkbeta(is))) - allocate(djl(kkbeta(is))) - allocate(jltmp(kkbeta(is))) - end if - allocate(fint(kkbeta(is))) - allocate(jl(kkbeta(is))) -! -! qqq and beta are now indexed and taken in the same order -! as vanderbilts ppot-code prints them out -! -! --------------------------------------------------------------- -! calculation of array qradx(igb,iv,jv,is) -! --------------------------------------------------------------- - WRITE( stdout,*) ' nlinit nh(is),ngb,is,kkbeta,lmaxq = ', & - & nh(is),ngb,is,kkbeta(is),nqlc(is) - do l=1,nqlc(is) - do il=1,mmx - xg=sqrt(refg*(il-1)) - call sph_bes (kkbeta(is), r(1,is), xg, l-1, jl) -! - if(tpre) then - ltmp=l-1 - ! - ! r(i0) is the first point such that r(i0) >0 - ! - i0 = 1 - if ( r(1,is) < 1.0d-8 ) i0 = 2 - ! special case q=0 - if (xg < 1.0d-8) then - if (l == 1) then - ! Note that dj_1/dx (x=0) = 1/3 - jltmp(:) = 1.0d0/3.d0 - else - jltmp(:) = 0.0d0 - end if - else - call sph_bes & - (kkbeta(is)+1-i0, r(i0,is), xg, ltmp-1, jltmp ) - end if - do ir=i0, kkbeta(is) - xrg=r(ir,is)*xg - djl(ir)=jltmp(ir)*xrg-l*jl(ir) - end do - if (i0.eq.2) djl(1) = djl(2) - endif -! - do iv= 1,nbeta(is) - do jv=iv,nbeta(is) -! -! note qrl(r)=r^2*q(r) -! - do ir=1,kkbeta(is) - fint(ir)=qrl(ir,iv,jv,l,is)*jl(ir) - end do - if (oldvan(is)) then - call herman_skillman_int & - & (kkbeta(is),cmesh(is),fint,qradx(il,iv,jv,l,is)) - else - call simpson_cp90 & - & (kkbeta(is),fint,rab(1,is),qradx(il,iv,jv,l,is)) - end if - qradx(il,jv,iv,l,is)=qradx(il,iv,jv,l,is) -! - if(tpre) then - do ir=1,kkbeta(is) - dfint(ir)=qrl(ir,iv,jv,l,is)*djl(ir) - end do - if (oldvan(is)) then - call herman_skillman_int & - & (kkbeta(is),cmesh(is),dfint, & - & dqradx(il,iv,jv,l,is)) - else - call simpson_cp90 & - & (kkbeta(is),dfint,rab(1,is), & - & dqradx(il,iv,jv,l,is)) - end if - end if -! - end do - end do - end do - end do -! - WRITE( stdout,*) - WRITE( stdout,'(20x,a)') ' qqq ' - do iv=1,nbeta(is) - WRITE( stdout,'(8f9.4)') (qqq(iv,jv,is),jv=1,nbeta(is)) - end do - WRITE( stdout,*) -! - deallocate(jl) - deallocate(fint) - if (tpre) then - deallocate(jltmp) - deallocate(djl) - deallocate(dfint) - end if -! - end do -! -! =============================================================== -! initialization that is common to all species -! =============================================================== + ! + ! initialization for vanderbilt species + ! + CALL compute_qradx( tpre ) + ! + ! initialization that is common to all species + ! WRITE( stdout, fmt="(//,3X,'Common initialization' )" ) - do is=1,nsp + + do is = 1, nsp WRITE( stdout, fmt="(/,3X,'Specie: ',I5)" ) is - if (.not.numeric(is)) then + if ( .not. numeric(is) ) then fac=1.0 else -! fac converts ry to hartree + ! fac converts ry to hartree fac=0.5 end if - if (tpre) then - allocate(dfint(kkbeta(is))) - allocate(djl(kkbeta(is))) - end if - allocate(fint(kkbeta(is))) - allocate(jl(kkbeta(is))) - allocate(jltmp(kkbeta(is))) -! --------------------------------------------------------------- -! calculation of array betagx(ig,iv,is) -! --------------------------------------------------------------- - WRITE( stdout,*) ' betagx ' - do iv=1,nh(is) - l=nhtol(iv,is)+1 - do il=1,mmx - xg=sqrt(refg*(il-1)) - call sph_bes (kkbeta(is), r(1,is), xg, l-1, jl ) -! - if(tpre)then - ltmp=l-1 - ! - ! r(i0) is the first point such that r(i0) >0 - ! - i0 = 1 - if ( r(1,is) < 1.0d-8 ) i0 = 2 - ! special case q=0 - if (xg < 1.0d-8) then - if (l == 1) then - ! Note that dj_1/dx (x=0) = 1/3 - jltmp(:) = 1.0d0/3.d0 - else - jltmp(:) = 0.0d0 - end if - else - call sph_bes & - (kkbeta(is)+1-i0, r(i0,is), xg, ltmp-1, jltmp ) - end if - do ir=i0, kkbeta(is) - xrg=r(ir,is)*xg - djl(ir)=jltmp(ir)*xrg-l*jl(ir) - end do - if (i0.eq.2) djl(1) = djl(2) -! - endif -! -! beta(ir)=r*beta(r) -! - do ir=1,kkbeta(is) - fint(ir)=r(ir,is)*betar(ir,indv(iv,is),is)*jl(ir) - end do - if (oldvan(is)) then - call herman_skillman_int & - & (kkbeta(is),cmesh(is),fint,betagx(il,iv,is)) - else - call simpson_cp90 & - & (kkbeta(is),fint,rab(1,is),betagx(il,iv,is)) - endif -! - if(tpre) then - do ir=1,kkbeta(is) - dfint(ir)=r(ir,is)*betar(ir,indv(iv,is),is)*djl(ir) - end do - if (oldvan(is)) then - call herman_skillman_int & - & (kkbeta(is),cmesh(is),dfint,dbetagx(il,iv,is)) - else - call simpson_cp90 & - & (kkbeta(is),dfint,rab(1,is),dbetagx(il,iv,is)) - end if - endif -! - end do - end do -! -! --------------------------------------------------------------- -! calculate array dvan(iv,jv,is) -! --------------------------------------------------------------- - do iv=1,nh(is) - do jv=1,nh(is) - if ( nhtolm(iv,is) == nhtolm(jv,is) ) then - dvan(iv,jv,is)=fac*dion(indv(iv,is),indv(jv,is),is) - endif - end do - end do -! - do iv=1,nh(is) - WRITE( stdout,901) iv,indv(iv,is),nhtol(iv,is) + do iv = 1, nh(is) + WRITE( stdout,901) iv, indv(iv,is), nhtol(iv,is) end do 901 format(2x,i2,' indv= ',i2,' ang. mom= ',i2) -! + ! 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, nbeta(is) + WRITE( stdout,'(8f9.4)') ( fac * dion(iv,jv,is), jv = 1, nbeta(is) ) end do -! - deallocate(jltmp) - deallocate(jl) - deallocate(fint) - if (tpre) then - deallocate(djl) - deallocate(dfint) - end if + ! end do -! -! newnlinit stores qgb and qq, calculates arrays beta qradb rhocb -! and derivatives wrt cell dbeta dqrad -! + ! + ! calculation of array betagx(ig,iv,is) + ! + call compute_betagx( tpre ) + ! + ! calculate array dvan(iv,jv,is) + ! + call compute_dvan() + ! + ! newnlinit stores qgb and qq, calculates arrays beta qradb rhocb + ! and derivatives wrt cell dbeta dqrad + ! call newnlinit return diff --git a/CPV/cptypes.f90 b/CPV/cptypes.f90 index b48266c63..3aed86605 100644 --- a/CPV/cptypes.f90 +++ b/CPV/cptypes.f90 @@ -41,11 +41,6 @@ !! ... pseudopotential TYPE pseudo -!! ... local part - REAL(dbl), POINTER :: vps(:,:) ! form factors - REAL(dbl), POINTER :: dvps(:,:) ! cell derivatives of form factors - REAL(dbl), POINTER :: rhops(:,:) ! Ewald pseudocharges form factors - REAL(dbl), POINTER :: drhops(:,:) ! Ewald pseudocharges form factors !! ... nonlocal part REAL(dbl), POINTER :: wsg(:,:) ! inverse of Kleinman-Bylander ! denominators @@ -91,14 +86,6 @@ LOGICAL, INTENT(IN) :: tcc(:) INTEGER :: ierr - ALLOCATE(ps%vps(ng,nsp), STAT=ierr) - IF( ierr /= 0 ) CALL errore(' allocate_pseudo ', ' allocating %vps ', ierr ) - ALLOCATE(ps%dvps(ng,nsp), STAT=ierr) - IF( ierr /= 0 ) CALL errore(' allocate_pseudo ', ' allocating %dvps ', ierr ) - ALLOCATE(ps%rhops(ng,nsp), STAT=ierr) - IF( ierr /= 0 ) CALL errore(' allocate_pseudo ', ' allocating %rhops ', ierr ) - ALLOCATE(ps%drhops(ng,nsp), STAT=ierr) - IF( ierr /= 0 ) CALL errore(' allocate_pseudo ', ' allocating %drhops ', ierr ) ALLOCATE(ps%wnl(ngw,lnlx,nsp,nk), STAT=ierr) IF( ierr /= 0 ) CALL errore(' allocate_pseudo ', ' allocating %wnl ', ierr ) ALLOCATE(ps%wsg(ngh,nsp), STAT=ierr) @@ -133,22 +120,6 @@ TYPE (pseudo) ps INTEGER :: ierr - IF(ASSOCIATED(ps%vps)) THEN - DEALLOCATE(ps%vps, STAT=ierr) - IF( ierr /= 0 ) CALL errore(' deallocate_pseudo ', ' deallocating %vps ', ierr ) - END IF - IF(ASSOCIATED(ps%dvps)) THEN - DEALLOCATE(ps%dvps, STAT=ierr) - IF( ierr /= 0 ) CALL errore(' deallocate_pseudo ', ' deallocating %dvps ', ierr ) - END IF - IF(ASSOCIATED(ps%rhops)) THEN - DEALLOCATE(ps%rhops, STAT=ierr) - IF( ierr /= 0 ) CALL errore(' deallocate_pseudo ', ' deallocating %rhops ', ierr ) - END IF - IF(ASSOCIATED(ps%drhops)) THEN - DEALLOCATE(ps%drhops, STAT=ierr) - IF( ierr /= 0 ) CALL errore(' deallocate_pseudo ', ' deallocating %drhops ', ierr ) - END IF IF(ASSOCIATED(ps%wnl)) THEN DEALLOCATE(ps%wnl, STAT=ierr) IF( ierr /= 0 ) CALL errore(' deallocate_pseudo ', ' deallocating %wnl ', ierr ) diff --git a/CPV/dealloc.f90 b/CPV/dealloc.f90 index c5459ed7d..55531c15b 100644 --- a/CPV/dealloc.f90 +++ b/CPV/dealloc.f90 @@ -22,7 +22,7 @@ use cg_module, only: deallocate_cg use reciprocal_vectors, only: deallocate_recvecs use recvecs_indexes, only: deallocate_recvecs_indexes - use pseu, only: deallocate_pseu + use local_pseudo, only: deallocate_local_pseudo use qgb_mod, only: deallocate_qgb_mod use dqgb_mod, only: deallocate_dqgb_mod use qradb_mod, only: deallocate_qradb_mod @@ -43,7 +43,6 @@ USE turbo, ONLY: deallocate_turbo USE cp_main_variables, ONLY: deallocate_mainvar USE derho, ONLY: deallocate_derho - USE dpseu, ONLY: deallocate_dpseu USE cdvan, ONLY: deallocate_cdvan USE pseudopotential, ONLY: deallocate_pseudopotential @@ -59,12 +58,11 @@ CALL deallocate_uspp() CALL deallocate_recvecs() CALL deallocate_recvecs_indexes() - CALL deallocate_pseu() + CALL deallocate_local_pseudo() CALL deallocate_qgb_mod() CALL deallocate_qradb_mod() CALL deallocate_derho() CALL deallocate_dqgb_mod() - CALL deallocate_dpseu() CALL deallocate_cdvan() CALL deallocate_dqrad_mod() CALL deallocate_betax() diff --git a/CPV/diis.f90 b/CPV/diis.f90 index 15119be7e..c0adb8b12 100644 --- a/CPV/diis.f90 +++ b/CPV/diis.f90 @@ -313,12 +313,15 @@ ! ... declare modules USE control_flags, ONLY: gamma_only USE cell_base, ONLY: tpiba2 - USE pseudopotential, ONLY: tl, l2ind, nspnl, lm1x + USE pseudopotential, ONLY: nspnl USE ions_base, ONLY: nsp, na USE mp_global, ONLY: group USE mp, ONLY: mp_sum, mp_max USE reciprocal_vectors, ONLY: gstart, gzero, ggp USE reciprocal_space_mesh, ONLY: gkmask_l, gkcutz_l + USE uspp_param, only: nh + USE uspp, only: nhtol, indv + IMPLICIT NONE @@ -332,7 +335,7 @@ ! ... declare other variables REAL(dbl) vp,ftpi,arg - INTEGER l,ll,i,ig,igh,is,m,j, ikk + INTEGER l,ll,i,ig,igh,is,m,j, ikk, ih, iv ! ... end of declarations ! ---------------------------------------------- @@ -360,19 +363,12 @@ vpp = vp ! ... nonlocal potential - igh = 0 - DO l = 0, lm1x - IF(tl(l)) THEN - DO m = -l, l - igh = igh + 1 - DO is = 1, nspnl - ll = l2ind( l+1, is) - IF( ll > 0 ) THEN - vpp(:) = vpp(:) + na(is) * wsg(igh,is) * wnl(:,ll,is,ikk)**2 - END IF - END DO - END DO - END IF + DO is = 1, nspnl + DO igh = 1, nh( is ) + ll = nhtol ( ih, is ) + 1 + iv = indv ( ih, is ) + vpp(:) = vpp(:) + na(is) * wsg(igh,is) * wnl(:,iv,is,ikk)**2 + END DO END DO ! ... kinetic energy diff --git a/CPV/electrons.f90 b/CPV/electrons.f90 index 0c53fb3ae..98e9ee466 100644 --- a/CPV/electrons.f90 +++ b/CPV/electrons.f90 @@ -10,6 +10,7 @@ !=----------------------------------------------------------------------------=! USE kinds + USE parameters, ONLY: nspinx USE parallel_toolkit, ONLY: pdspev_drv, dspev_drv, pzhpev_drv, zhpev_drv USE electrons_base, ONLY: nbnd, nbndx, nbsp, nbspx, nspin, nel, nelt, & nupdwn, iupdwn, telectrons_base_initval, f @@ -29,8 +30,8 @@ INTEGER :: n_emp = 0 ! number of empty states - INTEGER :: nb_l(2) = 0 ! local number of states ( for each spin components ) - INTEGER :: n_emp_l(2) = 0 + INTEGER :: nb_l(nspinx) = 0 ! local number of states ( for each spin components ) + INTEGER :: n_emp_l(nspinx) = 0 INTEGER, ALLOCATABLE :: ib_owner(:) INTEGER, ALLOCATABLE :: ib_local(:) @@ -109,13 +110,16 @@ DO ik = 1, nk occ( 1:nbnd, ik, 1 ) = f( 1:nbnd ) END DO - ELSE + ELSE IF( nspin == 2 ) THEN DO ik = 1, nk occ( 1:nupdwn(1), ik, 1 ) = f( 1:nupdwn(1) ) END DO DO ik = 1, nk occ( 1:nupdwn(2), ik, 2 ) = f( iupdwn(2) : ( iupdwn(2) + nupdwn(2) - 1 ) ) END DO + ELSE + WRITE( stdout, * ) ' nspin = ', nspin + CALL errore(' band_init ',' nspin not implemented ', 3) END IF IF( ionode ) THEN @@ -154,7 +158,7 @@ END IF DO i = 1, nspin - IF( i > 1 ) CALL errore( ' bmeshset ',' spin too large ', i) + IF( i > nspinx ) CALL errore( ' bmeshset ',' spin too large ', i) nb_l( i ) = nupdwn( i ) / nproc IF( mpime < MOD( nupdwn( i ), nproc ) ) nb_l( i ) = nb_l( i ) + 1 n_emp_l( i ) = n_emp / nproc diff --git a/CPV/emptystates.f90 b/CPV/emptystates.f90 index c531705a5..3b7ed2a11 100644 --- a/CPV/emptystates.f90 +++ b/CPV/emptystates.f90 @@ -427,7 +427,7 @@ USE wave_types, ONLY: wave_descriptor USE wave_functions, ONLY: cp_kinetic_energy, crot, dft_kinetic_energy, fixwave USE wave_base, ONLY: hpsi, converg_base, dotp - USE pseudopotential, ONLY: nsanl, ngh + USE pseudopotential, ONLY: nsanl USE constants, ONLY: au USE cell_base, ONLY: tpiba2 USE electrons_module, ONLY: pmss, n_emp @@ -449,6 +449,7 @@ USE control_flags, ONLY: force_pairing, gamma_only USE reciprocal_space_mesh, ONLY: gkmask_l, gkx_l, gk_l USE reciprocal_vectors, ONLY: g, gx + USE uspp_param, ONLY: nhm IMPLICIT NONE @@ -493,7 +494,7 @@ ekinc_old = 1.d+10 ekinc = 0.0d0 - CALL allocate_projector(fnle, nsanl, n_emp, ngh, gamma) + CALL allocate_projector(fnle, nsanl, n_emp, nhm, gamma) ALLOCATE( eforce( ngw, wempt%ldb, nk, nspin ) ) ALLOCATE( fi( wempt%ldb, nk, nspin ) ) ALLOCATE( cp_emp( SIZE(c_emp,1), SIZE(c_emp,2), SIZE(c_emp,3), SIZE(c_emp,4) ) ) diff --git a/CPV/exch_corr.f90 b/CPV/exch_corr.f90 index 302386fa4..9196ad443 100644 --- a/CPV/exch_corr.f90 +++ b/CPV/exch_corr.f90 @@ -256,13 +256,28 @@ REAL (dbl) :: v2xc(:,:,:,:,:) REAL(dbl) :: ddot - INTEGER :: nspin, nnr, ispin, j, k + INTEGER :: nspin, nnr, ispin, j, k, i ! vpot = vxc(rhoetr); vpot(r) <-- u(r) nnr = SIZE( rhoetr, 1 ) * SIZE( rhoetr, 2 ) * SIZE( rhoetr, 3 ) nspin = SIZE( rhoetr, 4 ) ! + IF( nnr /= nr3l * nr2l * nr1l ) THEN + DO ispin = 1, nspin + DO k = 1, SIZE( rhoetr, 3 ) + DO j = 1, SIZE( rhoetr, 2 ) + DO i = 1, SIZE( rhoetr, 1 ) + IF( i > nr1l .OR. j > nr2l .OR. k > nr3l ) THEN + rhoetr( i, j, k, ispin ) = 0.0d0 + grho ( i, j, k, :, ispin ) = 0.0d0 + END IF + END DO + END DO + END DO + END DO + END IF + ! CALL exch_corr_wrapper( nnr, nspin, grho(1,1,1,1,1), rhoetr(1,1,1,1), & sxc, vpot(1,1,1,1), v2xc(1,1,1,1,1) ) diff --git a/CPV/forces.f90 b/CPV/forces.f90 index ee2eb6f18..14b19d3de 100644 --- a/CPV/forces.f90 +++ b/CPV/forces.f90 @@ -102,7 +102,10 @@ ! ... declare modules USE spherical_harmonics USE ions_base, ONLY: na - USE pseudopotential, ONLY: l2ind, lm1x, nspnl, tl + USE pseudopotential, ONLY: nspnl + USE uspp_param, only: nh, lmaxkb + USE uspp, only: nhtol, nhtolm, indv + IMPLICIT NONE @@ -118,46 +121,40 @@ COMPLEX(dbl), ALLOCATABLE :: temp(:,:) REAL(dbl), ALLOCATABLE :: gwork(:,:) REAL(dbl) :: t1 - INTEGER :: igh, ll, is, isa, ig, l, m, ngw, nngw, iy + INTEGER :: igh, ll, is, isa, ig, l, m, ngw, nngw, iy, ih, iv ! end of declarations ! ---------------------------------------------- ngw = SIZE(df) nngw = 2*ngw - ALLOCATE(temp(ngw,2), gwork(ngw,(lm1x+1)**2)) + ALLOCATE(temp(ngw,2), gwork(ngw,(lmaxkb+1)**2)) - CALL ylmr2( (lm1x+1)**2, ngw, gx, hg, gwork ) + CALL ylmr2( (lmaxkb+1)**2, ngw, gx, hg, gwork ) - igh = 0 - iy = 0 - DO l = 0, lm1x - DO m = -l, l - iy = iy + 1 - IF(tl(l)) THEN - igh = igh + 1 - isa = 1 - DO is = 1, nspnl - ll = l2ind(l + 1,is) - IF(ll.GT.0) THEN - t1= - fio * wsg(igh,is) - CALL DGEMV('N', nngw, na(is), t1, eigr(1,isa), & - 2*SIZE(eigr,1), fnlo(isa,igh), 1, rzero, temp(1,1), 1) - t1= - fie * wsg(igh,is) - CALL DGEMV('N', nngw, na(is), t1, eigr(1,isa), & - 2*SIZE(eigr,1), fnle(isa,igh), 1, rzero, temp(1,2), 1) - CALL ZSCAL( nngw, cimgl(l), temp, 1) - DO ig=1,ngw - df(ig) = df(ig) + temp(ig,1) * wnl(ig,ll,is) * gwork(ig,iy) - END DO - DO ig=1,ngw - da(ig) = da(ig) + temp(ig,2) * wnl(ig,ll,is) * gwork(ig,iy) - END DO - END IF - isa=isa+na(is) - END DO - END IF + isa = 1 + DO is = 1, nspnl + DO ih = 1, nh( is ) + iv = indv ( ih, is ) + iy = nhtolm( ih, is ) + ll = nhtol ( ih, is ) + 1 + l = ll - 1 + igh = ih + t1= - fio * wsg(igh,is) + CALL DGEMV('N', nngw, na(is), t1, eigr(1,isa), & + 2*SIZE(eigr,1), fnlo(isa,igh), 1, rzero, temp(1,1), 1) + t1= - fie * wsg(igh,is) + CALL DGEMV('N', nngw, na(is), t1, eigr(1,isa), & + 2*SIZE(eigr,1), fnle(isa,igh), 1, rzero, temp(1,2), 1) + CALL ZSCAL( nngw, cimgl(l), temp, 1) + DO ig=1,ngw + df(ig) = df(ig) + temp(ig,1) * wnl(ig,iv,is) * gwork(ig,iy) + END DO + DO ig=1,ngw + da(ig) = da(ig) + temp(ig,2) * wnl(ig,iv,is) * gwork(ig,iy) + END DO END DO + isa=isa+na(is) END DO DEALLOCATE(temp, gwork) @@ -300,7 +297,9 @@ ! ... declare modules USE spherical_harmonics USE ions_base, ONLY: na - USE pseudopotential, ONLY: l2ind, lm1x, nspnl, tl + USE pseudopotential, ONLY: nspnl + USE uspp_param, only: nh, lmaxkb + USE uspp, only: nhtol, nhtolm, indv IMPLICIT NONE @@ -316,40 +315,36 @@ COMPLEX(dbl), ALLOCATABLE :: temp(:) REAL(dbl), ALLOCATABLE :: gwork(:,:) REAL(dbl) :: t1 - INTEGER :: igh, ll, is, isa, ig, l, m, iy + INTEGER :: igh, ll, is, isa, ig, l, m, iy, ih, iv INTEGER :: ngw, nngw ! end of declarations ngw = SIZE(df) nngw = 2*ngw - ALLOCATE( temp(ngw), gwork(ngw, (lm1x+1)**2) ) + ALLOCATE( temp(ngw), gwork(ngw, (lmaxkb+1)**2) ) - CALL ylmr2( (lm1x+1)**2, ngw, gx, hg, gwork ) + CALL ylmr2( (lmaxkb+1)**2, ngw, gx, hg, gwork ) - igh = 0 - iy = 0 - DO l = 0, lm1x - DO m = -l, l - iy = iy + 1 - IF(tl(l)) THEN - igh = igh + 1 - isa = 1 - DO is = 1, nspnl - ll = l2ind(l + 1,is) - IF(ll.GT.0) THEN - t1= - fi * wsg(igh,is) - CALL DGEMV('N', nngw, na(is), t1, eigr(1,isa), & - 2*SIZE(eigr,1), fnl(isa,igh), 1, rzero, temp(1), 1) - CALL ZSCAL(ngw, cimgl(l), temp(1), 1) - DO ig = 1, ngw - df(ig) = df(ig) + temp(ig) * wnl(ig,ll,is) * gwork(ig,iy) - END DO - END IF - isa = isa + na(is) - END DO - END IF + isa = 1 + DO is = 1, nspnl + DO ih = 1, nh( is ) + ! + iv = indv ( ih, is ) + iy = nhtolm( ih, is ) + ll = nhtol ( ih, is ) + 1 + l = ll - 1 + igh = ih + ! + t1= - fi * wsg(igh,is) + CALL DGEMV('N', nngw, na(is), t1, eigr(1,isa), & + 2*SIZE(eigr,1), fnl(isa,igh), 1, rzero, temp(1), 1) + CALL ZSCAL(ngw, cimgl(l), temp(1), 1) + DO ig = 1, ngw + df(ig) = df(ig) + temp(ig) * wnl(ig,iv,is) * gwork(ig,iy) + END DO END DO + isa = isa + na(is) END DO DEALLOCATE(temp, gwork) @@ -499,7 +494,9 @@ ! ... declare modules USE spherical_harmonics USE ions_base, ONLY: na - USE pseudopotential, ONLY: l2ind, lm1x, nspnl, tl + USE pseudopotential, ONLY: nspnl + USE uspp_param, only: nh, lmaxkb + USE uspp, only: nhtol, nhtolm, indv IMPLICIT NONE @@ -515,39 +512,35 @@ COMPLEX(dbl), ALLOCATABLE :: temp(:) REAL(dbl), ALLOCATABLE :: gwork(:,:) COMPLEX(dbl) fw - INTEGER :: igh, ll, is, isa, ig, l, m, ngw, iy + INTEGER :: igh, ll, is, isa, ig, l, m, ngw, iy, ih, iv ! end of declarations ngw = SIZE(df) - ALLOCATE( temp( SIZE(df) ), gwork( ngw, (lm1x+1)**2 ) ) + ALLOCATE( temp( SIZE(df) ), gwork( ngw, (lmaxkb+1)**2 ) ) - CALL ylmr2( (lm1x+1)**2, ngw, gx, hg, gwork ) + CALL ylmr2( (lmaxkb+1)**2, ngw, gx, hg, gwork ) - igh = 0 - iy = 0 - DO l = 0, lm1x - DO m = -l, l - iy = iy + 1 - IF(tl(l)) THEN - igh = igh + 1 - isa = 1 - DO is = 1, nspnl - ll = l2ind(l + 1,is) - IF(ll.GT.0) THEN -! ... TEMP(IG) = EIGR(IG,:) * FNLK(:,IGH) - fw = CMPLX( -fi * wsg(igh,is), 0.0d0) - CALL ZGEMV('N', ngw, na(is), fw, eigr(1,isa), & + isa = 1 + DO is = 1, nspnl + DO ih = 1, nh( is ) + ! + iv = indv ( ih, is ) + iy = nhtolm( ih, is ) + ll = nhtol ( ih, is ) + 1 + l = ll - 1 + igh = ih + ! + ! ... TEMP(IG) = EIGR(IG,:) * FNLK(:,IGH) + fw = CMPLX( -fi * wsg(igh,is), 0.0d0) + CALL ZGEMV('N', ngw, na(is), fw, eigr(1,isa), & size(eigr,1), fnlk(isa,igh), 1, czero, temp(1), 1) - CALL ZSCAL(ngw,cimgl(l),temp(1),1) - DO ig = 1, SIZE(df) - df(ig) = df(ig) + temp(ig) * wnl(ig,ll,is) * gwork(ig,iy) - END DO - END IF - isa = isa + na(is) - END DO - END IF + CALL ZSCAL(ngw,cimgl(l),temp(1),1) + DO ig = 1, SIZE(df) + df(ig) = df(ig) + temp(ig) * wnl(ig,iv,is) * gwork(ig,iy) + END DO END DO + isa = isa + na(is) END DO DEALLOCATE(temp, gwork) RETURN diff --git a/CPV/init.f90 b/CPV/init.f90 index 83162b107..df4101048 100644 --- a/CPV/init.f90 +++ b/CPV/init.f90 @@ -46,7 +46,9 @@ USE electrons_module, ONLY: electron_mass_init, band_init, n_emp USE electrons_base, ONLY: nspin, nupdwn USE reciprocal_vectors, ONLY: ngwt, gstart, gzero, ngm, ngmt, ngw, mill_l - USE pseudopotential, ONLY: formf, nsanl, ngh, pseudopotential_init + USE pseudopotential, ONLY: formf, nsanl, pseudopotential_init, & + pseudopotential_initval, & + pseudopotential_indexes USE ions_base, ONLY: nsp, na, nat USE ions_module, ONLY: atoms_init USE brillouin, ONLY: kpoints @@ -94,6 +96,8 @@ CALL atoms_init( atoms_m, atoms_0, atoms_p, taus, ind_srt, if_pos, atm, ht%hmat ) ! ... Allocate + Initialize pseudopotentials + CALL pseudopotential_indexes() + CALL pseudopotential_initval() CALL pseudopotential_init(ps, na, nsp, kp) s4 = cclock() diff --git a/CPV/input.f90 b/CPV/input.f90 index 6da20f181..fb33ef0e0 100644 --- a/CPV/input.f90 +++ b/CPV/input.f90 @@ -971,7 +971,7 @@ ! IF( program_name == 'FPMD' ) THEN - CALL pseudopotential_setup( ntyp, tpstab_inp, pstab_size_inp, ion_radius ) + CALL pseudopotential_setup( ntyp, tpstab_inp, pstab_size_inp ) ! o_diis_inp = .TRUE. oqnr_diis_inp = .TRUE. diff --git a/CPV/main.f90 b/CPV/main.f90 index a9a8ee2a2..b09401450 100644 --- a/CPV/main.f90 +++ b/CPV/main.f90 @@ -171,6 +171,10 @@ USE wave_init, ONLY: pw_atomic_init USE electrons_module, ONLY: electron_mass_init, band_init + ! + USE uspp_param , ONLY: nhm + ! + ! IMPLICIT NONE @@ -306,14 +310,17 @@ ALLOCATE( ei2( -nr2:nr2, nat ) ) ALLOCATE( ei3( -nr3:nr3, nat ) ) + ! initialize system geometry, cell and positions CALL init_geometry( ) + ! initialize variables and types ! CALL init0s(kp, ps, atoms_m, atoms_0, atoms_p, wfill, wempt, ht_m, ht_0 ) + lds_wfc = wfill%lds IF( force_pairing ) lds_wfc = 1 @@ -342,7 +349,7 @@ CALL electron_mass_init( alat, g, ngw ) ! ... initialize nonlocal pseudopotentials coefficients - CALL allocate_projector( fnl, nsanl, nbnd, ngh, kp%gamma_only) + CALL allocate_projector( fnl, nsanl, nbnd, nhm, kp%gamma_only) IF( t_diis ) THEN CALL allocate_diis( ngw, nbnd, kp%nkpt ) diff --git a/CPV/modules.f90 b/CPV/modules.f90 index 394fa7174..858e32717 100644 --- a/CPV/modules.f90 +++ b/CPV/modules.f90 @@ -90,18 +90,40 @@ end module core ! angular momentum (for r + + + SUBROUTINE compute_dvan() + use uspp, only: dvan, nhtolm, indv + use uspp_param, only: nhm, nh, dion + use ions_base, only: nsp + use atom, only: numeric + implicit none + integer :: is, iv, jv + real(dbl) :: fac + ! + if( allocated( dvan ) ) deallocate( dvan ) + allocate( dvan( nhm, nhm, nsp ) ) + dvan(:,:,:) =0.d0 + ! + do is = 1, nsp + if ( .not. numeric( is ) ) then + fac = 1.0d0 + else + ! fac converts ry to hartree + fac = 0.5d0 + end if + do iv=1,nh(is) + do jv=1,nh(is) + if ( nhtolm(iv,is) == nhtolm(jv,is) ) then + dvan( iv, jv, is ) = fac * dion( indv(iv,is), indv(jv,is), is ) + endif + end do + end do + end do + RETURN + END SUBROUTINE compute_dvan + + + +! ---------------------------------------------- + + + + SUBROUTINE pseudopotential_indexes() + + use parameters, only: lmaxx ! + use ions_base, only: nsp, & ! number of specie + na ! number of atoms for each specie + use cvan, only: ish ! + use uspp, only: nkb, & ! + nkbus ! + use core, only: nlcc_any ! + use uspp_param, only: nbeta, &! + lmaxkb, &! + lll, &! + nhm, &! + nh, &! + tvanp, &! + nqlc, &! + lmaxq ! + use uspp, only: nhtol, &! + nhtolm, &! + indv ! + use atom, only: nlcc ! + + + IMPLICIT NONE + + INTEGER :: is, iv, ind, il, lm + ! ------------------------------------------------------------------ + ! find number of beta functions per species, max dimensions, + ! total number of beta functions (all and Vanderbilt only) + ! ------------------------------------------------------------------ + lmaxkb = -1 + nhm = 0 + nkb = 0 + nkbus = 0 + nlcc_any = .false. + ! + do is = 1, nsp + ind = 0 + do iv = 1, nbeta(is) + lmaxkb = max( lmaxkb, lll( iv, is ) ) + ind = ind + 2 * lll( iv, is ) + 1 + end do + nh(is) = ind + nhm = max( nhm, nh(is) ) + ish(is)=nkb + nkb = nkb + na(is) * nh(is) + if( tvanp(is) ) nkbus = nkbus + na(is) * nh(is) + nlcc_any = nlcc_any .OR. nlcc(is) + end do + if (lmaxkb > lmaxx) call errore('nlinit ',' l > lmax ',lmaxkb) + lmaxq = 2*lmaxkb + 1 + ! + ! the following prevents an out-of-bound error: nqlc(is)=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 is=1,nsp + nqlc(is) = MIN ( nqlc(is), lmaxq ) + end do + if (nkb <= 0) call errore('nlinit ','not implemented ?',nkb) + + if( allocated( nhtol ) ) deallocate( nhtol ) + if( allocated( indv ) ) deallocate( indv ) + if( allocated( nhtolm ) ) deallocate( nhtolm ) + ! + allocate(nhtol(nhm,nsp)) + allocate(indv (nhm,nsp)) + allocate(nhtolm(nhm,nsp)) + + ! ------------------------------------------------------------------ + ! definition of indices nhtol, indv, nhtolm + ! ------------------------------------------------------------------ + ! + do is = 1, nsp + ind = 0 + do iv = 1, nbeta(is) + lm = lll(iv,is)**2 + do il = 1, 2*lll( iv, is ) + 1 + lm = lm + 1 + ind = ind + 1 + nhtolm( ind, is ) = lm + nhtol( ind, is ) = lll( iv, is ) + indv( ind, is ) = iv + end do + end do + end do + + + RETURN + END SUBROUTINE + + ! ---------------------------------------------- SUBROUTINE pseudopotential_init( ps, na, nsp, kp ) @@ -253,9 +333,11 @@ ! ... declare modules USE brillouin, ONLY: kpoints USE pseudo_types, ONLY: pseudo_ncpp, pseudo_upf + USE local_pseudo, ONLY: allocate_local_pseudo USE read_pseudo_module_fpmd, ONLY: ap USE gvecp, ONLY: ngm USE gvecw, ONLY: ngw + use uspp_param, only: lmaxkb, nhm IMPLICIT NONE @@ -272,7 +354,8 @@ tcc = .FALSE. WHERE ( ap(1:nsp)%tnlcc ) tcc = .TRUE. - CALL allocate_pseudo(ps, nsp, ngm, ngw, kp%nkpt, lnlx, ngh, tcc) + CALL allocate_local_pseudo( ngm, nsp ) + CALL allocate_pseudo(ps, nsp, ngm, ngw, kp%nkpt, nbetax, nhm, tcc) ps%ap => ap @@ -284,10 +367,11 @@ SUBROUTINE deallocate_pseudopotential USE splines, ONLY: kill_spline + USE uspp, ONLY: dvan INTEGER :: i, j - IF( ALLOCATED( wsginit ) ) DEALLOCATE( wsginit ) + IF( ALLOCATED( dvan ) ) DEALLOCATE( dvan ) IF( ALLOCATED( vps_sp ) ) THEN DO i = 1, size(vps_sp) CALL kill_spline(vps_sp(i),'a') @@ -349,7 +433,7 @@ ! ---------------------------------------------- ! ... declare modules - USE ions_base, ONLY: na, nsp + USE ions_base, ONLY: na, nsp, rcmax USE cell_module, ONLY: boxdimensions USE cell_base, ONLY: tpiba2 USE brillouin, ONLY: kpoints @@ -359,6 +443,10 @@ USE pseudo_base, ONLY: formfn USE pseudo_base, ONLY: compute_rhops, compute_rhocg USE reciprocal_vectors, ONLY: g, ngm + USE local_pseudo, ONLY: vps, dvps, rhops, drhops + use uspp_param, only: lmaxkb, nhm, nh + use uspp, only: dvan + USE constants, ONLY: pi IMPLICIT NONE @@ -368,9 +456,9 @@ TYPE (boxdimensions), INTENT(IN) :: ht ! ... declare other variables - INTEGER is, isc, ist, ig, igh, i + INTEGER is, isc, ist, ig, igh, i, l, iv REAL(dbl) omega, s1, s2, s3, s4 - REAL(dbl), ALLOCATABLE :: vps(:), dvps(:), vloc(:) + REAL(dbl), ALLOCATABLE :: vloc(:) ! end of declarations ! ---------------------------------------------- @@ -378,8 +466,8 @@ omega = ht%deth do is = 1, size(ps%wsg,2) - do igh = 1, size(ps%wsg,1) - ps%wsg(igh, is) = wsginit(igh, is) / omega + do igh = 1, nh( is ) + ps%wsg( igh, is) = 4.0d0 * ( 4.0d0 * pi ) ** 2 * dvan( igh, igh, is ) / omega end do end do @@ -398,15 +486,14 @@ ! ... handle local part DO is = 1, nsp - CALL compute_rhops( ps%rhops(:,is), ps%drhops(:,is), ps%ap(is)%zv, & - ps%ap(is)%raggio, g, omega, tpiba2, ngm, .true. ) + CALL compute_rhops( rhops(:,is), drhops(:,is), ps%ap(is)%zv, & + rcmax(is), g, omega, tpiba2, ngm, .true. ) IF( ps%ap(is)%tnlcc ) THEN IF(tpstab) THEN CALL corecortab_base(g, ps%rhoc1(:,is), ps%rhocp(:,is), & rhoc1_sp(is), rhocp_sp(is), xgtabmax, omega) ELSE - ! CALL corecor_base(ps%ap(is), g, ps%rhoc1(:,is), ps%rhocp(:,is), omega) CALL compute_rhocg( ps%rhoc1(:,is), ps%rhocp(:,is), ps%ap(is)%rw, & ps%ap(is)%rab, ps%ap(is)%rhoc, g, omega, tpiba2, ps%ap(is)%mesh, ngm, 1 ) @@ -415,34 +502,22 @@ ! ... numeric pseudopotential IF(tpstab) THEN - CALL formftab_base(g, ps%vps(:,is), ps%dvps(:,is), & + CALL formftab_base(g, vps(:,is), dvps(:,is), & vps_sp(is), dvps_sp(is), xgtabmax, omega ) ELSE ! ALLOCATE( vloc( ps%ap(is)%mesh ) ) ! vloc = ps%ap(is)%vloc * 2.0d0 - CALL formfn( ps%vps(:,is), ps%dvps(:,is), ps%ap(is)%rw, ps%ap(is)%rab, & - vloc, ps%ap(is)%zv, ps%ap(is)%raggio, g, omega, & + CALL formfn( vps(:,is), dvps(:,is), ps%ap(is)%rw, ps%ap(is)%rab, & + vloc, ps%ap(is)%zv, rcmax(is), g, omega, & tpiba2, 0.0d0, ps%ap(is)%mesh, ngm, .false., .true. ) - ps%dvps(:,is) = -ps%dvps(:,is) + dvps(:,is) = -dvps(:,is) ! DEALLOCATE( vloc ) END IF - ! DEBUG - ! IF( is == 2 ) THEN - ! DO i = 1, SIZE(ps%vps,1) - ! ps%dvps(1,is) = 0.0 - ! WRITE( stdout,fmt="(I5,2F14.6)" ) i,ps%vps(i,is),ps%dvps(i,is) - ! END DO - ! END IF - - ! WRITE( stdout,fmt="(/,'* FORMF TIMING: ',F18.8)" ) (s2-s1) - ! WRITE( stdout,fmt="(/,'* FORMF TIMING: ',F18.8)" ) (s3-s2) - ! WRITE( stdout,fmt="(/,'* FORMF TIMING: ',F18.8)" ) (s4-s3) - END DO @@ -459,7 +534,7 @@ SUBROUTINE build_pstab() - USE ions_base, ONLY: nsp + USE ions_base, ONLY: nsp, rcmax USE constants, ONLY: pi, fpi USE cell_base, ONLY: tpiba, tpiba2 USE splines, ONLY: init_spline, allocate_spline @@ -510,7 +585,7 @@ ALLOCATE( vloc( ap(is)%mesh ) ) vloc = ap(is)%vloc * 2.0d0 CALL formfn( vps_sp(is)%y, dvps_sp(is)%y, ap(is)%rw, ap(is)%rab, & - vloc, ap(is)%zv, ap(is)%raggio, xgtab, 1.0d0, & + vloc, ap(is)%zv, rcmax(is), xgtab, 1.0d0, & tpiba2, 0.0d0, ap(is)%mesh, pstab_size, .false., .true. ) dvps_sp(is)%y = -dvps_sp(is)%y DEALLOCATE( vloc ) @@ -533,23 +608,23 @@ NONLOCAL: IF( is <= nspnl ) THEN - DO l = 1, ap(is)%lnl + DO l = 1, ap(is)%nbeta CALL allocate_spline( wnl_sp(l,is), pstab_size, xgmin, xgmax ) CALL allocate_spline( wnla_sp(l,is), pstab_size, xgmin, xgmax ) END DO ALLOCATE( fintl( SIZE( xgtab ), SIZE( wnl_sp, 1) ) ) CALL nlin_base(ap(is), xgtab(:), fintl) - DO l = 1, ap(is)%lnl + DO l = 1, ap(is)%nbeta wnl_sp(l,is)%y = fintl(:,l) END DO CALL nlin_stress_base(ap(is), xgtab, fintl) - DO l = 1, ap(is)%lnl + DO l = 1, ap(is)%nbeta wnla_sp(l,is)%y = fintl(:,l) END DO DEALLOCATE(fintl) - DO l = 1, ap(is)%lnl + DO l = 1, ap(is)%nbeta CALL init_spline( wnl_sp(l,is) ) CALL init_spline( wnla_sp(l,is) ) END DO @@ -598,13 +673,13 @@ DO ik = 1, kp%nkpt IF( gamma_only ) THEN IF( tpstab ) THEN - CALL nlintab_base(g, wnl(:,:,is,ik), ap(is)%lnl, wnl_sp(:,is), xgtabmax) + CALL nlintab_base(g, wnl(:,:,is,ik), ap(is)%nbeta, wnl_sp(:,is), xgtabmax) ELSE CALL nlin_base(ap(is), g, wnl(:,:,is,ik)) END IF ELSE IF( tpstab ) THEN - CALL nlintab_base(gk_l(:,ik), wnl(:,:,is,ik), ap(is)%lnl, wnl_sp(:,is), xgtabmax) + CALL nlintab_base(gk_l(:,ik), wnl(:,:,is,ik), ap(is)%nbeta, wnl_sp(:,is), xgtabmax) ELSE CALL nlin_base(ap(is), gk_l(:,ik), wnl(:,:,is,ik)) END IF @@ -641,7 +716,7 @@ wnla = 0.0d0 DO IS = 1, NSPNL IF (tpstab) THEN - CALL nlintab_base(g, wnla(:,:,is), ap(is)%lnl, wnla_sp(:,is), xgtabmax) + CALL nlintab_base(g, wnla(:,:,is), ap(is)%nbeta, wnla_sp(:,is), xgtabmax) ELSE CALL nlin_stress_base(ap(is), g, wnla(:,:,is)) END IF @@ -692,6 +767,258 @@ RETURN END SUBROUTINE pseudo_wave_info + +! ---------------------------------------------- + +! calculation of array betagx(ig,iv,is) + +! ---------------------------------------------- + + + SUBROUTINE compute_betagx( tpre ) + ! + USE ions_base, ONLY: nsp + USE uspp_param, ONLY: nh, kkbeta, betar, nhm, nbeta + USE atom, ONLY: r, numeric, rab + USE uspp, ONLY: nhtol, indv + USE betax, only: refg, betagx, mmx, dbetagx + USE cvan, only: oldvan + USE qrl_mod, only: qrl, cmesh + ! + IMPLICIT NONE + ! + LOGICAL, INTENT(IN) :: tpre + ! + INTEGER :: is, iv, l, il, ltmp, i0, ir + REAL(dbl), ALLOCATABLE :: dfint(:), djl(:), fint(:), jl(:), jltmp(:) + REAL(dbl) :: xg, xrg + ! + IF( ALLOCATED( betagx ) ) DEALLOCATE( betagx ) + IF( ALLOCATED( dbetagx ) ) DEALLOCATE( dbetagx ) + ALLOCATE( betagx ( mmx, nhm, nsp ) ) + ALLOCATE( dbetagx( mmx, nhm, nsp ) ) + + ! + do is = 1, nsp + ! + if ( tpre ) then + allocate( dfint( kkbeta( is ) ) ) + allocate( djl ( kkbeta( is ) ) ) + allocate( jltmp( kkbeta( is ) ) ) + end if + allocate( fint ( kkbeta( is ) ) ) + allocate( jl ( kkbeta( is ) ) ) + ! + do iv = 1, nh(is) + ! + l = nhtol(iv,is) + 1 + ! + do il = 1, mmx + ! + xg = sqrt( refg * (il-1) ) + call sph_bes (kkbeta(is), r(1,is), xg, l-1, jl ) +! + if( tpre )then + ! + ltmp=l-1 + ! + ! r(i0) is the first point such that r(i0) >0 + ! + i0 = 1 + if ( r(1,is) < 1.0d-8 ) i0 = 2 + ! special case q=0 + if ( xg < 1.0d-8 ) then + if (l == 1) then + ! Note that dj_1/dx (x=0) = 1/3 + jltmp(:) = 1.0d0/3.d0 + else + jltmp(:) = 0.0d0 + end if + else + call sph_bes (kkbeta(is)+1-i0, r(i0,is), xg, ltmp-1, jltmp ) + end if + do ir = i0, kkbeta(is) + xrg = r(ir,is) * xg + djl(ir) = jltmp(ir) * xrg - l * jl(ir) + end do + if ( i0 == 2 ) djl(1) = djl(2) + ! + endif + ! + ! beta(ir)=r*beta(r) + ! + do ir = 1, kkbeta(is) + fint(ir) = r(ir,is) * betar( ir, indv(iv,is), is ) * jl(ir) + end do + if (oldvan(is)) then + call herman_skillman_int(kkbeta(is),cmesh(is),fint,betagx(il,iv,is)) + else + call simpson_cp90(kkbeta(is),fint,rab(1,is),betagx(il,iv,is)) + endif + ! + if(tpre) then + do ir = 1, kkbeta(is) + dfint(ir) = r(ir,is) * betar( ir, indv(iv,is), is ) * djl(ir) + end do + if (oldvan(is)) then + call herman_skillman_int(kkbeta(is),cmesh(is),dfint,dbetagx(il,iv,is)) + else + call simpson_cp90(kkbeta(is),dfint,rab(1,is),dbetagx(il,iv,is)) + end if + endif + ! + end do + end do +! + deallocate(jl) + deallocate(fint) + if (tpre) then + deallocate(jltmp) + deallocate(djl) + deallocate(dfint) + end if + ! + end do + + RETURN + END SUBROUTINE compute_betagx + + +! --------------------------------------------------------------- +! calculation of array qradx(igb,iv,jv,is) +! +! qradb(ig,l,k,is) = 4pi/omega int_0^r dr r^2 j_l(qr) q(r,l,k,is) +! +! --------------------------------------------------------------- + + + SUBROUTINE compute_qradx( tpre ) + ! + use io_global, only: stdout + USE ions_base, ONLY: nsp + USE uspp_param, ONLY: nh, kkbeta, betar, nhm, nqlc, qqq, nbrx, lmaxq, nbeta + USE atom, ONLY: r, numeric, rab + USE uspp, ONLY: nhtol, indv + USE betax, only: refg, qradx, mmx, dqradx + USE cvan, only: oldvan, ish, nvb + USE qrl_mod, only: qrl, cmesh + use gvecb, only: ngb + ! + IMPLICIT NONE + ! + LOGICAL, INTENT(IN) :: tpre + ! + INTEGER :: is, iv, l, il, ltmp, i0, ir, jv + REAL(dbl), ALLOCATABLE :: dfint(:), djl(:), fint(:), jl(:), jltmp(:) + REAL(dbl) :: xg, xrg + + IF( ALLOCATED( qradx ) ) DEALLOCATE( qradx ) + IF( ALLOCATED( dqradx ) ) DEALLOCATE( dqradx ) + ALLOCATE( qradx( mmx, nbrx, nbrx, lmaxq, nsp ) ) + ALLOCATE( dqradx( mmx, nbrx, nbrx, lmaxq, nsp ) ) + + DO is = 1, nvb + ! + IF ( tpre ) THEN + ALLOCATE( dfint( kkbeta(is) ) ) + ALLOCATE( djl ( kkbeta(is) ) ) + ALLOCATE( jltmp( kkbeta(is) ) ) + END IF + allocate( fint( kkbeta(is) ) ) + allocate( jl ( kkbeta(is) ) ) + ! + ! qqq and beta are now indexed and taken in the same order + ! as vanderbilts ppot-code prints them out + ! + WRITE( stdout,*) ' nlinit nh(is), ngb, is, kkbeta, lmaxq = ', & + & nh(is), ngb, is, kkbeta(is), nqlc(is) + ! + do l = 1, nqlc( is ) + ! + do il = 1, mmx + ! + xg = sqrt( refg * (il-1) ) + call sph_bes (kkbeta(is), r(1,is), xg, l-1, jl) + ! + if(tpre) then + ! + ltmp = l - 1 + ! + ! r(i0) is the first point such that r(i0) >0 + ! + i0 = 1 + if ( r(1,is) < 1.0d-8 ) i0 = 2 + ! special case q=0 + if ( xg < 1.0d-8 ) then + if (l == 1) then + ! Note that dj_1/dx (x=0) = 1/3 + jltmp(:) = 1.0d0/3.d0 + else + jltmp(:) = 0.0d0 + end if + else + call sph_bes (kkbeta(is)+1-i0, r(i0,is), xg, ltmp-1, jltmp ) + end if + do ir = i0, kkbeta(is) + xrg = r(ir,is) * xg + djl(ir) = jltmp(ir) * xrg - l * jl(ir) + end do + if (i0.eq.2) djl(1) = djl(2) + endif + ! + do iv = 1, nbeta(is) + do jv = iv, nbeta(is) + ! + ! note qrl(r)=r^2*q(r) + ! + do ir=1,kkbeta(is) + fint(ir)=qrl(ir,iv,jv,l,is)*jl(ir) + end do + if (oldvan(is)) then + call herman_skillman_int(kkbeta(is),cmesh(is),fint,qradx(il,iv,jv,l,is)) + else + call simpson_cp90(kkbeta(is),fint,rab(1,is),qradx(il,iv,jv,l,is)) + end if + ! + qradx(il,jv,iv,l,is)=qradx(il,iv,jv,l,is) + ! + if( tpre ) then + do ir = 1, kkbeta(is) + dfint(ir) = qrl(ir,iv,jv,l,is) * djl(ir) + end do + if ( oldvan(is) ) then + call herman_skillman_int(kkbeta(is),cmesh(is),dfint,dqradx(il,iv,jv,l,is)) + else + call simpson_cp90(kkbeta(is),dfint,rab(1,is),dqradx(il,iv,jv,l,is)) + end if + end if + ! + end do + end do + end do + end do + ! + WRITE( stdout,*) + WRITE( stdout,'(20x,a)') ' qqq ' + do iv=1,nbeta(is) + WRITE( stdout,'(8f9.4)') (qqq(iv,jv,is),jv=1,nbeta(is)) + end do + WRITE( stdout,*) + ! + deallocate(jl) + deallocate(fint) + if (tpre) then + deallocate(jltmp) + deallocate(djl) + deallocate(dfint) + end if + ! + end do + + RETURN + END SUBROUTINE compute_qradx + + ! ---------------------------------------------- END MODULE pseudopotential ! ---------------------------------------------- diff --git a/CPV/pseudotab_base.f90 b/CPV/pseudotab_base.f90 index d9f591eca..c134f5130 100644 --- a/CPV/pseudotab_base.f90 +++ b/CPV/pseudotab_base.f90 @@ -168,8 +168,4 @@ RETURN END SUBROUTINE corecortab_base -! ---------------------------------------------- -! ---------------------------------------------- - - END MODULE pseudotab_base - + END MODULE pseudotab_base diff --git a/CPV/read_pseudo.f90 b/CPV/read_pseudo.f90 index 4c1c9b3d2..1fd7dd5c1 100644 --- a/CPV/read_pseudo.f90 +++ b/CPV/read_pseudo.f90 @@ -37,7 +37,7 @@ TYPE (pseudo_ncpp), ALLOCATABLE, TARGET :: ap(:) TYPE (pseudo_upf), ALLOCATABLE, TARGET :: upf(:) - PUBLIC :: ap, upf, nspnl, l2ind, readpp + PUBLIC :: ap, upf, nspnl, readpp PUBLIC :: upf2internal, pseudo_filename, check_file_type !=----------------------------------------------------------------------------=! @@ -70,7 +70,6 @@ INTEGER FUNCTION check_file_type( is ) ! 0 file is unknow (guess: old CPV norm-conserving format) ! 1 file is *.vdb or *.van Vanderbilt US pseudopotential ! 2 file is *.RRKJ3 Andrea's US new code - ! 10 file is GIANNOZ (FPMD only) ! 11 file is NUMERIC (FPMD only) ! 12 file is ANALYTIC (FPMD only) ! 20 file is UPF @@ -100,9 +99,7 @@ INTEGER FUNCTION check_file_type( is ) READ ( pseudounit, *, iostat = ios, err = 300) dummy = ' ' READ ( pseudounit, *, iostat = ios, err = 300) dummy - IF( matches( "GIANNOZ", dummy ) ) THEN - info = 10 - ELSE IF( matches( "NUMERIC", dummy ) ) THEN + IF( matches( "NUMERIC", dummy ) ) THEN info = 11 ELSE IF( matches( "ANALYTIC", dummy ) ) THEN info = 12 @@ -134,28 +131,24 @@ SUBROUTINE check_types_order( ) INTEGER :: is, il LOGICAL :: tvanp ! + ! With Vanderbilt, only UPF are allowed + ! + IF( ANY( upf(1:nsp)%tvanp ) ) THEN + CALL errore( & + ' check_types_order ', ' vanderbilt pseudo, not yet implemened in FPMD ', 1 ) + END IF + ! ! non-local species must be ahead the local one, ! il = 0 DO is = 1, nsp - IF ( ap(is)%lnl == 0 ) THEN + IF ( ap(is)%nbeta == 0 ) THEN il = 1 ELSE IF ( il == 1 ) THEN CALL errore( & ' check_types_order ', ' Local pseudopotentials should follow non local ones ', 1 ) END IF END DO - ! - ! With Vanderbilt, only UPF are allowed - ! - IF( ANY( upf(1:nsp)%tvanp ) ) THEN - DO is = 1, nsp - IF ( ap(is)%pottyp /= 'UPF' ) THEN - CALL errore( & - ' check_types_order ', ' With vanderbilt pseudo, only UPF format is allowed ', 1 ) - END IF - END DO - END IF RETURN END SUBROUTINE check_types_order @@ -244,8 +237,7 @@ END FUNCTION calculate_dx CALL nullify_pseudo_upf( upf( is ) ) ! ap(is)%tnlcc = .FALSE. - ap(is)%raggio = 0.5d0 - ap(is)%lnl = 0 + ap(is)%nbeta = 0 upf(is)%tvanp = .FALSE. ! @@ -305,12 +297,6 @@ END FUNCTION calculate_dx call readAdC( is, pseudounit ) - ELSE IF( info == 10 ) THEN - - CALL read_head_pp( pseudounit, ap(is), error_msg, ierr) - CALL read_giannoz(pseudounit, ap(is), ierr) - CALL ncpp2internal ( ap(is), is, xc_type, ierr ) - ELSE IF( info == 11 ) THEN CALL read_head_pp( pseudounit, ap(is), error_msg, ierr) @@ -342,7 +328,7 @@ END FUNCTION calculate_dx IF( program_name == 'FPMD' ) THEN ! - IF( ap(is)%lnl > 0 ) nspnl = nspnl + 1 + IF( ap(is)%nbeta > 0 ) nspnl = nspnl + 1 IF( upf(is)%tvanp ) nvb = nvb + 1 IF( ionode ) THEN CALL ap_info( ap(is) ) @@ -402,27 +388,6 @@ END FUNCTION calculate_dx RETURN END SUBROUTINE readpp -!=----------------------------------------------------------------------------=! - - INTEGER FUNCTION l2ind( lw, is ) - - ! this function returns the index of the wanted channel lw - ! ( 1 = s, 2 = p, 3 = d, 4 = f ) within the non local pseudopotential - ! array WNL. The p component of the pseudopotential is stored, for the specie - ! is, in wnl(:,l2ind(2,is),is) - ! ---------------------------------------------- - - INTEGER, INTENT(IN) :: lw, is - INTEGER :: l - - l2ind = 0 - DO l = 1, ap(is)%lnl - IF ( ap(is)%indl(l) == lw ) l2ind = l - END DO - - RETURN - END FUNCTION l2ind - !=----------------------------------------------------------------------------=! SUBROUTINE analytic_to_numeric(ap) @@ -477,8 +442,8 @@ END FUNCTION calculate_dx ! ... Copy local component to a separate array ap%vloc(:) = ap%vnl(:,ap%lloc) - DO l = 1, ap%lnl - ll=ap%indl(l) ! find out the angular momentum (ll-1) of the component stored + DO l = 1, ap%nbeta + ll=ap%lll(l) + 1 ! find out the angular momentum (ll-1) of the component stored ! in position l ap%vrps(:,l) = ( ap%vnl(:,ll) - ap%vloc(:) ) * ap%rps(:,ll) END DO @@ -488,135 +453,6 @@ END FUNCTION calculate_dx !=----------------------------------------------------------------------------=! - SUBROUTINE read_giannoz(uni, ap, ierr) - USE constants, ONLY: pi - USE pseudo_types, ONLY: pseudo_ncpp - IMPLICIT NONE - TYPE (pseudo_ncpp), INTENT(INOUT) :: ap - INTEGER, INTENT(IN) :: uni - INTEGER, INTENT(OUT) :: ierr - REAL(dbl) :: chi( SIZE(ap%rps, 1), SIZE(ap%rps, 2) ) - REAL(dbl) :: vnl( SIZE(ap%vnl, 1), SIZE(ap%vnl, 2) ) - REAL(dbl) :: rho_core( SIZE(ap%rhoc, 1) ) - REAL(dbl) :: r, ra, rb, fac - REAL(dbl) :: oc( SIZE(ap%rps, 2) ) - REAL(dbl) :: enl( SIZE(ap%rps, 2) ) - REAL(dbl) :: zmesh, xmin, dx, etot - REAL(dbl) :: zval - INTEGER :: nn(SIZE(ap%rps, 2)), ll(SIZE(ap%rps, 2)) - INTEGER :: nwf, mesh, i, j, in1, in2, in3, in4, m - INTEGER :: lmax, nlc, nnl, lloc, l, il - LOGICAL :: nlcc - CHARACTER(len=80) :: dft - CHARACTER(len=4) :: atom - CHARACTER(len=2) :: el( SIZE(ap%rps, 2) ) - CHARACTER(len=80) :: ppinfo - CHARACTER(len=80) :: strdum - CHARACTER(len=2) :: sdum1, sdum2 - -! - ierr = 0 - - READ(uni,fmt='(a)') dft - READ(uni,fmt='(a4,f5.1,3i2,a2,l1,a2,i2,a)') & - atom, zval, lmax, nlc, nnl, sdum1, nlcc, sdum2, lloc, ppinfo - - ! WRITE( stdout,*) ' DEBUG ', atom, zval,lmax, nlc, nnl, nlcc, lloc, ppinfo - - IF( (lmax+1) > SIZE(ap%vnl, 2) ) THEN - ierr = 1 - RETURN - END IF - IF( (nlcc .AND. .NOT.ap%tnlcc) .OR. (.NOT.nlcc .AND. ap%tnlcc) ) THEN - ierr = 2 - RETURN - END IF - - READ(uni,fmt='(f8.2,f8.4,f10.6,2i6)') zmesh, xmin, dx, mesh, nwf - - IF( mesh > SIZE(ap%rps, 1) ) THEN - ierr = 3 - RETURN - END IF - IF( nwf > SIZE(ap%rps, 2) ) THEN - ierr = 4 - RETURN - END IF - - DO j = 0, lmax - READ(uni,fmt="(A16,i1)") strdum, l - READ(uni,'(4e16.8)') (vnl(i,j+1), i=1,mesh) - END DO - IF (nlcc) THEN - READ(uni,fmt='(4e16.8)') (rho_core(i), i=1,mesh) - END IF - DO j = 1, nwf - READ(uni,fmt="(A16,a2)") strdum,el(j) - READ(uni,fmt='(i5,f6.2)') ll(j),oc(j) - READ(uni,fmt='(4e16.8)') (chi(i,j), i=1,mesh) - END DO - - ap%zv = zval - ap%nchan = lmax+1 - ap%mesh = mesh - ap%rw = 0.0d0 - ap%vnl = 0.0d0 - ap%vrps = 0.0d0 - fac = 0.5d0 - - ! WRITE( stdout,*) ' DEBUG ', ap%lloc, ap%numeric, ap%lnl, ap%raggio, ap%zv - - DO i = 1, mesh - r = EXP(xmin+REAL(i-1)*dx)/zmesh - ap%rw(i) = r - DO j = 1, lmax+1 - ap%vnl(i,j) = vnl(i,j) * fac - END DO - END DO - IF( MINVAL( ap%rw(1:mesh) ) <= 0.0d0 ) THEN - ierr = 5 - RETURN - END IF - ap%dx = dx - ap%rab = ap%dx * ap%rw - ap%vloc(:) = ap%vnl(:,ap%lloc) - - ap%lrps(1:nwf) = ll(1:nwf) - ap%oc = 0.0d0 - ap%nrps = nwf - ap%mesh = mesh - ap%rps = 0.0d0 - fac = 1.0d0/SQRT(4.0d0*pi) - fac = 1.0d0 - DO i = 1, mesh - r = EXP(xmin+REAL(i-1)*dx)/zmesh - DO j = 1, nwf - ap%rps(i,j) = chi(i,j) * fac - END DO - END DO - - DO l = 1, ap%lnl - il=ap%indl(l) ! find out the angular momentum (il-1) of the component stored - ! in position l - DO i = 1, mesh - ap%vrps(i,l) = ( ap%vnl(i,il) - ap%vloc(i) ) * ap%rps(i,il) - END DO - END DO - - IF( nlcc ) THEN - ap%rhoc = 0.0d0 - DO i = 1, mesh - r = EXP(xmin+REAL(i-1)*dx)/zmesh - ap%rhoc(i) = rho_core(i) - END DO - END IF - - RETURN - END SUBROUTINE read_giannoz - -!=----------------------------------------------------------------------------=! - - SUBROUTINE ap_info( ap ) USE pseudo_types, ONLY: pseudo_ncpp USE io_global, ONLY: stdout @@ -625,16 +461,16 @@ END FUNCTION calculate_dx INTEGER :: in1, in2, in3, in4, m, il, ib, l, i WRITE( stdout, * ) - IF (ap%lnl > 0) THEN + IF (ap%nbeta > 0) THEN WRITE( stdout,10) ap%pottyp IF (ap%tmix) THEN WRITE( stdout,107) - WRITE( stdout,106) (ap%indl(l),l=1,ap%lnl) - WRITE( stdout,105) (ap%wgv(l),l=1,ap%lnl) + WRITE( stdout,106) (ap%lll(l),l=1,ap%nbeta) + WRITE( stdout,105) (ap%wgv(l),l=1,ap%nbeta) ELSE WRITE( stdout,50) ap%lloc END IF - WRITE( stdout,60) (ap%indl(l),l=1,ap%lnl) + WRITE( stdout,60) (ap%lll(l),l=1,ap%nbeta) ELSE ! ... A local pseudopotential has been read. WRITE( stdout,11) ap%pottyp @@ -652,9 +488,9 @@ END FUNCTION calculate_dx 60 FORMAT( 3X,'Non local components are : ',4I3) 11 FORMAT( 3X,'Type is ',A10,' and LOCAL. ') 12 FORMAT( 3X,'Using non local core corcorrections for this pseudo') - 20 FORMAT( 3X,'Pseudo charge : ',F8.3,', pseudo radius : ',F8.3) + 20 FORMAT( 3X,'Pseudo charge : ',F8.3) - WRITE( stdout,20) ap%zv, ap%raggio + WRITE( stdout,20) ap%zv IF( ap%pottyp /= 'ANALYTIC' ) THEN @@ -664,10 +500,10 @@ END FUNCTION calculate_dx in3=ap%mesh/2 in4=ap%mesh WRITE( stdout,132) - WRITE( stdout,120) in1,ap%rw(in1),ap%vloc(in1),(ap%vrps(in1,m),m=1,ap%lnl) - WRITE( stdout,120) in2,ap%rw(in2),ap%vloc(in2),(ap%vrps(in2,m),m=1,ap%lnl) - WRITE( stdout,120) in3,ap%rw(in3),ap%vloc(in3),(ap%vrps(in3,m),m=1,ap%lnl) - WRITE( stdout,120) in4,ap%rw(in4),ap%vloc(in4),(ap%vrps(in4,m),m=1,ap%lnl) + WRITE( stdout,120) in1,ap%rw(in1),ap%vloc(in1),(ap%vrps(in1,m),m=1,ap%nbeta) + WRITE( stdout,120) in2,ap%rw(in2),ap%vloc(in2),(ap%vrps(in2,m),m=1,ap%nbeta) + WRITE( stdout,120) in3,ap%rw(in3),ap%vloc(in3),(ap%vrps(in3,m),m=1,ap%nbeta) + WRITE( stdout,120) in4,ap%rw(in4),ap%vloc(in4),(ap%vrps(in4,m),m=1,ap%nbeta) 131 FORMAT(/, 3X,'Pseudopotentials Grid : Channels = ',I2,& ', Mesh = ',I5,/,30X,'dx = ',F16.14) 132 FORMAT( 3X,'point radius vloc ( vnl - vloc )') @@ -763,7 +599,7 @@ SUBROUTINE read_atomic_wf( iunit, ap, err_msg, ierr) ap%lrps = 0 ! this is for local pseudopotentials - IF( ap%lnl == 0 ) RETURN + IF( ap%nbeta == 0 ) RETURN READ(iunit,'(A80)',end=100) input_line CALL field_count(nf, input_line) @@ -829,7 +665,7 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr) err_msg = ' error while reading atomic numeric pseudo ' IF(ap%tmix) THEN - READ(iunit,*) (ap%wgv(l),l=1,ap%lnl) + READ(iunit,*) (ap%wgv(l),l=1,ap%nbeta) END IF READ(iunit,*,IOSTAT=ierr) ap%zv @@ -878,8 +714,8 @@ SUBROUTINE read_numeric_pp( iunit, ap, err_msg, ierr) CALL read_atomic_wf( iunit, ap, err_msg, ierr) IF( ierr /= 0 ) GO TO 110 - DO l = 1, ap%lnl - ll=ap%indl(l) + DO l = 1, ap%nbeta + ll=ap%lll(l) + 1 ap%vrps(:,l) = ( ap%vnl(:,ll) - ap%vloc(:) ) * ap%rps(:,ll) END DO @@ -1024,7 +860,7 @@ subroutine ncpp2internal ( ap, is, xc_type, ierr ) nlcc(is) = ap%tnlcc mesh(is) = ap%mesh nchi(is) = ap%nrps - nbeta(is) = ap%lnl + nbeta(is) = ap%nbeta kkbeta(is) = mesh(is) nqlc(is) = 0 ! upf%nqlc nqf (is) = 0 ! upf%nqf @@ -1037,9 +873,9 @@ subroutine ncpp2internal ( ap, is, xc_type, ierr ) lchi( 1 : ap%nrps, is ) = ap%lrps( 1 : ap%nrps ) chi ( 1 : ap%mesh, 1 : ap%nrps, is ) = ap%rps( 1 : ap%mesh, 1 : ap%nrps ) ! - betar( 1 : ap%mesh, 1 : ap%lnl, is ) = 2.0d0 * ap%vrps( 1 : ap%mesh, 1 : ap%lnl ) + betar( 1 : ap%mesh, 1 : ap%nbeta, is ) = 2.0d0 * ap%vrps( 1 : ap%mesh, 1 : ap%nbeta ) ! - lll ( 1 : ap%lnl, is ) = ap%indl( 1 : ap%lnl ) - 1 ! = upf%lll( 1:upf%nbeta ) + lll ( 1 : ap%nbeta, is ) = ap%lll( 1 : ap%nbeta ) ! = upf%lll( 1:upf%nbeta ) rinner(:,is) = 0.0d0 qqq(:,:,is) = 0.0d0 @@ -1109,19 +945,17 @@ SUBROUTINE upf2ncpp( upf, ap ) ap%tnlcc = upf%nlcc ! ap%zv = upf%zp - ap%lnl = upf%nbeta + ap%nbeta = upf%nbeta - ! angular momentum indl: S = 1, P = 2, ecc ... while lll: S = 0, P = 1, ecc ... - - ap%indl( 1:upf%nbeta ) = upf%lll( 1:upf%nbeta ) + 1 + ap%lll( 1:upf%nbeta ) = upf%lll( 1:upf%nbeta ) ! Calculate lloc ap%lloc = upf%nbeta + 1 which_lloc = 0 - DO l = 1, ap%lnl - which_lloc( ap%indl( l ) ) = 1 + DO l = 1, ap%nbeta + which_lloc( ap%lll( l ) + 1 ) = 1 END DO - DO l = 1, ap%lnl + 1 + DO l = 1, ap%nbeta + 1 IF( which_lloc( l ) == 0 ) ap%lloc = l END DO @@ -1140,11 +974,9 @@ SUBROUTINE upf2ncpp( upf, ap ) ap%lrps( 1:upf%nwfc ) = upf%lchi( 1:upf%nwfc ) ap%rps( 1:upf%mesh, 1:upf%nwfc ) = upf%chi( 1:upf%mesh, 1:upf%nwfc ) - DO l = 1, ap%lnl + DO l = 1, ap%nbeta - ! find out the angular momentum (il-1) of the component stored in position l - ! il = ap%indl(l) - ! vrps(i, il) = ( vnl(i, il) - vloc(i) ) * rps(i, il) + ! vrps(i, l) = ( vnl(i, l) - vloc(i) ) * rps(i, l) ap%vrps( 1:upf%mesh, l ) = upf%beta( 1:upf%mesh, l ) / 2.0d0 @@ -1175,13 +1007,14 @@ SUBROUTINE read_head_pp( iunit, ap, err_msg, ierr) ierr = 0 err_msg = ' error while reading header pseudo ' - ap%indl = 0 + ap%lll = 0 READ(iunit, *) ap%tnlcc, ap%tmix - READ(iunit, *) ap%pottyp, ap%lloc, ap%lnl, (ap%indl(l), l = 1, MIN(ap%lnl, SIZE(ap%indl)) ) + READ(iunit, *) ap%pottyp, ap%lloc, ap%nbeta, (ap%lll(l), l = 1, MIN(ap%nbeta, SIZE(ap%lll)) ) + ap%lll = ap%lll - 1 - IF( ap%lnl > SIZE(ap%indl) .OR. ap%lnl < 0 ) THEN + IF( ap%nbeta > SIZE(ap%lll) .OR. ap%nbeta < 0 ) THEN ierr = 1 - err_msg = 'LNL out of range' + err_msg = 'nbeta out of range' GO TO 110 END IF IF( ap%lloc < 0 .OR. ap%lloc > SIZE(ap%vnl,2) ) THEN @@ -1194,15 +1027,15 @@ SUBROUTINE read_head_pp( iunit, ap, err_msg, ierr) err_msg = 'tmix not implemented for pseudo ' // ap%pottyp GO TO 110 END IF - DO l = 2, ap%lnl - IF( ap%indl(l) <= ap%indl(l-1)) THEN + DO l = 2, ap%nbeta + IF( ap%lll(l) <= ap%lll(l-1)) THEN ierr = 5 err_msg =' NONLOCAL COMPONENTS MUST BE GIVEN IN ASCENDING ORDER' GO TO 110 END IF END DO - DO l = 1, ap%lnl - IF( ap%indl(l) == ap%lloc) THEN + DO l = 1, ap%nbeta + IF( ap%lll(l) + 1 == ap%lloc) THEN ierr = 6 err_msg = ' LLOC.EQ.L NON LOCAL!!' GO TO 110 diff --git a/CPV/rundiis.f90 b/CPV/rundiis.f90 index 9241cfda6..f0a888533 100644 --- a/CPV/rundiis.f90 +++ b/CPV/rundiis.f90 @@ -101,6 +101,7 @@ use grid_dimensions, only: nr1, nr2, nr3 USE reciprocal_vectors, ONLY: mill_l USE gvecp, ONLY: ngm + USE local_pseudo, ONLY: vps IMPLICIT NONE @@ -337,11 +338,11 @@ svar3_0 = delt * delt / pmss(1) IF(.NOT.kp%gamma_only) THEN CALL simupd_kp(ekinc,doions,kp,c0(:,:,:,1),cgrad(:,:,:,1),cdesc, svar1,svar2, & - svar3_0,edft%etot,fs(:,:,1),eigr,sfac,ps%vps,ps%wsg,ps%wnl,treset_diis, & + svar3_0,edft%etot,fs(:,:,1),eigr,sfac,vps,ps%wsg,ps%wnl,treset_diis, & istate,cnorm,eold,ndiis,nowv) ELSE CALL simupd(ekinc,doions,c0(:,:,:,1),cgrad(:,:,:,1),cdesc, svar1,svar2, & - svar3_0,edft%etot,fs(:,1,1),eigr,sfac,ps%vps,ps%wsg,ps%wnl, & + svar3_0,edft%etot,fs(:,1,1),eigr,sfac,vps,ps%wsg,ps%wnl, & treset_diis,istate,cnorm,eold,ndiis,nowv) END IF @@ -464,6 +465,7 @@ USE pseudo_projector, ONLY: projector USE atoms_type_module, ONLY: atoms_type USE charge_types, ONLY: charge_descriptor + USE local_pseudo, ONLY: vps IMPLICIT NONE @@ -606,12 +608,12 @@ IF(.NOT.kp%gamma_only) THEN CALL simupd_kp(ekinc(ispin), ddoions(ispin), kp, c0(:,:,:,ispin), & cgrad(:,:,:,ispin), cdesc, svar1, svar2, svar3_0, edft%etot, fi(:,:,ispin), & - eigr, sfac, ps%vps, ps%wsg, ps%wnl, ttreset_diis(ispin), istate, & + eigr, sfac, vps, ps%wsg, ps%wnl, ttreset_diis(ispin), istate, & cnorm, eold, ndiis, nowv) ELSE CALL simupd(ekinc(ispin), ddoions(ispin), c0(:,:,:,ispin), cgrad(:,:,:,ispin), cdesc, & svar1, svar2, svar3_0, edft%etot, fi(:,1,ispin), eigr, sfac, & - ps%vps, ps%wsg, ps%wnl, ttreset_diis(ispin), istate, cnorm, & + vps, ps%wsg, ps%wnl, ttreset_diis(ispin), istate, cnorm, & eold, ndiis, nowv) END IF CALL gram( ispin, c0(:,:,:,ispin), cdesc) diff --git a/CPV/smcp.f90 b/CPV/smcp.f90 index d6d39471b..1598749ad 100644 --- a/CPV/smcp.f90 +++ b/CPV/smcp.f90 @@ -82,14 +82,13 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out ) use cell_base, only: thdiag, tpiba2 use smooth_grid_dimensions, only: nnrsx, nr1s, nr2s, nr3s use smallbox_grid_dimensions, only: nnrb => nnrbx, nr1b, nr2b, nr3b - use pseu, only: vps, rhops + use local_pseudo, only: allocate_local_pseudo use io_global, ONLY: io_global_start, stdout, ionode use mp_global, ONLY: mp_global_start use mp, ONLY: mp_sum use para_mod use dener use derho - use dpseu use cdvan use stre use gvecw, only: ggp @@ -445,13 +444,11 @@ subroutine smdmain( tau, fion_out, etot_out, nat_out ) nlcc_any, smd = .TRUE. ) ! ! - allocate(rhops(ngs,nsp)) - allocate(vps(ngs,nsp)) + CALL allocate_local_pseudo( ngs, nsp ) + ! allocate(betae(ngw,nhsa)) allocate(deeq(nhm,nhm,nat,nspin)) allocate(dbec (nhsa,n,3,3)) - allocate(dvps(ngs,nsp)) - allocate(drhops(ngs,nsp)) allocate(drhog(ngm,nspin,3,3)) allocate(drhor(nnr,nspin,3,3)) allocate(drhovan(nhm*(nhm+1)/2,nat,nspin,3,3)) diff --git a/CPV/stress.f90 b/CPV/stress.f90 index 00556f60d..d6a9139f4 100644 --- a/CPV/stress.f90 +++ b/CPV/stress.f90 @@ -78,6 +78,7 @@ USE control_flags, ONLY: iprsta USE reciprocal_vectors, ONLY: gx USE gvecp, ONLY: ngm + USE local_pseudo, ONLY: dvps IMPLICIT NONE @@ -151,7 +152,7 @@ IF(timing) s5 = cclock() - CALL pseudo_stress(deps, edft%epseu, gagx_l, sfac, ps%dvps, rhoeg, box) + CALL pseudo_stress(deps, edft%epseu, gagx_l, sfac, dvps, rhoeg, box) IF(timing) s6 = cclock() @@ -239,8 +240,8 @@ ! ... declare modules - USE pseudopotential, ONLY: nlin_stress, ngh, & - l2ind,nspnl,nsanl, lnlx, ts,tp,td, tf, lm1x + USE pseudopotential, ONLY: nlin_stress, & + nspnl,nsanl USE ions_base, ONLY: nsp, na USE spherical_harmonics, ONLY: set_dmqm, set_fmrm, set_pmtm USE mp_global, ONLY: mpime, nproc @@ -249,6 +250,8 @@ USE cell_base, ONLY: tpiba2 USE control_flags, ONLY: force_pairing USE reciprocal_vectors, ONLY: gstart, gzero, g, gx + USE uspp_param, only: nh, lmaxkb, nbeta + USE uspp, only: nhtol, nhtolm, indv IMPLICIT NONE @@ -269,19 +272,19 @@ ! ... declare other variables INTEGER :: is, l, ll, me, al, be, s, k - INTEGER :: ir, kk, m, mm, isa, ighp, ig, iy - INTEGER :: igh, ia, ighd, ighf, in, i, iss, nx, ispin, nspin, ngw - INTEGER :: ispin_wfc, im(7) + INTEGER :: ir, kk, m, mm, isa, ig, iy, iv, iyy, ih, ihh + INTEGER :: ia, in, i, iss, nx, ispin, nspin, ngw + INTEGER :: ispin_wfc, mi(16), igh(0:3) REAL(dbl) xg,xrg,arg,wnd,wnd1,wnd2,temp,tt1,fac,tt2 - REAL(dbl) temp2, fg, gmod + REAL(dbl) temp2, fg, gmod, anm REAL(dbl) pm(3,3), pmtm(6,3,3) REAL(dbl) dm(6,5), dmqm(6,5,5) REAL(dbl) fm(3,3,3,7), fmrm(6,7,7) + REAL(dbl) facty(16) COMPLEX(dbl), ALLOCATABLE :: auxc(:,:) REAL(dbl), ALLOCATABLE :: wnla(:,:,:) REAL(dbl), ALLOCATABLE :: fnls(:,:) - REAL(dbl), ALLOCATABLE :: gwork(:) REAL(dbl), ALLOCATABLE :: gspha(:,:) REAL(dbl), ALLOCATABLE :: gwtmp(:) REAL(dbl), PARAMETER :: twothird = 2.0d0/3.0d0 @@ -306,14 +309,44 @@ ngw = cdesc%ngwl ! ... initialize array wnla - ALLOCATE(wnla(ngw, lnlx, nsp)) - CALL nlin_stress(wnla) + ALLOCATE( wnla( ngw, MAXVAL( nbeta( 1:nsp ) ), nsp) ) + CALL nlin_stress( wnla ) - ALLOCATE(gwork(ngw)) - ALLOCATE(gwtmp(ngw)) - ALLOCATE(gspha(ngw, (lm1x+1)**2 )) + ALLOCATE( gwtmp( ngw ) ) + ALLOCATE( gspha( ngw, (lmaxkb+1)**2 ) ) + + CALL ylmr2( (lmaxkb+1)**2, ngw, gx, g, gspha ) + + DO iy = 1, (lmaxkb+1)**2 + DO ig = gstart, ngw + gspha(ig,iy) = gspha(ig,iy) / (g(ig)*tpiba2) + END DO + END DO + + CALL set_pmtm( pm, pmtm ) + CALL set_dmqm( dm, dmqm ) + CALL set_fmrm( fm, fmrm ) + + mi( 1 ) = 1 + + mi( 2 ) = 2 ! im( 1 ) = 3 + mi( 3 ) = 3 ! im( 2 ) = 1 + mi( 4 ) = 1 ! im( 3 ) = 2 + + mi( 5 ) = 3 ! im( 1 ) = 5 + mi( 6 ) = 4 ! im( 2 ) = 3 + mi( 7 ) = 2 ! im( 3 ) = 1 + mi( 8 ) = 5 ! im( 4 ) = 2 + mi( 9 ) = 1 ! im( 5 ) = 4 + + mi( 10 ) = 4 ! im( 1 ) = 7 + mi( 11 ) = 5 ! im( 2 ) = 5 + mi( 12 ) = 3 ! im( 3 ) = 3 + mi( 13 ) = 6 ! im( 4 ) = 1 + mi( 14 ) = 2 ! im( 5 ) = 2 + mi( 15 ) = 7 ! im( 6 ) = 4 + mi( 16 ) = 1 ! im( 7 ) = 6 - CALL ylmr2( (lm1x+1)**2, ngw, gx, g, gspha ) SPIN_LOOP: DO ispin = 1, nspin @@ -324,357 +357,107 @@ IF( nx < 1 ) CYCLE SPIN_LOOP - ALLOCATE(fnls(nsanl,nx)) + iss = 1 + + SPECIES: DO is = 1, nspnl + + ALLOCATE(fnls(na(is),nx)) + ALLOCATE(auxc(ngw,na(is))) - igh = 0 - ! - IF (ts) THEN - ! - igh = igh + 1 - iy = 1 - ! - DO ig = gstart, ngw - gspha(ig,iy) = gspha(ig,iy) / (g(ig)*tpiba2) - END DO DO kk = 1, 6 - fnls = 0.0d0 - iss = 1 - gwork(1) = 0.0d0 - DO ig = gstart, ngw - gwork(ig) = gagx_l(kk,ig) * gspha(ig,iy) - END DO - DO is = 1, nspnl - ll = l2ind(1,is) - IF(ll.GT.0) THEN - - ALLOCATE(auxc(ngw,na(is))) - - gwtmp(1) = 0.0d0 - DO ig = gstart, ngw - gwtmp(ig) = gwork(ig) * (wnl(ig,ll,is)-wnla(ig,ll,is)) - END DO - - DO ia = 1, na(is) - auxc(1,ia) = CMPLX(0.0d0,0.0d0) - DO ig = gstart, ngw - auxc(ig,ia) = csign(0) * gwtmp(ig) * eigr(ig,ia+iss-1) - END DO - END DO - - CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, auxc(1,1), & - 2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, 0.0d0, fnls(iss,1), nsanl ) - - DEALLOCATE(auxc) - - END IF - iss = iss + na(is) - END DO - - CALL DSCAL(size(fnls),2.d0,fnls,1) - DO i = 1, nx - isa = 1 - DO is = 1, nspnl - tt1 = DDOT(na(is), fnl(1,ispin)%r(isa,igh,i), 1, fnls(isa,i), 1) - denl(kk) = denl(kk) + 2.0d0 * occ(i,1,ispin) * wsg(igh,is) * tt1 - isa = isa + na(is) - END DO - END DO - - END DO - END IF - - - ! - IF(tp) THEN - ! - ighp = igh - ! - im( 1 ) = 3 - im( 2 ) = 1 - im( 3 ) = 2 - ! - CALL set_pmtm( pm, pmtm ) - ! - DO m = 1, 3 - - iy = 1 + im( m ) - - DO ig = gstart, ngw - gspha(ig,iy) = gspha(ig,iy) / (g(ig)*tpiba2) - END DO - - DO kk=1,6 - - fnls = 0.0d0 - iss = 1 - gwork(1) = 0.0d0 - DO ig = gstart, ngw - gwork(ig)= gagx_l(kk,ig) * gspha(ig,iy) - END DO - - DO is=1,nspnl - - ll = l2ind(2,is) - IF(ll.GT.0) THEN - ALLOCATE(auxc(ngw,na(is))) - - gwtmp(1) = 0.0d0 - DO ig = gstart, ngw - gwtmp(ig) = gwork(ig) * ( 3.d0 * wnl(ig,ll,is) - wnla(ig,ll,is) ) - END DO - - DO ia = 1, na(is) - auxc(1,ia) = CMPLX(0.0d0,0.0d0) - DO ig = gstart, ngw - auxc(ig,ia) = csign(1) * gwtmp(ig) * eigr(ig,ia+iss-1) - END DO - END DO - - CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, & - auxc(1,1),2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, & - 0.0d0, fnls(iss,1), nsanl ) - - DEALLOCATE(auxc) - END IF - iss = iss + na(is) - END DO - - isa=0 - DO is=1,nspnl - temp = 2.d0 * wsg( ighp + im( m ), is) - DO ia=1,na(is) - isa=isa+1 - DO in=1,nx - IF(me.EQ.1) THEN - temp2=0.d0 - DO mm=1,3 - temp2=temp2 + pmtm(kk, m, mm )* & - fnl(1,ispin)%r(isa, ighp + im( mm ),in) - END DO - fnls(isa,in)= -temp2+2.d0*fnls(isa,in) - ELSE - fnls(isa,in)= 2.d0*fnls(isa,in) - END IF - tt1 = fnl(1,ispin)%r(isa, ighp + im( m ), in ) - fac = occ(in,1,ispin) * temp - tt2 = fnls(isa,in) - denl(kk) = denl(kk) + fac * tt1 * tt2 - END DO - END DO - END DO - - END DO - END DO - ! - igh = ighp + 3 - ! - END IF - - ! ... d-nonlocality - - IF(td) THEN - - ighd = igh - ! - im( 1 ) = 5 - im( 2 ) = 3 - im( 3 ) = 1 - im( 4 ) = 2 - im( 5 ) = 4 - - CALL set_dmqm( dm, dmqm ) - - DO m = 1, 5 - - iy = 4 + im( m ) - ! - DO ig = gstart, ngw - gspha(ig,iy) = gspha(ig,iy) / (g(ig)*tpiba2) - END DO + igh(0:3) = -1 - DO kk = 1, 6 + DO ih = 1, nh( is ) + + iy = nhtolm( ih, is ) + iv = indv ( ih, is ) + l = nhtol ( ih, is ) + anm = 2*l + 1 - fnls = 0.0d0 - iss = 1 - gwork(1) = 0.0d0 - DO ig=gstart,ngw - gwork(ig)= gagx_l(kk,ig) * gspha(ig,iy) - END DO + ! WRITE(6,*) 'DEBUG ih, iy, iv, l = ', ih, iy, iv, l - DO is = 1, nspnl - - ll = l2ind(3,is) - IF(ll.GT.0) THEN - ALLOCATE(auxc(ngw,na(is))) - - gwtmp(1) = 0.0d0 - DO ig = gstart, ngw - gwtmp(ig) = gwork(ig) * ( 5.d0 * wnl(ig,ll,is) - wnla(ig,ll,is) ) - END DO - - DO ig = gstart, ngw - gwtmp(ig) = gwtmp(ig) - 2.0d0/3.0d0 * dm( kk, m ) * wnl(ig,ll,is) - END DO - - DO ia= 1 , na(is) - auxc(1,ia) = CMPLX(0.0d0,0.0d0) - DO ig = gstart, ngw - auxc(ig,ia) = csign(2) * gwtmp(ig) * eigr(ig,ia+iss-1) - END DO - END DO - - CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, & - auxc(1,1), 2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, & - 0.0d0, fnls(iss,1), nsanl ) - - DEALLOCATE(auxc) - - END IF - iss = iss + na(is) - END DO - - isa=0 - DO is=1,nspnl - temp = 2.d0 * wsg(ighd+im(m),is) - DO ia=1,na(is) - isa=isa+1 - DO in=1,nx - IF(me.EQ.1) THEN - temp2=0.d0 - DO mm=1,5 - temp2=temp2 + dmqm(kk, m, mm )* & - fnl(1,ispin)%r(isa, ighd + im(mm),in) - END DO - fnls(isa,in)= -2.d0*temp2+2.d0*fnls(isa,in) - ELSE - fnls(isa,in)= 2.d0*fnls(isa,in) - END IF - tt1 = fnl(1,ispin)%r(isa,ighd+im(m),in) - fac = occ(in,1,ispin) * temp - tt2 = fnls(isa,in) - denl(kk) = denl(kk) + fac * tt1 * tt2 - END DO - END DO - END DO - - END DO - END DO - ! - igh = ighd + 5 - ! - END IF - - ! ... f-nonlocality - - IF(tf) THEN - - ighf = igh - - CALL set_fmrm(fm, fmrm) - - im( 1 ) = 7 - im( 2 ) = 5 - im( 3 ) = 3 - im( 4 ) = 1 - im( 5 ) = 2 - im( 6 ) = 4 - im( 7 ) = 6 - - DO m=1,7 - - iy = 9 + im( m ) - - DO ig = gstart, ngw - gspha(ig,iy) = gspha(ig,iy) / (g(ig)*tpiba2) - END DO - - DO kk=1,6 - - fnls = 0.0d0 - iss = 1 - gwork(1) = 0.0d0 + gwtmp(1) = 0.0d0 DO ig = gstart, ngw - gwork(ig)= gagx_l(kk,ig) * gspha(ig,iy) + gwtmp( ig ) = gagx_l(kk,ig) * gspha(ig,iy) * ( anm * wnl(ig,iv,is) - wnla(ig,iv,is) ) END DO - DO is=1,nspnl - - ll = l2ind(4,is) - - IF(ll > 0) THEN - ALLOCATE(auxc(ngw,na(is))) - - gwtmp(1) = 0.0d0 - DO ig = gstart, ngw - gwtmp(ig) = gwork(ig) * ( 7.d0 * wnl(ig,ll,is) - wnla(ig,ll,is) ) + IF( igh(l) < 0 ) igh(l) = ih + IF ( l == 1 ) THEN + ELSE IF ( l == 2 ) THEN + DO ig = gstart, ngw + gwtmp(ig) = gwtmp(ig) - 2.0d0/3.0d0 * dm( kk, mi( iy ) ) * wnl(ig,iv,is) + END DO + ELSE IF ( l == 3 ) THEN + al = alpha(kk) + be = beta(kk) + DO ig = gstart, ngw + fg = 0.0d0 + gmod = SQRT( g(ig) ) + DO s = 1, 3 + fg = fg + 3.0d0/5.0d0 * fm(be,s,s,mi(iy)) * gx(al,ig) / gmod END DO - - al = alpha(kk) - be = beta(kk) - DO ig = gstart, ngw - fg = 0.0d0 - gmod = SQRT( g(ig) ) - DO s = 1, 3 - fg = fg + 3.0d0/5.0d0 * fm(be,s,s,m) * gx(al,ig) / gmod - END DO - DO s = 1, 3 - fg = fg + 6.0d0/5.0d0 * fm(be,s,al,m) * gx(s,ig) / gmod - END DO - gwtmp(ig) = gwtmp(ig) - fg * wnl(ig,ll,is) + DO s = 1, 3 + fg = fg + 6.0d0/5.0d0 * fm(be,s,al,mi(iy)) * gx(s,ig) / gmod END DO - - DO ia= 1 , na(is) - auxc(1,ia) = CMPLX(0.0d0,0.0d0) - DO ig = gstart, ngw - auxc(ig,ia) = csign(3) * gwtmp(ig) * eigr(ig,ia+iss-1) - END DO - END DO - - CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, & - auxc(1,1), 2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, & - 0.0d0, fnls(iss,1), nsanl) - - DEALLOCATE(auxc) - + gwtmp(ig) = gwtmp(ig) - fg * wnl(ig,iv,is) + END DO + END IF + DO ihh = igh(l), igh(l) + 2*l + iyy = nhtolm( ihh, is ) + IF ( l == 0 ) THEN + facty( ihh ) = 0.0d0 + ELSE IF( l == 1 ) THEN + facty( ihh ) = pmtm(kk, mi( iy ), mi( iyy ) ) + ELSE IF( l == 2 ) THEN + facty( ihh ) = dmqm(kk, mi( iy ), mi( iyy ) ) + ELSE IF( l == 3 ) THEN + facty( ihh ) = fmrm(kk, mi( iy ), mi( iyy ) ) END IF - - iss = iss + na(is) END DO + ! + DO ia = 1, na(is) + auxc(1,ia) = CMPLX(0.0d0,0.0d0) + DO ig = gstart, ngw + auxc(ig,ia) = csign(l) * gwtmp(ig) * eigr(ig,ia+iss-1) + END DO + END DO + + CALL DGEMM( 'T', 'N', na(is), nx, 2*ngw, 1.0d0, auxc(1,1), & + 2*ngw, c0(1,1,1,ispin_wfc), 2 * cdesc%ldg, 0.0d0, fnls(1,1), na(is) ) - isa=0 - DO is=1,nspnl - temp = 2.d0 * wsg( ighf + im(m), is ) - DO ia=1,na(is) - isa=isa+1 - DO in=1,nx - IF(me == 1) THEN - temp2=0.d0 - DO mm=1,7 - temp2=temp2 + fmrm(kk,m,mm) * fnl(1,ispin)%r(isa,ighf+im(mm),in) - END DO - fnls(isa,in)= -3.d0*temp2+2.d0*fnls(isa,in) - ELSE - fnls(isa,in)= 2.d0*fnls(isa,in) - END IF - tt1 = fnl(1,ispin)%r(isa,ighf+im(m),in) - fac = occ(in,1,ispin) * temp - tt2 = fnls(isa,in) - denl(kk) = denl(kk) + fac * tt1 * tt2 - END DO + DO in = 1, nx + ! + fac = 2.0d0 * occ( in, 1, ispin ) * wsg( ih, is) + ! + DO ia = 1, na(is) + isa = iss + ia - 1 + temp2 = 0.d0 + IF( me == 1 ) THEN + DO ihh = igh(l), igh(l) + 2*l + temp2 = temp2 + facty( ihh ) * fnl(1,ispin)%r( isa, ihh, in ) + END DO + END IF + tt1 = fnl(1,ispin)%r(isa, ih, in ) + tt2 = - l * temp2 + 2.d0 * fnls( ia, in ) + denl(kk) = denl(kk) + fac * tt1 * tt2 END DO END DO END DO END DO - END IF + ! + DEALLOCATE(auxc) + DEALLOCATE(fnls) + + iss = iss + na(is) + ! + END DO SPECIES - DEALLOCATE(fnls) END DO SPIN_LOOP - DEALLOCATE(gwork) DEALLOCATE(gwtmp) DEALLOCATE(gspha) DEALLOCATE(wnla) @@ -827,7 +610,7 @@ SUBROUTINE stress_har(deht, ehr, sfac, ps, rhoeg, gagx_l, box ) - use ions_base, only: nsp + use ions_base, only: nsp, rcmax USE cell_module, only: boxdimensions use mp_global, ONLY: mpime, nproc USE constants, ONLY: fpi @@ -835,6 +618,7 @@ USE cp_types, ONLY: pseudo, pseudo_ncpp USE reciprocal_vectors, ONLY: gstart, g USE gvecp, ONLY: ngm + USE local_pseudo, ONLY: rhops IMPLICIT NONE @@ -874,9 +658,9 @@ RHOP = (0.D0,0.D0) RHOPR= (0.D0,0.D0) DO IS = 1, NSP - RHOP = RHOP + sfac( IG, is ) * ps%RHOPS(IG,is) + RHOP = RHOP + sfac( IG, is ) * RHOPS(IG,is) ! RHOPR = RHOPR + sfac( IG, is ) * ps%DRHOPS(IG,is) - RHOPR = RHOPR + sfac( IG, is ) * ps%RHOPS(IG,is) * ps%ap(is)%RAGGIO**2 * 0.5D0 + RHOPR = RHOPR + sfac( IG, is ) * RHOPS(IG,is) * rcmax(is)**2 * 0.5D0 END DO HGM1 = 1.D0 / g(IG) / TPIBA2 RHET = 0.0_dbl diff --git a/CPV/vofrho2.f90 b/CPV/vofrho2.f90 index 03ae3647c..b353dfa28 100644 --- a/CPV/vofrho2.f90 +++ b/CPV/vofrho2.f90 @@ -42,13 +42,12 @@ use energies, only: etot, eself, enl, ekin, epseu, esr, eht, exc, & & atot, egrand, entropy - use pseu + use local_pseudo, only: vps, rhops use core use gvecb ! use dener use derho - use dpseu use mp, only: mp_sum ! implicit none diff --git a/CPV/wf.f90 b/CPV/wf.f90 index 646b31444..878249bc0 100644 --- a/CPV/wf.f90 +++ b/CPV/wf.f90 @@ -4104,7 +4104,6 @@ subroutine tric_wts(rp1,rp2,rp3,alat,wts) nr1sx, nr2sx, nr3sx, nnrsx use electrons_base, only: nx => nbspx, n => nbsp, nspin, f, ispin => fspin use constants, only: pi, fpi - use pseu use wfparm, only : iwf ! use cdvan diff --git a/Modules/pseudo_types.f90 b/Modules/pseudo_types.f90 index bc3050f3c..f51f8609f 100644 --- a/Modules/pseudo_types.f90 +++ b/Modules/pseudo_types.f90 @@ -85,12 +85,11 @@ LOGICAL :: tnlcc INTEGER :: igau INTEGER :: lloc - INTEGER :: lnl - INTEGER :: indl(ndmx) + INTEGER :: nbeta + INTEGER :: lll(cp_lmax) INTEGER :: nchan INTEGER :: mesh REAL(dbl) :: zv - REAL(dbl) :: raggio REAL(dbl) :: dx ! r(i) = cost * EXP( xmin + dx * (i-1) ) REAL(dbl) :: rab(ndmx) REAL(dbl) :: rw(ndmx)