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

- pseudopotential base modules uspp_param and uspp now used 

also in FPMD, for norm-conserving pseudo (like in CP)
- Few clean-ups in pseudopotential parameters initialization


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1880 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
cavazzon 2005-05-18 21:01:05 +00:00
parent ff0b81c104
commit d3aabc5972
31 changed files with 988 additions and 1374 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
CPV/bessel.f90
View File

@ -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) ) )

4
CPV/cg_sub.f90
View File

@ -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

7
CPV/cp_fpmd.f90
View File

@ -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')

18
CPV/cplib.f90
View File

@ -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
!

10
CPV/cpr.f90
View File

@ -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 ) )

332
CPV/cprsub.f90
View File

@ -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

29
CPV/cptypes.f90
View File

@ -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 )

6
CPV/dealloc.f90
View File

@ -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()

26
CPV/diis.f90
View File

@ -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

12
CPV/electrons.f90
View File

@ -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

5
CPV/emptystates.f90
View File

@ -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) ) )

17
CPV/exch_corr.f90
View File

@ -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) )

161
CPV/forces.f90
View File

@ -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

6
CPV/init.f90
View File

@ -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()

2
CPV/input.f90
View File

@ -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.

9
CPV/main.f90
View File

@ -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 )

41
CPV/modules.f90
View File

@ -90,18 +90,40 @@ end module core
! angular momentum (for r<rinner)
! vloc_at local potential for each atom
module pseu
module local_pseudo
implicit none
save
!
! rhops = ionic pseudocharges (for Ewald term)
! vps = local pseudopotential in G space for each species
!
real(kind=8), allocatable:: rhops(:,:), vps(:,:)
!
! drhops = derivative of rhops respect to G^2
! dvps = derivative of vps respect to G^2
!
real(kind=8),allocatable:: dvps(:,:), drhops(:,:)
!
contains
subroutine deallocate_pseu
!
subroutine allocate_local_pseudo( ng, nsp )
integer, intent(in) :: ng, nsp
call deallocate_local_pseudo()
ALLOCATE( rhops( ng, nsp ) )
ALLOCATE( vps( ng, nsp ) )
ALLOCATE( drhops( ng, nsp ) )
ALLOCATE( dvps( ng, nsp ) )
end subroutine
!
subroutine deallocate_local_pseudo
IF( ALLOCATED( rhops ) ) DEALLOCATE( rhops )
IF( ALLOCATED( vps ) ) DEALLOCATE( vps )
end subroutine deallocate_pseu
end module pseu
IF( ALLOCATED( dvps ) ) DEALLOCATE( dvps )
IF( ALLOCATED( drhops ) ) DEALLOCATE( drhops )
end subroutine
!
end module local_pseudo
module qgb_mod
implicit none
@ -153,17 +175,6 @@ contains
end subroutine deallocate_dqgb_mod
end module dqgb_mod
module dpseu
implicit none
save
real(kind=8),allocatable:: dvps(:,:), drhops(:,:)
contains
subroutine deallocate_dpseu
IF( ALLOCATED( dvps ) ) DEALLOCATE( dvps )
IF( ALLOCATED( drhops ) ) DEALLOCATE( drhops )
end subroutine deallocate_dpseu
end module dpseu
module cdvan
implicit none
save

146
CPV/nl.f90
View File

@ -90,7 +90,6 @@
CALL nl_projector_m(c0, cdesc, atoms, kp, fnl, wsg, wnl, eigr)
nlrh_m = nl_energy_m(fnl, atoms, occ, kp, wsg)
! WRITE( stdout,*) 'DEBUG nlrh ', SUM( fnl(1,1)%r ), SUM( wsg ), SUM( wnl )
IF( tforce ) THEN
DO ispin = 1, cdesc%nspin
ispin_wfc = ispin
@ -99,7 +98,6 @@
kp, fnl(:,ispin), wsg, wnl, eigr)
END DO
END IF
! .. WRITE( stdout,*) 'DEBUG NLRH:',atoms%for(:,1:atoms%nat)
RETURN
END FUNCTION nlrh_m
@ -285,7 +283,7 @@
! END manual
! ... include modules
USE pseudopotential, ONLY: nspnl, nsanl, ngh
USE pseudopotential, ONLY: nspnl, nsanl
USE brillouin, ONLY: kpoints
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector, allocate_projector, deallocate_projector
@ -293,6 +291,7 @@
USE reciprocal_space_mesh, ONLY: gkx_l, gk_l
USE control_flags, ONLY: gamma_only
USE reciprocal_vectors, ONLY: gx, g
USE uspp_param, ONLY: nhm
IMPLICIT NONE
@ -321,14 +320,14 @@
! ... compute nonlocal contribution to forces on ions
KAPPA: DO ik = 1, cdesc%nkl
CARTE: DO k = 1, 3 ! x,y,z directions
CALL allocate_projector( dfnl, nsanl, cdesc%nbl( ispin ), ngh, cdesc%gamma )
CALL allocate_projector( dfnl, nsanl, cdesc%nbl( ispin ), nhm, cdesc%gamma )
IF( gamma_only ) THEN
CALL nlsm2_s( ispin, wnl(:,:,:,ik), atoms, eigr, c0(:,:,ik), cdesc, g, gx, dfnl, k)
ELSE
CALL nlsm2_s( ispin, wnl(:,:,:,ik), atoms, eigr, c0(:,:,ik), cdesc, &
gk_l(:,ik), gkx_l(:,:,ik), dfnl, k)
END IF
CHANN: DO igh = 1, ngh
CHANN: DO igh = 1, nhm
BANDE: DO ib = 1, cdesc%nbl( ispin )
isa=0
SPECS: DO is = 1, nspnl
@ -381,12 +380,14 @@
! ... declare modules
USE pseudopotential, ONLY: l2ind, nspnl, tl, lm1x
USE pseudopotential, ONLY: nspnl
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE mp
USE mp_global, ONLY: nproc, mpime, group
USE atoms_type_module, ONLY: atoms_type
USE uspp_param, only: nh, lmaxkb
USE uspp, only: nhtol, nhtolm, indv
IMPLICIT NONE
@ -402,7 +403,7 @@
! ... declare other variables
INTEGER :: is, igh, isa, iss, ia, ig, nb, iy
INTEGER :: l, ll, m, ngw, lda, ldw, ldf
INTEGER :: l, ll, m, ngw, lda, ldw, ldf, ih, iv
REAL(dbl), ALLOCATABLE :: gwork(:,:)
COMPLEX(dbl), ALLOCATABLE :: gxtmp(:)
COMPLEX(dbl), ALLOCATABLE :: auxc(:,:)
@ -429,9 +430,8 @@
fnl%c = 0.0d0
END IF
! WRITE( stdout,fmt="('DEBUG nlsm1 ',10I6)" ) ngw, nb, lda, ldw, ldf
ALLOCATE( gwork( ngw, (lm1x+1)**2 ), gxtmp(ngw) )
ALLOCATE( gwork( ngw, (lmaxkb+1)**2 ), gxtmp(ngw) )
! ... angular momentum l = 0
! ... orbital: s
@ -440,51 +440,44 @@
! ... angular momentum l = 2
! ... orbitals: d_{z^2}, d_{x^2-y^2}, d_{xy}, d_{yz}, d_{zx}
CALL ylmr2( (lm1x+1)**2, ngw, gx, g2, gwork )
CALL ylmr2( (lmaxkb+1)**2, ngw, gx, g2, gwork )
igh = 0
iy = 0
DO l = 0, lm1x
DO m = -l, l
iy = iy + 1
IF(tl(l)) THEN
igh = igh + 1
iss = 1
DO is = 1, nspnl
ll = l2ind( l + 1, is )
IF(ll.gt.0) THEN
ALLOCATE( auxc( ngw, atoms%na(is) ) )
gxtmp(1:ngw) = csign(l) * wnl(1:ngw,ll,is) * gwork(1:ngw, iy )
IF( cdesc%gamma .AND. cdesc%gzero ) gxtmp(1) = gxtmp(1) * 0.5d0
DO ia = 1, atoms%na(is)
auxc(1:ngw,ia) = gxtmp(1:ngw) * eigr(1:ngw,iss+ia-1)
END DO
IF ( cdesc%gamma ) THEN
CALL DGEMM( 'T', 'N', atoms%na(is), nb, 2*ngw, 1.0d0, &
auxc(1,1), lda, c(1,1), ldw, 0.0d0, fnl%r(iss,igh,1), ldf)
ELSE
CALL ZGEMM( 'C', 'N', atoms%na(is), nb, ngw, one, auxc(1,1), lda, &
c(1,1), ldw, zero, fnl%c(iss,igh,1), ldf )
END IF
DEALLOCATE(auxc)
! WRITE( stdout,*) 'DEBUG nlsm_s', tl(l), l, SUM( fnl%r )
END IF
iss = iss + atoms%na(is)
END DO
iss = 1
DO is = 1, nspnl
ALLOCATE( auxc( ngw, atoms%na(is) ) )
DO ih = 1, nh( is )
iy = nhtolm( ih, is )
iv = indv ( ih, is )
ll = nhtol ( ih, is ) + 1
l = ll - 1
igh = ih
gxtmp(1:ngw) = csign(l) * wnl(1:ngw, iv, is) * gwork(1:ngw, iy )
! WRITE(6,* ) 'debug is, igh, ll, iy, iv = ', is, igh, ll, iy, iv ! debug
IF( cdesc%gamma .AND. cdesc%gzero ) gxtmp(1) = gxtmp(1) * 0.5d0
DO ia = 1, atoms%na(is)
auxc(1:ngw,ia) = gxtmp(1:ngw) * eigr(1:ngw,iss+ia-1)
END DO
IF ( cdesc%gamma ) THEN
CALL DGEMM( 'T', 'N', atoms%na(is), nb, 2*ngw, 1.0d0, &
auxc(1,1), lda, c(1,1), ldw, 0.0d0, fnl%r(iss,igh,1), ldf)
ELSE
CALL ZGEMM( 'C', 'N', atoms%na(is), nb, ngw, one, auxc(1,1), lda, &
c(1,1), ldw, zero, fnl%c(iss,igh,1), ldf )
END IF
END DO
iss = iss + atoms%na(is)
DEALLOCATE(auxc)
END DO
DEALLOCATE(gwork, gxtmp)
! ... since G vectors only span half space, multiply results by two
IF ( cdesc%gamma ) THEN
CALL DSCAL( size( fnl%r ), 2.0d0, fnl%r(1,1,1), 1 )
! WRITE( stdout,*) ' DEBUG nlsm1: ', SIZE(fnl%r,1), SIZE(fnl%r,2), SIZE(fnl%r,3) ! DEBUG
CALL mp_sum( fnl%r, group )
ELSE
CALL mp_sum( fnl%c, group )
END IF
RETURN
END SUBROUTINE nlsm1_s
@ -548,11 +541,13 @@
! END manual
! ... declare modules
USE pseudopotential, ONLY: l2ind, lm1x, nspnl, tl
USE pseudopotential, ONLY: nspnl
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector
USE atoms_type_module, ONLY: atoms_type
USE cell_base, ONLY: tpiba
USE uspp_param, only: nh, lmaxkb
USE uspp, only: nhtol, nhtolm, indv
IMPLICIT NONE
@ -569,7 +564,7 @@
! ... declare other variables
REAL(dbl), ALLOCATABLE :: gwork(:,:)
INTEGER :: is, ia, igh, isa, ig, iss, ll, l, m, ngw, nb, lda, ldw, ldf
INTEGER :: iy
INTEGER :: iy, ih, iv
COMPLEX(dbl), ALLOCATABLE :: auxc(:,:), gxtmp(:)
COMPLEX(dbl), PARAMETER :: ONE = (1.0d0,0.0d0), ZERO = (0.0d0,0.0d0)
! ... (-i) * i^l
@ -593,45 +588,42 @@
dfnl%c = 0.0d0
END IF
ALLOCATE(gwork(ngw, (lm1x+1)**2 ), gxtmp(ngw))
ALLOCATE(gwork(ngw, (lmaxkb+1)**2 ), gxtmp(ngw))
CALL ylmr2( (lm1x+1)**2, ngw, gx, g2, gwork )
CALL ylmr2( (lmaxkb+1)**2, ngw, gx, g2, gwork )
!
DO iy = 1, (lmaxkb+1)**2
gwork(1:ngw,iy) = tpiba * gx(kk,1:ngw) * gwork(1:ngw,iy)
END DO
igh = 0
iy = 0
ANGMO: DO l = 0, lm1x
MAGNE: DO m = -l, l
iy = iy + 1
IF(tl(l)) THEN
igh = igh + 1
iss = 1
gwork(1:ngw,iy) = tpiba * gx(kk,1:ngw) * gwork(1:ngw,iy)
SPECS: DO is = 1, nspnl
ll = l2ind(l + 1,is)
IF(ll.gt.0) THEN
ALLOCATE(auxc(ngw,atoms%na(is)))
gxtmp(1:ngw) = csign(l) * wnl(1:ngw,ll,is) * gwork(1:ngw,iy)
DO ia = 1, atoms%na(is)
auxc(1:ngw,ia) = gxtmp(1:ngw) * eigr(1:ngw,iss + ia - 1)
END DO
IF( cdesc%gamma ) THEN
CALL DGEMM('T', 'N', atoms%na(is), nb, 2*ngw, 1.0d0, auxc(1,1), lda, &
c(1,1), ldw, 0.0d0, dfnl%r(iss,igh,1), ldf)
ELSE
CALL ZGEMM('C', 'N', atoms%na(is), nb, ngw, one, auxc(1,1), lda, &
c(1,1), ldw, zero, dfnl%c(iss,igh,1), ldf)
END IF
DEALLOCATE(auxc)
END IF
iss = iss + atoms%na(is)
END DO SPECS
END IF
END DO MAGNE
END DO ANGMO
iss = 1
SPECS: DO is = 1, nspnl
ALLOCATE(auxc(ngw,atoms%na(is)))
LM: DO ih = 1, nh( is )
iv = indv ( ih, is )
iy = nhtolm( ih, is )
ll = nhtol ( ih, is ) + 1
l = ll - 1
igh = ih
gxtmp(1:ngw) = csign(l) * wnl(1:ngw,iv,is) * gwork(1:ngw,iy)
DO ia = 1, atoms%na(is)
auxc(1:ngw,ia) = gxtmp(1:ngw) * eigr(1:ngw,iss + ia - 1)
END DO
IF( cdesc%gamma ) THEN
CALL DGEMM('T', 'N', atoms%na(is), nb, 2*ngw, 1.0d0, auxc(1,1), lda, &
c(1,1), ldw, 0.0d0, dfnl%r(iss,igh,1), ldf)
ELSE
CALL ZGEMM('C', 'N', atoms%na(is), nb, ngw, one, auxc(1,1), lda, &
c(1,1), ldw, zero, dfnl%c(iss,igh,1), ldf)
END IF
END DO LM
DEALLOCATE(auxc)
iss = iss + atoms%na(is)
END DO SPECS
IF( cdesc%gamma ) CALL DSCAL(size(dfnl%r),2.0d0,dfnl%r(1,1,1),1)
DEALLOCATE(gwork, gxtmp)
RETURN
END SUBROUTINE nlsm2_s
END MODULE nl
END MODULE nl

5
CPV/optical.f90
View File

@ -63,7 +63,7 @@
USE cell_module, ONLY: boxdimensions
USE wave_types, ONLY: wave_descriptor
USE pseudo_projector, ONLY: projector, allocate_projector, deallocate_projector
USE pseudopotential, ONLY: nsanl, ngh
USE pseudopotential, ONLY: nsanl
USE nl, ONLY: nlsm1
USE forces, ONLY: dforce_all
USE brillouin, ONLY: kpoints
@ -79,6 +79,7 @@
USE control_flags, ONLY: force_pairing
USE reciprocal_vectors, ONLY: gx, g
USE reciprocal_space_mesh, ONLY: gkx_l
USE uspp_param, ONLY: nhm
INTEGER, INTENT(IN) :: nfi
TYPE(boxdimensions), INTENT(IN) :: box
@ -149,7 +150,7 @@
END DO
ALLOCATE( eforce( ngw, nb_l, nk ) )
CALL allocate_projector(fnle, nsanl, nb_l, ngh, kp%gamma_only)
CALL allocate_projector(fnle, nsanl, nb_l, nhm, kp%gamma_only)
CALL nlsm1( ispin, ps%wnl(:,:,:,1), atoms, eigr, ce(:,:,1,ispin), wempt, g, &
gx, fnle(1))
CALL dforce_all( ispin, ce(:,:,:,ispin), wempt, ff, eforce, vpot(:,:,:,ispin), &

16
CPV/potentials.f90
View File

@ -274,6 +274,8 @@
USE reciprocal_vectors, ONLY: gx
USE gvecp, ONLY: ngm
USE pseudo_base, ONLY: compute_eself
USE local_pseudo, ONLY: vps, rhops
USE ions_base, ONLY: rcmax
IMPLICIT NONE
@ -335,7 +337,7 @@
REAL(dbl) :: omega, desr(6), pesum(16)
REAL(dbl) :: s0, s1, s2, s3, s4, s5, s6, s7, s8
REAL(dbl) :: rsum( SIZE( rhoe, 4 ) )
REAL(dbl) :: zv( atoms%nsp ), raggio( atoms%nsp )
REAL(dbl) :: zv( atoms%nsp )
LOGICAL :: ttstress, ttforce, ttscreen, ttsic, tgc
@ -628,7 +630,7 @@
ALLOCATE( vloc( ngm ) )
CALL vofloc(ttscreen, ttforce, edft%ehte, edft%ehti, ehp, &
eps, vloc, rhoeg, fion, atoms, ps%rhops, ps%vps, kp, eigr, &
eps, vloc, rhoeg, fion, atoms, rhops, vps, kp, eigr, &
ei1, ei2, ei3, sfac, box, desc, ps%ap )
! edft%ehte = REAL ( ehtep )
@ -723,9 +725,8 @@
! ... self interaction energy of the pseudocharges
do is = 1, atoms%nsp
zv( is ) = ps%ap(is)%zv
raggio( is ) = ps%ap(is)%raggio
end do
edft%eself = compute_eself( atoms%na, zv, raggio, atoms%nsp )
edft%eself = compute_eself( atoms%na, zv, rcmax, atoms%nsp )
IF( ttsic .and. ionode ) THEN
@ -752,11 +753,11 @@
IF(tprint.AND.tvhmean) THEN
IF(nspin.GT.2) CALL errore(' vofrho ',' spin + vmean ',nspin)
IF(nspin==1) CALL vofmean( sfac, ps%rhops, rhoeg(:,1) )
IF(nspin==1) CALL vofmean( sfac, rhops, rhoeg(:,1) )
IF(nspin==2) THEN
ALLOCATE(rho12(ngm))
rho12 (:) = rhoeg(:,1)+rhoeg(:,2)
CALL vofmean( sfac, ps%rhops, rho12)
CALL vofmean( sfac, rhops, rho12)
DEALLOCATE(rho12)
END IF
END IF
@ -1150,6 +1151,7 @@
USE descriptors_module, ONLY: global_index, local_dimension
USE atoms_type_module, ONLY: atoms_type
USE cp_types, ONLY: pseudo
USE ions_base, ONLY: rcmax
IMPLICIT NONE
@ -1224,7 +1226,7 @@
DO k = 1, atoms%nsp
DO j = k, atoms%nsp
zv2( k, j ) = ps%ap(k)%zv * ps%ap(j)%zv
rc ( k, j ) = SQRT( ps%ap(k)%raggio**2 + ps%ap(j)%raggio**2 )
rc ( k, j ) = SQRT( rcmax(k)**2 + rcmax(j)**2 )
END DO
END DO

20
CPV/problem_size.f90
View File

@ -24,14 +24,14 @@
use electrons_module, only: n_emp
use brillouin, only: get_kpoints_number
use reciprocal_vectors, only: ngwx, ngmx, ngmt
use pseudopotential, only: lmax, ngh
use uspp_param, only: nhm, lmaxkb
USE io_global, ONLY: ionode
USE io_global, ONLY: stdout
USE fft_base, ONLY: dfftp, dffts
implicit none
integer NGWXM_EMP,ngh_EMP,NSAX_EMP,NCHAINX, NPROC
integer NGWXM_EMP,nhm_EMP,NSAX_EMP,NCHAINX, NPROC
integer NNR1X, NNR2X,NNR3X,LM1X
integer NR1X, NR2X, NR3X, nr1_l,nr2_l,nr3_l
integer NSAX,NAMX,NAFX
@ -55,7 +55,7 @@
nk = get_kpoints_number()
ngh_EMP = ngh
nhm_EMP = nhm
NGWXM_EMP = ngwx
NNR1X = 2*NR1X-1
NNR2X = 2*NR2X-1
@ -63,7 +63,7 @@
NSAX = NSP*NAX
NAMX = NSAX
NAFX = NSAX
LM1X = LMAX -1
LM1X = lmaxkb
nbyte = 0
nbyte_alloc = 0
NSAX_EMP = NSAX
@ -80,11 +80,11 @@
nbyte = nbyte + 8 * ngmx * 3
! ... FNL
nbyte = nbyte + 8 * NSAX * NX * ngh * nspin
nbyte = nbyte + 8 * NSAX * NX * nhm * nspin
! ... WNL RHOPS VPS GNL RW RPS VR
nbyte = nbyte + 8*ndmx*NSP*3
nbyte = nbyte + 8*ngwx*ngh*NSP
nbyte = nbyte + 8*ngwx*nhm*NSP
nbyte = nbyte + 8*NSP*ngmx
nbyte = nbyte + 8*NSP*ngmx
nbyte = nbyte + 8*ndmx*NSP
@ -113,10 +113,10 @@
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... nlsm1 e nlsm2
itmp = 8 * 2*ngwx*NSAX + 8*NSAX*NX*ngh
itmp = 8 * 2*ngwx*NSAX + 8*NSAX*NX*nhm
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... eigsnew
itmp = 8 * ( 3*NX + NX*NX + N_EMP + NSAX_EMP*N_EMP*ngh_EMP &
itmp = 8 * ( 3*NX + NX*NX + N_EMP + NSAX_EMP*N_EMP*nhm_EMP &
& + NX*(NX+1)/2 + N_EMP*N_EMP + N_EMP + NX )
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... phfac
@ -124,8 +124,8 @@
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... pvofrho & pstress
itmp = 8 * (2*3*NAX*NSP+3*NAX*NSP+NR1_L*NR2_L*NR3_L*8 + &
& NSAX*NX*ngh*6 + 6*ngmx + ndmx + 6*ngwx + &
& ngwx*ngh*NSP + 2*ngwx*NSAX)
& NSAX*NX*nhm*6 + 6*ngmx + ndmx + 6*ngwx + &
& ngwx*nhm*NSP + 2*ngwx*NSAX)
if(itmp.gt.nbyte_alloc) nbyte_alloc = itmp
! ... formf
itmp = 8 * 3 * ndmx + 8 * 3 * ngmx

102
CPV/pseudo_base.f90
View File

@ -23,7 +23,7 @@
PRIVATE
PUBLIC :: nlin_base, nlin_stress_base, nlset_base
PUBLIC :: nlin_base, nlin_stress_base
PUBLIC :: compute_rhops, formfn, formfa
PUBLIC :: compute_eself, compute_rhocg
@ -68,25 +68,25 @@
! ... declare other variables
REAL(dbl), ALLOCATABLE :: fint(:,:)
REAL(dbl) :: xg, dx
INTEGER :: ig, mmax, lnl, l, ind
INTEGER :: ig, mmax, nbeta, l, ind
REAL(dbl), ALLOCATABLE :: ftest(:,:)
! ... end of declarations
! ----------------------------------------------
lnl = ap%lnl
nbeta = ap%nbeta
mmax = ap%mesh
dx = ap%dx
ALLOCATE(fint(mmax,lnl))
ALLOCATE(fint(mmax,nbeta))
! WRITE( stdout,*) 'DEBUG nlin_base, tpiba = ', tpiba
DO ig = 1, SIZE( wnl, 1 )
IF( hg(ig) < gsmall ) THEN
DO l = 1,lnl
DO l = 1,nbeta
! ... G=0 (Only if l=1, since otherwise the radial Bessel function jl=0)
IF( ap%indl(l) == 1 ) THEN
IF( ap%lll(l) == 0 ) THEN
fint(1:mmax,l) = ap%rw(1:mmax)**2 * ap%vrps(1:mmax,l)
call simpson_fpmd(mmax, fint(:,l), dx, wnl(ig,l))
! call simpson(mmax, fint(:,l), ap%rab, wnl(ig,l))
@ -97,8 +97,8 @@
ELSE
! ... Bessel functions: j_0(0)=1, j_n(0)=0 for n>0
xg = SQRT( hg(ig) ) * tpiba
CALL bessel2(xg, ap%rw, fint, lnl, ap%indl, mmax)
DO l = 1, lnl
CALL bessel2(xg, ap%rw, fint, nbeta, ap%lll, mmax)
DO l = 1, nbeta
fint(1:mmax,l) = fint(1:mmax,l) * ap%rw(1:mmax)**2 * ap%vrps(1:mmax,l)
call simpson_fpmd(mmax, fint(:,l), dx, wnl(ig,l))
! call simpson(mmax, fint(:,l), ap%rab, wnl(ig,l))
@ -133,24 +133,24 @@
! ... declare other variables
REAL(dbl), ALLOCATABLE :: fint(:,:)
REAL(dbl) xg, dx
INTEGER ig,mmax,gstart,lnl
INTEGER ig,mmax,gstart,nbeta
INTEGER l,ll
INTEGER ir
! ... end of declarations
! ----------------------------------------------
lnl = ap%lnl
nbeta = ap%nbeta
mmax = ap%mesh
dx = ap%dx
ALLOCATE( fint(mmax, lnl) )
ALLOCATE( fint(mmax, nbeta) )
DO ig = 1, SIZE( wnla, 1 )
IF( hg(ig) < gsmall ) THEN
! ... G=0 (Only if L=1, since otherwise the radial Bessel function JL=0)
DO l = 1, lnl
IF(ap%indl(l).EQ.1) THEN
DO l = 1, nbeta
IF( ap%lll(l) == 0 ) THEN
fint(1:mmax,l) = ap%rw(1:mmax)**2 * ap%vrps(1:mmax,l)
call simpson_fpmd(mmax, fint(:,l), dx, wnla(ig, l))
ELSE
@ -159,8 +159,8 @@
END DO
ELSE
xg = SQRT(hg(ig)) * tpiba
CALL bessel3(xg, ap%rw, fint, lnl, ap%indl, mmax)
DO l = 1, lnl
CALL bessel3(xg, ap%rw, fint, nbeta, ap%lll, mmax)
DO l = 1, nbeta
fint(1:mmax,l) = fint(1:mmax,l) * ap%rw(1:mmax)**2 * ap%vrps(1:mmax,l)
call simpson_fpmd(mmax, fint(:,l), dx, wnla(ig,l))
END DO
@ -172,78 +172,6 @@
RETURN
END SUBROUTINE nlin_stress_base
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE nlset_base(ap, wsgset)
! This subroutine computes the pseudopotential
! constants
! ap%rw**2 * ap%vrps = [ ( Vpsnl(r) - Vpsloc(r) )* Rps(r) * r^2 ]
! = [ DVpsnl(r) * Rps(r) * r^2 ]
! wsgset = 4 PI (2l+1) / < Rps(r) | DVpsnl(r) | Rps(r) >
!
! for all l and m
!
! This subroutine is executed only for non-local
! pseudopotentials.
! ----------------------------------------------
USE pseudo_types, ONLY: pseudo_ncpp
! ... declare subroutine arguments
TYPE (pseudo_ncpp), INTENT(IN) :: ap
REAL(dbl), INTENT(OUT) :: wsgset(:)
! ... declare other variables
INTEGER :: ir,i,igh2,l,ll,igh1,igh,mesh,igau,lnl,ltru
REAL(dbl) :: alt,blt,one_by_rcl
REAL(dbl) :: wsgl, dx
REAL(dbl), ALLOCATABLE :: gloc(:), fint(:)
! ... end of declarations
! ----------------------------------------------
lnl = ap%lnl
igau = ap%igau
dx = ap%dx
mesh = ap%mesh
wsgset = 0.0d0
ALLOCATE(fint(mesh))
! ... Set the normalizing factor "wsg"
igh2 = 0
DO l = 1, lnl
! find out the angular momentum (ll-1) of the component stored in position l
ll = ap%indl( l )
fint( 1:mesh ) = ap%rps( 1:mesh, ll ) * ap%vrps( 1:mesh, l ) * ap%rw( 1:mesh )
call simpson_fpmd( mesh, fint, dx, wsgl )
! ltru is the true angular momentum quantum number
ltru = ll - 1
igh1 = igh2 + 1
igh2 = igh1 + ( 2*ltru + 1 ) - 1
DO igh = igh1, igh2
! wsgset( igh ) = 4.0d0 * pi * ( 2*ltru + 1 ) / wsgl
wsgset( igh ) = 4.0d0 * pi / wsgl * ( 4.0d0 * pi )
END DO
END DO
DEALLOCATE(fint)
RETURN
END SUBROUTINE nlset_base
!-----------------------------------------------------------------------
subroutine compute_rhocg( rhocb, drhocb, r, rab, rho_atc, gb, omegab, &

601
CPV/pseudopot.f90
View File

@ -15,45 +15,22 @@
#include "f_defs.h"
! BEGIN manual
MODULE pseudopotential
! (describe briefly what this module does...)
! ----------------------------------------------
! routines in this module:
! SUBROUTINE readpot(na,nsp,tcc)
! SUBROUTINE defngh(nsp)
! INTEGER FUNCTION l2ind( lw, is )
! SUBROUTINE pseudopotential_setup(nsp,psfile)
! SUBROUTINE deallocate_pseudopotential
! SUBROUTINE formf(ht,ps)
! SUBROUTINE nlin(kp,wnl)
! SUBROUTINE nlin_stress(wnla)
! SUBROUTINE pseudo_wave_info(oc_out,nchan_out,mesh_out,dx_out, &
! rw_out,rps_out)
! REAL(dbl) FUNCTION zvpseudo(is)
! ----------------------------------------------
! END manual
! ... declare modules
USE kinds
USE parameters, ONLY: cp_lmax
USE cp_types, ONLY: pseudo, allocate_pseudo
USE splines, ONLY: spline_data
USE read_pseudo_module_fpmd, ONLY: nspnl, l2ind
USE read_pseudo_module_fpmd, ONLY: nspnl
IMPLICIT NONE
SAVE
! declare module-scope variables
INTEGER :: nsanl ! number of atoms of the non local species
INTEGER :: lnlx ! maximum number of different non local components
INTEGER :: lmax ! maximum value of the angular momentum
INTEGER :: lm1x ! maximum value of the angular momentum, of non local components
INTEGER :: ngh ! actual number of spherical harmonics
REAL(dbl), ALLOCATABLE :: wsginit(:,:)
INTEGER :: nbetax ! number of atoms of the non local species
TYPE (spline_data), ALLOCATABLE :: vps_sp(:)
TYPE (spline_data), ALLOCATABLE :: dvps_sp(:)
@ -67,53 +44,33 @@
LOGICAL :: tpstab, tpstab_first
LOGICAL :: ts, tp, td, tf, tl(0:(cp_lmax-1))
PRIVATE
PUBLIC :: pseudopotential_setup, formf, nlin, nlin_stress
PUBLIC :: pseudo_wave_info, pseudopotential_init
PUBLIC :: deallocate_pseudopotential
PUBLIC :: l2ind
PUBLIC :: nspnl, tl, lm1x, ngh, nsanl, ts, tp, td, tf, lnlx, lmax
PUBLIC :: nspnl, nsanl
PUBLIC :: pseudopotential_initval, pseudopotential_indexes
PUBLIC :: compute_dvan, compute_betagx, compute_qradx
! end of module-scope declarations
! ----------------------------------------------
CONTAINS
! subroutines
! ----------------------------------------------
! ----------------------------------------------
SUBROUTINE pseudopotential_setup(nsp, tpstab_inp, pstab_size_inp, raggio_inp)
SUBROUTINE pseudopotential_setup( nsp, tpstab_inp, pstab_size_inp )
! (describe briefly what this routine does...)
! ----------------------------------------------
USE splines, ONLY: nullify_spline
USE pseudo_base, ONLY: nlset_base
USE ions_base, ONLY: zv
USE pseudo_types, ONLY: pseudo_ncpp, pseudo_upf
USE pseudo_types, ONLY: pseudo_ncpp
USE read_pseudo_module_fpmd, ONLY: ap
USE splines, ONLY: kill_spline
INTEGER, INTENT(IN) :: nsp, pstab_size_inp
LOGICAL, INTENT(IN) :: tpstab_inp
REAL(dbl), INTENT(IN) :: raggio_inp(:)
INTEGER :: i, is, il, l, j
! end of declarations
! ----------------------------------------------
DO i = 1, nsp
ap(i)%raggio = raggio_inp(i)
IF( ap(i)%raggio <= 0.0d0 ) THEN
CALL errore(' pseudopotential_setup ',' ion_radius less than 0 ',-1)
END IF
END DO
INTEGER :: i
! initialize the global array zv containing the
! value of the ionic charges
@ -122,49 +79,36 @@
zv( i ) = ap( i )%zv
END DO
! ... set flags selecting angular momentum
lnlx = 0
ts=.FALSE.
tp=.FALSE.
td=.FALSE.
tf=.FALSE.
tl=.FALSE.
DO is = 1, nspnl
lnlx = MAX( lnlx, ap(is)%lnl )
DO l = 1, ap(is)%lnl
IF (ap(is)%indl(l).EQ.1) ts = .TRUE. ! state = s
IF (ap(is)%indl(l).EQ.2) tp = .TRUE. ! state = p
IF (ap(is)%indl(l).EQ.3) td = .TRUE. ! state = d
IF (ap(is)%indl(l).EQ.4) tf = .TRUE. ! state = f
END DO
END DO
tl(0) = ts
tl(1) = tp
tl(2) = td
tl(3) = tf
! set the sizes for the spline tables
! ... count orbitals
ngh = 0
IF (ts) ngh = ngh + 1
IF (tp) ngh = ngh + 3
IF (td) ngh = ngh + 5
IF (tf) ngh = ngh + 7
tpstab = tpstab_inp
pstab_size = pstab_size_inp
! maximum value of l (lmax = 1(s),2(p),3(d),4(f) )
! lmax - 1, number of non-local components (Projectors)
! one of the component is added to the local part
RETURN
END SUBROUTINE pseudopotential_setup
lm1x = 0
IF (ts) lm1x = 0
IF (tp) lm1x = 1
IF (td) lm1x = 2
IF (tf) lm1x = 3
lmax = lm1x
DO is = 1, nspnl
lmax = MAX( lmax, ( ap(is)%lloc - 1 ) )
END DO
tpstab = tpstab_inp
! ----------------------------------------------
SUBROUTINE pseudopotential_initval()
USE splines, ONLY: nullify_spline
USE ions_base, ONLY: zv, nsp
USE pseudo_types, ONLY: pseudo_ncpp, pseudo_upf
USE read_pseudo_module_fpmd, ONLY: ap
USE splines, ONLY: kill_spline
USE constants, ONLY: pi
use uspp_param, only: nbeta
INTEGER :: i, is, il, j
! ... set flags selecting angular momentum
nbetax = MAXVAL( nbeta( 1:nsp ) )
IF( tpstab ) THEN
!
IF( ALLOCATED( vps_sp ) ) THEN
@ -207,7 +151,7 @@
END DO
DEALLOCATE(wnl_sp)
END IF
ALLOCATE( wnl_sp(lnlx,nsp))
ALLOCATE( wnl_sp(nbetax,nsp))
!
IF(ALLOCATED(wnla_sp)) THEN
DO i = 1, size(wnla_sp,2)
@ -217,35 +161,171 @@
END DO
DEALLOCATE(wnla_sp)
END IF
ALLOCATE( wnla_sp(lnlx,nsp))
ALLOCATE( wnla_sp(nbetax,nsp))
!
DO is = 1, nsp
CALL nullify_spline( vps_sp( is ) )
CALL nullify_spline( dvps_sp( is ) )
CALL nullify_spline( rhoc1_sp( is ) )
CALL nullify_spline( rhocp_sp( is ) )
DO il = 1, lnlx
DO il = 1, nbetax
CALL nullify_spline( wnl_sp( il, is ) )
CALL nullify_spline( wnla_sp( il, is ) )
END DO
END DO
!
tpstab_first = .TRUE.
pstab_size = pstab_size_inp
!
END IF
IF( ALLOCATED( wsginit ) ) DEALLOCATE( wsginit )
ALLOCATE(wsginit(ngh,nsp))
wsginit = 0.0d0
DO is = 1, nspnl
CALL nlset_base(ap(is), wsginit(:,is))
END DO
call compute_dvan()
RETURN
END SUBROUTINE pseudopotential_initval
RETURN
END SUBROUTINE pseudopotential_setup
! ----------------------------------------------
!
! calculate array dvan(iv,jv,is)
!
! rw**2 * vrps = [ ( Vpsnl(r) - Vpsloc(r) )* Rps(r) * r^2 ]
! = [ DVpsnl(r) * Rps(r) * r^2 ]
! dion = (2l+1) / < Rps(r) | DVpsnl(r) | Rps(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
! ----------------------------------------------

6
CPV/pseudotab_base.f90
View File

@ -168,8 +168,4 @@
RETURN
END SUBROUTINE corecortab_base
! ----------------------------------------------
! ----------------------------------------------
END MODULE pseudotab_base
END MODULE pseudotab_base

261
CPV/read_pseudo.f90
View File

@ -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

10
CPV/rundiis.f90
View File

@ -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)

9
CPV/smcp.f90
View File

@ -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))

472
CPV/stress.f90
View File

@ -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

3
CPV/vofrho2.f90
View File

@ -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

1
CPV/wf.f90
View File

@ -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

5
Modules/pseudo_types.f90
View File

@ -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)