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

Step 1 of many: paw augmentation functions moved to upf structure 

(I apologise for the many check-ins I'll do, but I have to keep track of changes.)
LP


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@4424 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
paulatto 2007-11-12 16:20:55 +00:00
parent 3ee15febfe
commit 513379419b
12 changed files with 233 additions and 199 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
Modules/paw_variables.f90
View File

@ -15,8 +15,6 @@ MODULE paw_variables
!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!
!!!! Control flags: !!!! !!!! Control flags: !!!!
! debug: if .false. don't apply PAW correction
LOGICAL, PARAMETER :: really_do_paw = .true.
! Set to true after initialization, to prevent double allocs: ! Set to true after initialization, to prevent double allocs:
LOGICAL :: is_init = .false. LOGICAL :: is_init = .false.
! Analogous to okvan in "uspp_param" (Modules/uspp.f90) ! Analogous to okvan in "uspp_param" (Modules/uspp.f90)
@ -24,8 +22,8 @@ MODULE paw_variables
okpaw ! if .TRUE. at least one pseudo is PAW okpaw ! if .TRUE. at least one pseudo is PAW
! Analogous to tvanp in "uspp_param" (Modules/uspp.f90) ! Analogous to tvanp in "uspp_param" (Modules/uspp.f90)
LOGICAL :: & ! LOGICAL :: &
tpawp(npsx) = .false. ! if .TRUE. the atom is of PAW type ! tpawp(npsx) = .false. ! if .TRUE. the atom is of PAW type
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!! Initialization data: !!!! !!!! Initialization data: !!!!
@ -57,13 +55,13 @@ MODULE paw_variables
ptfunc(ndmx,nbrx,nbrx,npsx) ! PS: \tilde{\phi}_{mu}(|r|)-\tilde{\phi}_{nu}(|r|) ptfunc(ndmx,nbrx,nbrx,npsx) ! PS: \tilde{\phi}_{mu}(|r|)-\tilde{\phi}_{nu}(|r|)
! Augmentation on radial grid: ! Augmentation on radial grid:
TYPE augfun_t ! TYPE augfun_t
REAL(DP), ALLOCATABLE :: fun(:,:,:,:) ! REAL(DP), ALLOCATABLE :: fun(:,:,:,:)
END TYPE ! END TYPE
TYPE(augfun_t) :: aug(npsx) ! TYPE(augfun_t) :: aug(npsx)
! Moments of the augmentation functions ! Moments of the augmentation functions
REAL (DP) :: & REAL (DP) :: &
augmom(nbrx,nbrx,0:6,npsx) ! moments of PAW augm. functions augmom(nbrx,nbrx,0:6,npsx) ! moments of PAW augm. functions
INTEGER :: & INTEGER :: &
nraug(npsx) ! augm. functions cutoff parameter nraug(npsx) ! augm. functions cutoff parameter
@ -71,7 +69,7 @@ MODULE paw_variables
REAL(DP), TARGET :: & REAL(DP), TARGET :: &
aerho_atc(ndmx,npsx), &! radial AE core charge density aerho_atc(ndmx,npsx), &! radial AE core charge density
psrho_atc(ndmx,npsx) ! radial PS core charge density psrho_atc(ndmx,npsx) ! radial PS core charge density
! Analogous to vloc_at in "uspp_param" (Modules/uspp.f90) ! Analogous to vloc_at in "uspp_param" (Modules/uspp.f90)
! actually pseudopotential (AE and PS) on radial grid. ! actually pseudopotential (AE and PS) on radial grid.
REAL(DP), TARGET :: & REAL(DP), TARGET :: &
@ -100,13 +98,13 @@ MODULE paw_variables
saved(:) ! allocated in PAW_rad_init saved(:) ! allocated in PAW_rad_init
! This type contains some useful data that has to be passed to all ! This type contains some useful data that has to be passed to all
! the functions, but cannot stay in global variables for parallel: ! functions, but cannot stay in global variables for parallel:
TYPE paw_info TYPE paw_info
INTEGER :: a ! atom index INTEGER :: a ! atom index
INTEGER :: t ! atom type index INTEGER :: t ! atom type index
INTEGER :: m ! atom mesh = g(nt)%mesh INTEGER :: m ! atom mesh = g(nt)%mesh
INTEGER :: w ! w=1 --> all electron, w=2 --> pseudo INTEGER :: w ! number of atomic wavefunctions
! (used only for gradient correction) INTEGER :: l ! max angular index l
END TYPE END TYPE
! Analogous to deeq in "uspp_param" (Modules/uspp.f90) ! Analogous to deeq in "uspp_param" (Modules/uspp.f90)

29
Modules/pseudo_types.f90
View File

@ -60,9 +60,21 @@ TYPE :: paw_t
END TYPE paw_t END TYPE paw_t
! !
! Additional data to make a PAW setup out of an US pseudo,
! they are all stored on a radial grid:
TYPE paw_in_upf
REAL(DP),ALLOCATABLE :: aug(:,:,:,:) ! Augmentation charge
REAL(DP),ALLOCATABLE :: ae_core(:), &! AE core charge
ps_core(:) ! PS core charge
REAL(DP),ALLOCATABLE :: ae_vloc(:), &! AE local potential
ps_vloc(:) ! PS local potential
REAL(DP),ALLOCATABLE :: pfunc(:,:,:),&! Psi_i(r)*Psi_j(r)
ptfunc(:,:,:) ! as above, but for pseudo
END TYPE paw_in_upf
TYPE pseudo_upf TYPE pseudo_upf
CHARACTER(LEN=80):: generated ! CHARACTER(LEN=80):: generated !
CHARACTER(LEN=80):: date_author ! Misc info CHARACTER(LEN=80):: date_author ! Misc info
CHARACTER(LEN=80):: comment ! CHARACTER(LEN=80):: comment !
CHARACTER(LEN=2) :: psd ! Element label CHARACTER(LEN=2) :: psd ! Element label
@ -119,13 +131,18 @@ END TYPE paw_t
REAL(DP), POINTER :: qfcoef(:,:,:,:) ! qfcoef(nqf,0:2*lmax,nbeta,nbeta) REAL(DP), POINTER :: qfcoef(:,:,:,:) ! qfcoef(nqf,0:2*lmax,nbeta,nbeta)
! coefficients for Q for |r|<r_L ! coefficients for Q for |r|<r_L
REAL(DP), POINTER :: chi(:,:) ! chi(mesh,nwfc) atomic wavefcts REAL(DP), POINTER :: chi(:,:) ! chi(mesh,nwfc) atomic wavefcts
REAL(DP), POINTER :: rho_at(:) ! rho_at(mesh) atomic charge REAL(DP), POINTER :: rho_at(:) ! rho_at(mesh) atomic charge
LOGICAL :: has_paw ! Whether PAW data is included ! PAW:
LOGICAL :: has_paw ! Whether PAW data is included
REAL(DP) :: paw_data_format ! The version of the format REAL(DP) :: paw_data_format ! The version of the format
LOGICAL :: has_gipaw ! Whether GIPAW data is included LOGICAL :: tpawp ! true if atom is PAW
TYPE(paw_in_upf) :: paw ! additional data for PAW (see above)
! GIPAW:
LOGICAL :: has_gipaw ! Whether GIPAW data is included
REAL(DP) :: gipaw_data_format ! The version of the format REAL(DP) :: gipaw_data_format ! The version of the format
INTEGER :: gipaw_ncore_orbitals INTEGER :: gipaw_ncore_orbitals
REAL(DP), POINTER :: gipaw_core_orbital_n(:) REAL(DP), POINTER :: gipaw_core_orbital_n(:)
REAL(DP), POINTER :: gipaw_core_orbital_l(:) REAL(DP), POINTER :: gipaw_core_orbital_l(:)
CHARACTER(LEN=2), POINTER :: gipaw_core_orbital_el(:) CHARACTER(LEN=2), POINTER :: gipaw_core_orbital_el(:)

191
Modules/read_paw.f90
View File

@ -292,9 +292,9 @@ subroutine set_pseudo_paw (is, pawset)
USE pseudo_types USE pseudo_types
USE constants, ONLY: FPI USE constants, ONLY: FPI
! !
USE paw_variables, ONLY : tpawp, pfunc, ptfunc, aevloc_at, psvloc_at, & USE paw_variables, ONLY : pfunc, ptfunc, aevloc_at, psvloc_at, &
aerho_atc, psrho_atc, kdiff, & aerho_atc, psrho_atc, kdiff, &
augmom, nraug, step_f, aug augmom, nraug, step_f
! !
implicit none implicit none
! !
@ -306,6 +306,7 @@ subroutine set_pseudo_paw (is, pawset)
integer :: nb, mb, ijv integer :: nb, mb, ijv
TYPE (paw_t) :: pawset TYPE (paw_t) :: pawset
integer :: i,j, l, nrc, nrs integer :: i,j, l, nrc, nrs
integer :: nwfc, mesh
real (DP) :: pow real (DP) :: pow
real (DP), ALLOCATABLE :: aux (:) real (DP), ALLOCATABLE :: aux (:)
! !
@ -330,8 +331,7 @@ subroutine set_pseudo_paw (is, pawset)
upf(is)%psd = pawset%symbol upf(is)%psd = pawset%symbol
upf(is)%tvanp=.true. upf(is)%tvanp=.true.
upf(is)%nlcc = pawset%nlcc upf(is)%nlcc = pawset%nlcc
!tvanp(is)=.true. upf(is)%tpawp = .true.
tpawp(is)=.true.
call set_dft_from_name( pawset%dft ) call set_dft_from_name( pawset%dft )
! !
IF ( dft_is_meta() ) & IF ( dft_is_meta() ) &
@ -343,50 +343,52 @@ subroutine set_pseudo_paw (is, pawset)
#endif #endif
! !
rgrid(is)%mesh = pawset%grid%mesh rgrid(is)%mesh = pawset%grid%mesh
mesh = pawset%grid%mesh ! to make this file barely readable
nwfc = pawset%nwfc ! " " " " " "
! !
! ... Copy wavefunctions used for PAW construction. ! ... Copy wavefunctions used for PAW construction.
! ... Copy also the unoccupied ones, e.g. ! ... Copy also the unoccupied ones, e.g.
! ... corresponding to second energy for the same channel ! ... corresponding to second energy for the same channel
! ... (necessary to set starting occupations correctly) ! ... (necessary to set starting occupations correctly)
! !
upf(is)%nwfc = pawset%nwfc upf(is)%nwfc = nwfc
ALLOCATE ( upf(is)%lchi(pawset%nwfc), upf(is)%oc(pawset%nwfc) ) ALLOCATE ( upf(is)%lchi(nwfc), upf(is)%oc(nwfc) )
ALLOCATE ( upf(is)%chi( pawset%grid%mesh, pawset%nwfc) ) ALLOCATE ( upf(is)%chi( mesh, nwfc) )
do i=1, pawset%nwfc do i=1, nwfc
upf(is)%lchi(i)=pawset%l(i) upf(is)%lchi(i)=pawset%l(i)
upf(is)%oc(i)=MAX(pawset%oc(i),0._DP) upf(is)%oc(i)=MAX(pawset%oc(i),0._DP)
upf(is)%chi(1:pawset%grid%mesh, i) = pawset%pswfc(1:pawset%grid%mesh, i) upf(is)%chi(1:mesh, i) = pawset%pswfc(1:mesh, i)
end do end do
! !
upf(is)%nbeta= pawset%nwfc upf(is)%nbeta= nwfc
allocate ( upf(is)%kbeta(pawset%nwfc) ) allocate ( upf(is)%kbeta(nwfc) )
do nb=1,pawset%nwfc do nb=1,nwfc
upf(is)%kbeta(nb)=pawset%ikk(nb) upf(is)%kbeta(nb)=pawset%ikk(nb)
end do end do
! kkbeta is the maximum distance from nucleus where AE ! kkbeta is the maximum distance from nucleus where AE
! and PS wavefunction may not match: ! and PS wavefunction may not match:
upf(is)%kkbeta=MAXVAL (upf(is)%kbeta(:)) upf(is)%kkbeta=MAXVAL (upf(is)%kbeta(:))
allocate (upf(is)%beta(1:pawset%grid%mesh, 1:pawset%nwfc)) allocate (upf(is)%beta(1:mesh, 1:nwfc))
upf(is)%beta(1:pawset%grid%mesh, 1:pawset%nwfc) = & upf(is)%beta(1:mesh, 1:nwfc) = &
pawset%proj(1:pawset%grid%mesh, 1:pawset%nwfc) pawset%proj(1:mesh, 1:nwfc)
allocate(upf(is)%dion(1:pawset%nwfc, 1:pawset%nwfc)) allocate(upf(is)%dion(1:nwfc, 1:nwfc))
upf(is)%dion(1:pawset%nwfc, 1:pawset%nwfc) = pawset%dion(1:pawset%nwfc, 1:pawset%nwfc) upf(is)%dion(1:nwfc, 1:nwfc) = pawset%dion(1:nwfc, 1:nwfc)
kdiff(1:pawset%nwfc, 1:pawset%nwfc, is) = pawset%kdiff(1:pawset%nwfc, 1:pawset%nwfc) kdiff(1:nwfc, 1:nwfc, is) = pawset%kdiff(1:nwfc, 1:nwfc)
upf(is)%nqlc = 2*pawset%lmax+1 upf(is)%nqlc = 2*pawset%lmax+1
upf(is)%nqf = 0 !! no rinner, all numeric upf(is)%nqf = 0 !! no rinner, all numeric
allocate (upf(is)%lll(pawset%nwfc) ) allocate (upf(is)%lll(nwfc) )
upf(is)%lll(1:pawset%nwfc) = pawset%l(1:pawset%nwfc) upf(is)%lll(1:nwfc) = pawset%l(1:nwfc)
allocate (upf(is)%rinner(upf(is)%nqlc)) allocate (upf(is)%rinner(upf(is)%nqlc))
upf(is)%rinner(1:upf(is)%nqlc) = 0._dp !! no rinner, all numeric upf(is)%rinner(1:upf(is)%nqlc) = 0._dp !! no rinner, all numeric
! integral of augmentation charges vanishes for different values of l ! integral of augmentation charges vanishes for different values of l
allocate ( upf(is)%qqq(pawset%nwfc,pawset%nwfc)) allocate ( upf(is)%qqq(nwfc,nwfc))
do i = 1, pawset%nwfc do i = 1, nwfc
do j = 1, pawset%nwfc do j = 1, nwfc
if (pawset%l(i)==pawset%l(j)) then if (pawset%l(i)==pawset%l(j)) then
upf(is)%qqq(i,j) = pawset%augmom(i,j,0) !!gf spherical approximation upf(is)%qqq(i,j) = pawset%augmom(i,j,0) !!gf spherical approximation
else else
@ -398,45 +400,53 @@ subroutine set_pseudo_paw (is, pawset)
! import augmentation charge: ! import augmentation charge:
nraug(is) = pawset%irc nraug(is) = pawset%irc
rgrid(is)%dx = pawset%grid%dx rgrid(is)%dx = pawset%grid%dx
rgrid(is)%r2 (1:pawset%grid%mesh) = pawset%grid%r2 (1:pawset%grid%mesh) rgrid(is)%r2 (1:mesh) = pawset%grid%r2 (1:mesh)
do l = 0, 2*pawset%lmax do l = 0, 2*pawset%lmax
do i = 1, pawset%nwfc do i = 1, nwfc
do j = 1, pawset%nwfc do j = 1, nwfc
augmom(i,j,l,is) = pawset%augmom(i,j,l) augmom(i,j,l,is) = pawset%augmom(i,j,l)
end do end do
end do end do
end do end do
! FIXME!!! Probably this loop is unnecessary as the qfunc are reconstraucted from
! the full augfuncs !!
! triangularize matrix of qfunc's ! triangularize matrix of qfunc's
allocate ( upf(is)%qfunc(1:pawset%grid%mesh,pawset%nwfc*(pawset%nwfc+1)/2) ) allocate ( upf(is)%qfunc(1:mesh,nwfc*(nwfc+1)/2) )
do nb = 1, pawset%nwfc do nb = 1, nwfc
do mb = nb, pawset%nwfc do mb = nb, nwfc
ijv = mb * (mb-1) / 2 + nb ijv = mb * (mb-1) / 2 + nb
upf(is)%qfunc (1:pawset%grid%mesh,ijv) = & upf(is)%qfunc (1:mesh,ijv) = &
pawset%augfun(1:pawset%grid%mesh,nb,mb,0) pawset%augfun(1:mesh,nb,mb,0)
enddo enddo
enddo enddo
! new sparse allocation for augmentation functions ! new sparse allocation for augmentation functions
if (allocated(aug(is)%fun)) deallocate(aug(is)%fun) !if (allocated(aug(is)%fun)) deallocate(aug(is)%fun)
allocate(aug(is)%fun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax)) ! allocate(aug(is)%fun(1:mesh,1:nwfc,1:nwfc,0:2*pawset%lmax))
! ! !
aug(is)%fun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax) & ! aug(is)%fun(1:mesh,1:nwfc,1:nwfc,0:2*pawset%lmax) &
= & ! = &
pawset%augfun(1:pawset%grid%mesh,1:pawset%nwfc,1:pawset%nwfc,0:2*pawset%lmax) ! pawset%augfun(1:mesh,1:nwfc,1:nwfc,0:2*pawset%lmax)
! Augmentation functions for PAW depend on angular momentum!
allocate(upf(is)%paw%aug(1:mesh,1:nwfc,1:nwfc,0:2*pawset%lmax))
upf(is)%paw%aug(1:mesh,1:nwfc,1:nwfc,0:2*pawset%lmax) &
= pawset%augfun(1:mesh,1:nwfc,1:nwfc,0:2*pawset%lmax)
! !
#if ! defined __DO_NOT_CUTOFF_PAW_FUNC #if ! defined __DO_NOT_CUTOFF_PAW_FUNC
ALLOCATE (aux(1:pawset%grid%mesh)) ALLOCATE (aux(1:mesh))
#endif #endif
do i=1,pawset%nwfc do i=1,nwfc
do j=1,pawset%nwfc do j=1,nwfc
#if defined __DO_NOT_CUTOFF_PAW_FUNC #if defined __DO_NOT_CUTOFF_PAW_FUNC
pfunc (1:pawset%grid%mesh, i, j, is) = & pfunc (1:mesh, i, j, is) = &
pawset%aewfc(1:pawset%grid%mesh, i) * pawset%aewfc(1:pawset%grid%mesh, j) pawset%aewfc(1:mesh, i) * pawset%aewfc(1:mesh, j)
ptfunc (1:pawset%grid%mesh, i, j, is) = & ptfunc (1:mesh, i, j, is) = &
pawset%pswfc(1:pawset%grid%mesh, i) * pawset%pswfc(1:pawset%grid%mesh, j) pawset%pswfc(1:mesh, i) * pawset%pswfc(1:mesh, j)
#else #else
#if defined __SHARP_CUTOFF_PAW_FUNC #if defined __SHARP_CUTOFF_PAW_FUNC
pfunc(:,i,j,is) = 0._dp pfunc(:,i,j,is) = 0._dp
@ -446,14 +456,14 @@ subroutine set_pseudo_paw (is, pawset)
ptfunc (1:pawset%irc, i, j, is) = & ptfunc (1:pawset%irc, i, j, is) = &
pawset%pswfc(1:pawset%irc, i) * pawset%pswfc(1:pawset%irc, j) pawset%pswfc(1:pawset%irc, i) * pawset%pswfc(1:pawset%irc, j)
#else #else
aux(1:pawset%grid%mesh) = pawset%aewfc(1:pawset%grid%mesh, i) * & aux(1:mesh) = pawset%aewfc(1:mesh, i) * &
pawset%aewfc(1:pawset%grid%mesh, j) pawset%aewfc(1:mesh, j)
CALL step_f( pfunc(1:pawset%grid%mesh,i,j,is), aux(1:pawset%grid%mesh), & CALL step_f( pfunc(1:mesh,i,j,is), aux(1:mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh) pawset%grid%r(1:mesh), nrs, nrc, pow, mesh)
aux(1:pawset%grid%mesh) = pawset%pswfc(1:pawset%grid%mesh, i) * & aux(1:mesh) = pawset%pswfc(1:mesh, i) * &
pawset%pswfc(1:pawset%grid%mesh, j) pawset%pswfc(1:mesh, j)
CALL step_f( ptfunc(1:pawset%grid%mesh,i,j,is), aux(1:pawset%grid%mesh), & CALL step_f( ptfunc(1:mesh,i,j,is), aux(1:mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh) pawset%grid%r(1:mesh), nrs, nrc, pow, mesh)
#endif #endif
#endif #endif
end do end do
@ -463,80 +473,79 @@ subroutine set_pseudo_paw (is, pawset)
! ... One should not need \tilde{n}^1 alone in any case. ! ... One should not need \tilde{n}^1 alone in any case.
! !
! nqf is always 0 for this PAW format ! nqf is always 0 for this PAW format
! qfcoef(1:pawset%nqf, 1:pawset%nqlc, 1:pawset%nwfc, 1:pawset%nwfc, is ) = 0._dp ! qfcoef(1:pawset%nqf, 1:pawset%nqlc, 1:nwfc, 1:nwfc, is ) = 0._dp
! !
rgrid(is)%r (1:pawset%grid%mesh) = pawset%grid%r (1:pawset%grid%mesh) rgrid(is)%r (1:mesh) = pawset%grid%r (1:mesh)
rgrid(is)%rab(1:pawset%grid%mesh) = pawset%grid%r (1:pawset%grid%mesh)*pawset%grid%dx rgrid(is)%rab(1:mesh) = pawset%grid%r (1:mesh)*pawset%grid%dx
! !
! set radial grid data ! set radial grid data
! !
rgrid(is)%mesh = pawset%grid%mesh rgrid(is)%mesh = mesh
rgrid(is)%xmin = pawset%grid%xmin rgrid(is)%xmin = pawset%grid%xmin
rgrid(is)%rmax = pawset%grid%rmax rgrid(is)%rmax = pawset%grid%rmax
rgrid(is)%zmesh= pawset%grid%zmesh rgrid(is)%zmesh= pawset%grid%zmesh
rgrid(is)%dx = pawset%grid%dx rgrid(is)%dx = pawset%grid%dx
rgrid(is)%r(1:pawset%grid%mesh) = pawset%grid%r(1:pawset%grid%mesh) rgrid(is)%r(1:mesh) = pawset%grid%r(1:mesh)
rgrid(is)%r2(1:pawset%grid%mesh) = pawset%grid%r2(1:pawset%grid%mesh) rgrid(is)%r2(1:mesh) = pawset%grid%r2(1:mesh)
rgrid(is)%rab(1:pawset%grid%mesh) = pawset%grid%rab(1:pawset%grid%mesh) rgrid(is)%rab(1:mesh) = pawset%grid%rab(1:mesh)
rgrid(is)%sqr(1:pawset%grid%mesh) = sqrt(pawset%grid%r(1:pawset%grid%mesh)) rgrid(is)%sqr(1:mesh) = sqrt(pawset%grid%r(1:mesh))
! these speed up a lot a few calculations: ! these speed up a lot a few calculations (paw XC and GCXC):
rgrid(is)%rm1(1:pawset%grid%mesh) = pawset%grid%rm1(1:pawset%grid%mesh) rgrid(is)%rm1(1:mesh) = pawset%grid%rm1(1:mesh)
rgrid(is)%rm2(1:pawset%grid%mesh) = pawset%grid%rm2(1:pawset%grid%mesh) rgrid(is)%rm2(1:mesh) = pawset%grid%rm2(1:mesh)
rgrid(is)%rm3(1:pawset%grid%mesh) = pawset%grid%rm3(1:pawset%grid%mesh) rgrid(is)%rm3(1:mesh) = pawset%grid%rm3(1:mesh)
! NO spin orbit PAW implemented right now (oct 2007) ! NO spin orbit PAW implemented right now (oct 2007)
!!$ if (lspinorb.and..not.pawset%has_so) & !!$ if (lspinorb.and..not.pawset%has_so) &
!!$ call infomsg ('pawset_to_internal','At least one non s.o. pseudo') !!$ call infomsg ('pawset_to_internal','At least one non s.o. pseudo')
allocate (upf(is)%jchi(1:pawset%nwfc)) allocate (upf(is)%jchi(1:nwfc))
allocate (upf(is)%jjj(1:pawset%nwfc)) allocate (upf(is)%jjj(1:nwfc))
!!$ lspinorb=lspinorb.and.pawset%has_so !!$ lspinorb=lspinorb.and.pawset%has_so
!!$ if (pawset%has_so) then !!$ if (pawset%has_so) then
!!$ jchi(1:pawset%nwfc, is) = pawset%jchi(1:pawset%nwfc) !!$ jchi(1:nwfc, is) = pawset%jchi(1:nwfc)
!!$ jjj(1:pawset%nbeta, is) = pawset%jjj(1:pawset%nbeta) !!$ jjj(1:pawset%nbeta, is) = pawset%jjj(1:pawset%nbeta)
!!$ else !!$ else
upf(is)%jchi(1:pawset%nwfc) = 0._dp upf(is)%jchi(1:nwfc) = 0._dp
upf(is)%jjj(1:pawset%nwfc) = 0._dp upf(is)%jjj(1:nwfc) = 0._dp
!!$ endif !!$ endif
! !
if ( pawset%nlcc) then if ( pawset%nlcc) then
allocate ( upf(is)%rho_atc(pawset%grid%mesh) ) allocate ( upf(is)%rho_atc(mesh) )
upf(is)%rho_atc(1:pawset%grid%mesh) = & upf(is)%rho_atc(1:mesh) = &
pawset%psccharge(1:pawset%grid%mesh) & pawset%psccharge(1:mesh) / FPI / pawset%grid%r2(1:mesh)
/ FPI / pawset%grid%r2(1:pawset%grid%mesh)
end if end if
aerho_atc(1:pawset%grid%mesh, is) = pawset%aeccharge(1:pawset%grid%mesh) & aerho_atc(1:mesh, is) = pawset%aeccharge(1:mesh) &
& / FPI / pawset%grid%r2(1:pawset%grid%mesh) & / FPI / pawset%grid%r2(1:mesh)
if ( pawset%nlcc) then if ( pawset%nlcc) then
psrho_atc(1:pawset%grid%mesh, is) = pawset%psccharge(1:pawset%grid%mesh) & psrho_atc(1:mesh, is) = pawset%psccharge(1:mesh) &
& / FPI / pawset%grid%r2(1:pawset%grid%mesh) & / FPI / pawset%grid%r2(1:mesh)
else else
psrho_atc(:,is) = 0._dp psrho_atc(:,is) = 0._dp
end if end if
! !
allocate ( upf(is)%rho_at(pawset%grid%mesh) ) allocate ( upf(is)%rho_at(mesh) )
upf(is)%rho_at (1:pawset%grid%mesh) = pawset%pscharge(1:pawset%grid%mesh) upf(is)%rho_at (1:mesh) = pawset%pscharge(1:mesh)
allocate (upf(is)%vloc(1:pawset%grid%mesh)) allocate (upf(is)%vloc(1:mesh))
upf(is)%vloc(1:pawset%grid%mesh) = pawset%psloc(1:pawset%grid%mesh) upf(is)%vloc(1:mesh) = pawset%psloc(1:mesh)
#if defined __DO_NOT_CUTOFF_PAW_FUNC #if defined __DO_NOT_CUTOFF_PAW_FUNC
aevloc_at(1:pawset%grid%mesh,is) = pawset%aeloc(1:pawset%grid%mesh) aevloc_at(1:mesh,is) = pawset%aeloc(1:mesh)
psvloc_at(1:pawset%grid%mesh,is) = pawset%psloc(1:pawset%grid%mesh) psvloc_at(1:mesh,is) = pawset%psloc(1:mesh)
#else #else
#if defined __SHARP_CUTOFF_PAW_FUNC #if defined __SHARP_CUTOFF_PAW_FUNC
aevloc_at(1:pawset%grid%mesh,is) = pawset%aeloc(1:pawset%grid%mesh) aevloc_at(1:mesh,is) = pawset%aeloc(1:mesh)
psvloc_at(1:pawset%grid%mesh,is) = pawset%psloc(1:pawset%grid%mesh) psvloc_at(1:mesh,is) = pawset%psloc(1:mesh)
aux(1:pawset%grid%mesh) = pawset%aeloc(1:pawset%grid%mesh) aux(1:mesh) = pawset%aeloc(1:mesh)
#else #else
CALL step_f( aevloc_at(1:pawset%grid%mesh,is), aux(1:pawset%grid%mesh), & CALL step_f( aevloc_at(1:mesh,is), aux(1:mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh) pawset%grid%r(1:mesh), nrs, nrc, pow, mesh)
aux(1:pawset%grid%mesh) = pawset%psloc(1:pawset%grid%mesh) aux(1:mesh) = pawset%psloc(1:mesh)
CALL step_f( psvloc_at(1:pawset%grid%mesh,is), aux(1:pawset%grid%mesh), & CALL step_f( psvloc_at(1:mesh,is), aux(1:mesh), &
pawset%grid%r(1:pawset%grid%mesh), nrs, nrc, pow, pawset%grid%mesh) pawset%grid%r(1:mesh), nrs, nrc, pow, mesh)
DEALLOCATE (aux) DEALLOCATE (aux)
#endif #endif
#endif #endif

23
PW/electrons.f90
View File

@ -63,7 +63,7 @@ SUBROUTINE electrons()
USE mp_global, ONLY : intra_pool_comm, npool USE mp_global, ONLY : intra_pool_comm, npool
USE mp, ONLY : mp_sum USE mp, ONLY : mp_sum
! !
USE paw_variables, ONLY : really_do_paw, okpaw, tpawp, becnew USE paw_variables, ONLY : okpaw, becnew
USE paw_onecenter, ONLY : PAW_potential,PAW_integrate USE paw_onecenter, ONLY : PAW_potential,PAW_integrate
USE uspp, ONLY : becsum ! used for PAW USE uspp, ONLY : becsum ! used for PAW
USE uspp_param, ONLY : nhm ! used for PAW USE uspp_param, ONLY : nhm ! used for PAW
@ -175,11 +175,8 @@ SUBROUTINE electrons()
! !
IF (okpaw) THEN IF (okpaw) THEN
IF ( .not. ALLOCATED(becstep) ) ALLOCATE (becstep(nhm*(nhm+1)/2,nat,nspin)) IF ( .not. ALLOCATED(becstep) ) ALLOCATE (becstep(nhm*(nhm+1)/2,nat,nspin))
becstep (:,:,:) = 0._dp !becstep (:,:,:) = 0._dp
! DO na = 1, nat becstep(:,:,:) = becsum(:,:,:)
! IF (tpawp(ityp(na))) becstep(:,na,:) = becsum(:,na,:)
! END DO
becstep(:,:,:) = becsum(:,:,:)
END IF END IF
call create_scf_type ( rhoin ) call create_scf_type ( rhoin )
! !
@ -459,14 +456,12 @@ SUBROUTINE electrons()
etot = eband + ( etxc - etxcc ) + ewld + ehart + deband + demet + descf etot = eband + ( etxc - etxcc ) + ewld + ehart + deband + demet + descf
! !
IF (okpaw) THEN IF (okpaw) THEN
correction1c = (deband_PAW + descf_PAW + e_PAW) correction1c = (deband_PAW + descf_PAW + e_PAW)
PRINT '(5x,A,f12.6,A)', 'PAW correction: ',correction1c, ', composed of: ' PRINT '(5x,A,f12.6,A)', 'PAW correction: ',correction1c, ', composed of: '
PRINT '(5x,A,f10.6,A,f10.6,A,f12.6)',& PRINT '(5x,A,f10.6,A,f10.6,A,f12.6)',&
'de_band: ',deband_PAW,', de_scf: ',descf_PAW,', 1-center E: ', e_PAW 'de_band: ',deband_PAW,', de_scf: ',descf_PAW,', 1-center E: ', e_PAW
IF (really_do_paw) THEN etot = etot + correction1c
etot = etot + correction1c hwf_energy = hwf_energy + correction1c
hwf_energy = hwf_energy + correction1c
ENDIF
END IF END IF
! !
#if defined (EXX) #if defined (EXX)

41
PW/init_us_1.f90
View File

@ -27,27 +27,27 @@ subroutine init_us_1
! g) It computes the q terms which define the S matrix. ! g) It computes the q terms which define the S matrix.
! h) It fills the interpolation table for the beta functions ! h) It fills the interpolation table for the beta functions
! !
USE kinds, ONLY : DP USE kinds, ONLY : DP
USE parameters, ONLY : lmaxx USE parameters, ONLY : lmaxx
USE constants, ONLY : fpi USE constants, ONLY : fpi
USE atom, ONLY : rgrid USE atom, ONLY : rgrid
USE ions_base, ONLY : ntyp => nsp USE ions_base, ONLY : ntyp => nsp
USE cell_base, ONLY : omega, tpiba USE cell_base, ONLY : omega, tpiba
USE gvect, ONLY : g, gg USE gvect, ONLY : g, gg
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE us, ONLY : nqxq, dq, nqx, tab, tab_d2y, qrad, spline_ps USE us, ONLY : nqxq, dq, nqx, tab, tab_d2y, qrad, spline_ps
USE splinelib USE splinelib
USE uspp, ONLY : nhtol, nhtoj, nhtolm, dvan, qq, indv, ap, aainit, & USE uspp, ONLY : nhtol, nhtoj, nhtolm, dvan, qq, indv,&
qq_so, dvan_so, okvan ap, aainit, qq_so, dvan_so, okvan
USE uspp_param, ONLY : upf, lmaxq, nbetam, nh, nhm, lmaxkb USE uspp_param, ONLY : upf, lmaxq, nbetam, nh, nhm, lmaxkb
USE spin_orb, ONLY : lspinorb, so, rot_ylm, fcoef USE spin_orb, ONLY : lspinorb, so, rot_ylm, fcoef
USE paw_variables, ONLY : really_do_paw, okpaw, tpawp, aug USE paw_variables,ONLY : okpaw
USE paw_init, ONLY : PAW_post_init
! !
implicit none implicit none
! !
! here a few local variables ! here a few local variables
! !
integer :: nt, ih, jh, nb, mb, ijv, l, m, ir, iq, is, startq, & integer :: nt, ih, jh, nb, mb, ijv, l, m, ir, iq, is, startq, &
lastq, ilast, ndm lastq, ilast, ndm
! various counters ! various counters
@ -242,8 +242,8 @@ subroutine init_us_1
do mb = nb, upf(nt)%nbeta do mb = nb, upf(nt)%nbeta
ijv = mb * (mb-1) / 2 + nb ijv = mb * (mb-1) / 2 + nb
paw : & ! in PAW formalism aug. charge is computed elsewhere paw : & ! in PAW formalism aug. charge is computed elsewhere
if (tpawp(nt)) then if (upf(nt)%tpawp) then
qtot(1:upf(nt)%kkbeta,ijv) = aug(nt)%fun(1:upf(nt)%kkbeta,nb,mb,l) qtot(1:upf(nt)%kkbeta,ijv) = upf(nt)%paw%aug(1:upf(nt)%kkbeta,nb,mb,l)
else else
if ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .and. & if ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .and. &
( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .and. & ( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .and. &
@ -281,7 +281,7 @@ subroutine init_us_1
if ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .and. & if ( ( l >= abs(upf(nt)%lll(nb) - upf(nt)%lll(mb)) ) .and. &
( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .and. & ( l <= upf(nt)%lll(nb) + upf(nt)%lll(mb) ) .and. &
(mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) .or.& (mod (l+upf(nt)%lll(nb)+upf(nt)%lll(mb), 2) == 0) .or.&
tpawp(nt) ) then upf(nt)%tpawp ) then
do ir = 1, upf(nt)%kkbeta do ir = 1, upf(nt)%kkbeta
aux1 (ir) = aux (ir) * qtot (ir, ijv) aux1 (ir) = aux (ir) * qtot (ir, ijv)
enddo enddo
@ -405,6 +405,11 @@ subroutine init_us_1
deallocate (besr) deallocate (besr)
deallocate (aux1) deallocate (aux1)
deallocate (aux) deallocate (aux)
! Some paw variables are used on all nodes only for initialization,
! then they can be deallocated on all nodes but one.
call PAW_post_init()
call stop_clock ('init_us_1') call stop_clock ('init_us_1')
return return
end subroutine init_us_1 end subroutine init_us_1

3
PW/newd.f90
View File

@ -39,8 +39,7 @@ SUBROUTINE newd_g()
USE spin_orb, ONLY : lspinorb, so, domag USE spin_orb, ONLY : lspinorb, so, domag
USE noncollin_module, ONLY : noncolin USE noncollin_module, ONLY : noncolin
! !
USE paw_variables, ONLY : really_do_paw, okpaw, tpawp, & USE paw_variables, ONLY : okpaw, kdiff, dpaw_ae, dpaw_ps
kdiff, dpaw_ae, dpaw_ps
USE paw_onecenter, ONLY : PAW_newd USE paw_onecenter, ONLY : PAW_newd
USE uspp, ONLY : nhtol, nhtolm USE uspp, ONLY : nhtol, nhtolm
! !

28
PW/paw_init.f90
View File

@ -16,6 +16,8 @@ MODULE paw_init
PUBLIC :: PAW_init_becsum PUBLIC :: PAW_init_becsum
PUBLIC :: PAW_init_onecenter PUBLIC :: PAW_init_onecenter
PUBLIC :: PAW_post_init
PUBLIC :: allocate_paw_internals, deallocate_paw_internals PUBLIC :: allocate_paw_internals, deallocate_paw_internals
!!!========================================================================= !!!=========================================================================
@ -42,7 +44,8 @@ MODULE paw_init
! Called from clean_pw ! Called from clean_pw
SUBROUTINE deallocate_paw_internals SUBROUTINE deallocate_paw_internals
USE paw_variables USE paw_variables
USE ions_base, ONLY : ntyp => nsp USE uspp_param, ONLY : upf
USE ions_base, ONLY : ntyp => nsp
! !
IMPLICIT NONE IMPLICIT NONE
INTEGER :: nt INTEGER :: nt
@ -54,22 +57,29 @@ MODULE paw_init
! !
! Allocated in read_paw: ! Allocated in read_paw:
DO nt = 1,ntyp DO nt = 1,ntyp
IF(allocated(aug(nt)%fun)) DEALLOCATE (aug(nt)%fun) IF(allocated(upf(nt)%paw%aug)) DEALLOCATE (upf(nt)%paw%aug)
ENDDO ENDDO
! !
END SUBROUTINE deallocate_paw_internals END SUBROUTINE deallocate_paw_internals
! Deallocate variables that are used only at init and then no more necessary.
SUBROUTINE PAW_post_init()
! this routine does nothing at this moment...
RETURN
END SUBROUTINE PAW_post_init
! Initialize becsum with atomic occupations (for PAW atoms only) ! Initialize becsum with atomic occupations (for PAW atoms only)
! Notice: requires exact correspondence chi <--> beta in the atom, ! Notice: requires exact correspondence chi <--> beta in the atom,
! that is that all wavefunctions considered for PAW generation are ! that is that all wavefunctions considered for PAW generation are
! counted in chi (otherwise the array "oc" does not correspond to beta) ! counted in chi (otherwise the array "oc" does not correspond to beta)
SUBROUTINE PAW_init_becsum() SUBROUTINE PAW_init_becsum()
USE kinds, ONLY : DP
USE uspp, ONLY : becsum, nhtol, indv USE uspp, ONLY : becsum, nhtol, indv
USE uspp_param, ONLY : upf, nh USE uspp_param, ONLY : upf, nh, upf
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE lsda_mod, ONLY : nspin, starting_magnetization USE lsda_mod, ONLY : nspin, starting_magnetization
USE paw_variables, ONLY : tpawp, okpaw USE paw_variables, ONLY : okpaw
USE paw_onecenter, ONLY : PAW_symmetrize USE paw_onecenter, ONLY : PAW_symmetrize
IMPLICIT NONE IMPLICIT NONE
INTEGER :: ispin, na, nt, ijh, ih, jh, nb, mb INTEGER :: ispin, na, nt, ijh, ih, jh, nb, mb
@ -80,7 +90,7 @@ SUBROUTINE PAW_init_becsum()
! !
na_loop: DO na = 1, nat na_loop: DO na = 1, nat
nt = ityp(na) nt = ityp(na)
is_paw: IF (tpawp(nt)) THEN is_paw: IF (upf(nt)%tpawp) THEN
! !
ijh = 1 ijh = 1
ih_loop: DO ih = 1, nh(nt) ih_loop: DO ih = 1, nh(nt)
@ -121,9 +131,9 @@ END SUBROUTINE PAW_init_becsum
! calls PAW_rad_init to initialize onecenter integration. ! calls PAW_rad_init to initialize onecenter integration.
SUBROUTINE PAW_init_onecenter() SUBROUTINE PAW_init_onecenter()
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE paw_variables, ONLY : tpawp, xlm, saved USE paw_variables, ONLY : xlm, saved
USE atom, ONLY : g => rgrid USE atom, ONLY : g => rgrid
USE uspp_param, ONLY : lmaxq USE uspp_param, ONLY : lmaxq, upf
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE funct, ONLY : dft_is_gradient USE funct, ONLY : dft_is_gradient
@ -138,7 +148,7 @@ SUBROUTINE PAW_init_onecenter()
DO na = first_nat, last_nat DO na = first_nat, last_nat
nt = ityp(na) nt = ityp(na)
! note that if the atom is not paw it is left unallocated ! note that if the atom is not paw it is left unallocated
IF ( tpawp(nt) ) THEN IF ( upf(nt)%tpawp ) THEN
IF (allocated(saved(na)%v)) DEALLOCATE(saved(na)%v) IF (allocated(saved(na)%v)) DEALLOCATE(saved(na)%v)
ALLOCATE( saved(na)%v(g(nt)%mesh, lmaxq**2, nspin, 2 ) ) ALLOCATE( saved(na)%v(g(nt)%mesh, lmaxq**2, nspin, 2 ) )
! {AE|PS} ! {AE|PS}

80
PW/paw_onecenter.f90
View File

@ -52,6 +52,7 @@ MODULE paw_onecenter
USE parameters, ONLY : ntypx USE parameters, ONLY : ntypx
USE paw_variables, ONLY : paw_info, saved, ww, nx, lm_max, ylm, & USE paw_variables, ONLY : paw_info, saved, ww, nx, lm_max, ylm, &
xlm, dylmp, dylmt, sin_th, cos_th xlm, dylmp, dylmt, sin_th, cos_th
USE uspp_param, ONLY : upf
! !
IMPLICIT NONE IMPLICIT NONE
@ -80,9 +81,8 @@ SUBROUTINE PAW_potential(becsum, energy, e_cmp)
USE atom, ONLY : g => rgrid USE atom, ONLY : g => rgrid
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE uspp_param, ONLY : nhm, lmaxq USE uspp_param, ONLY : nhm, lmaxq, upf
USE paw_variables, ONLY : pfunc, ptfunc, tpawp, & USE paw_variables, ONLY : pfunc, ptfunc, aerho_atc, psrho_atc
aerho_atc, psrho_atc, aug
REAL(DP), INTENT(IN) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupations REAL(DP), INTENT(IN) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupations
REAL(DP),INTENT(OUT),OPTIONAL :: energy ! if present compute E[rho] REAL(DP),INTENT(OUT),OPTIONAL :: energy ! if present compute E[rho]
@ -112,8 +112,10 @@ SUBROUTINE PAW_potential(becsum, energy, e_cmp)
i%a = na ! the index of the atom i%a = na ! the index of the atom
i%t = ityp(na) ! the type of atom na i%t = ityp(na) ! the type of atom na
i%m = g(i%t)%mesh! radial mesh size for atom na i%m = g(i%t)%mesh! radial mesh size for atom na
i%w = upf(i%t)%nwfc
i%l = upf(i%t)%nqlc-1
! !
ifpaw: IF (tpawp(i%t)) THEN ifpaw: IF (upf(i%t)%tpawp) THEN
! !
! Arrays are allocated inside the cycle to allow reduced ! Arrays are allocated inside the cycle to allow reduced
! memory usage as differnt atoms have different meshes (they ! memory usage as differnt atoms have different meshes (they
@ -122,7 +124,6 @@ SUBROUTINE PAW_potential(becsum, energy, e_cmp)
ALLOCATE(v_lm(i%m,lmaxq**2,nspin)) ALLOCATE(v_lm(i%m,lmaxq**2,nspin))
! !
whattodo: DO i_what = AE, PS whattodo: DO i_what = AE, PS
i%w = i_what ! spherical_gradient likes to know
! STEP: 1 [ build rho_lm (PAW_rho_lm) ] ! STEP: 1 [ build rho_lm (PAW_rho_lm) ]
NULLIFY(rho_core) NULLIFY(rho_core)
IF (i_what == AE) THEN IF (i_what == AE) THEN
@ -134,31 +135,31 @@ SUBROUTINE PAW_potential(becsum, energy, e_cmp)
! sign to sum up the enrgy ! sign to sum up the enrgy
sgn = +1._dp sgn = +1._dp
ELSE ELSE
CALL PAW_rho_lm(i, becsum, ptfunc, rho_lm, aug) CALL PAW_rho_lm(i, becsum, ptfunc, rho_lm, upf(i%t)%paw%aug)
! optional argument for pseudo part --> ^^^ ! optional argument for pseudo part --> ^^^
rho_core => psrho_atc ! as before rho_core => psrho_atc ! as before
sgn = -1._dp ! as before sgn = -1._dp ! as before
ENDIF ENDIF
! cleanup previously stored potentials ! cleanup previously stored potentials
saved(i%a)%v(:,:,:,i%w) = 0._dp saved(i%a)%v(:,:,:,i_what) = 0._dp
! First compute the Hartree potential (it does not depend on spin...): ! First compute the Hartree potential (it does not depend on spin...):
CALL PAW_h_potential(i, rho_lm, v_lm(:,:,1), energy) CALL PAW_h_potential(i, rho_lm, v_lm(:,:,1), energy)
! using "energy" as the in/out parameter I save a double call, but I have to do this: ! using "energy" as the in/out parameter I save a double call, but I have to do this:
IF (present(energy)) energy_h = energy IF (present(energy)) energy_h = energy
DO is = 1,nspin ! ... so it has to be copied to all spin components DO is = 1,nspin ! ... so it has to be copied to all spin components
saved(i%a)%v(:,:,is,i%w) = v_lm(:,:,1) saved(i%a)%v(:,:,is,i_what) = v_lm(:,:,1)
ENDDO ENDDO
! Than the XC one: ! Than the XC one:
CALL PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy) CALL PAW_xc_potential(i, rho_lm, rho_core, v_lm, energy)
IF (present(energy)) energy_xc = energy IF (present(energy)) energy_xc = energy
saved(i%a)%v(:,:,:,i%w) = saved(i%a)%v(:,:,:,i_what) & saved(i%a)%v(:,:,:,i_what) = saved(i%a)%v(:,:,:,i_what) &
+ v_lm(:,:,:) + v_lm(:,:,:)
IF (present(energy)) energy_tot = energy_tot + sgn*(energy_xc + energy_h) IF (present(energy)) energy_tot = energy_tot + sgn*(energy_xc + energy_h)
IF (present(e_cmp)) THEN IF (present(e_cmp)) THEN
e_cmp(na, 1, i%w) = energy_xc e_cmp(na, 1, i_what) = energy_xc
e_cmp(na, 2, i%w) = energy_h e_cmp(na, 2, i_what) = energy_h
ENDIF ENDIF
ENDDO whattodo ENDDO whattodo
DEALLOCATE(rho_lm, v_lm) DEALLOCATE(rho_lm, v_lm)
@ -186,8 +187,7 @@ SUBROUTINE PAW_symmetrize(becsum)
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE symme, ONLY : nsym, irt, d1, d2, d3 USE symme, ONLY : nsym, irt, d1, d2, d3
USE uspp, ONLY : nhtolm,nhtol USE uspp, ONLY : nhtolm,nhtol
USE uspp_param, ONLY : nh USE uspp_param, ONLY : nh, upf
USE paw_variables, ONLY : tpawp
USE wvfct, ONLY : gamma_only USE wvfct, ONLY : gamma_only
USE control_flags, ONLY : nosym USE control_flags, ONLY : nosym
@ -260,7 +260,7 @@ SUBROUTINE PAW_symmetrize(becsum)
DO na = 1, nat DO na = 1, nat
nt = ityp(na) nt = ityp(na)
! No need to symmetrize non-PAW atoms ! No need to symmetrize non-PAW atoms
IF ( .not. tpawp(nt) ) CYCLE IF ( .not. upf(nt)%tpawp ) CYCLE
! !
DO ih = 1, nh(nt) DO ih = 1, nh(nt)
DO jh = ih, nh(nt) ! note: jh >= ih DO jh = ih, nh(nt) ! note: jh >= ih
@ -352,9 +352,9 @@ SUBROUTINE PAW_newd(d_ae, d_ps)
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE atom, ONLY : g => rgrid USE atom, ONLY : g => rgrid
USE paw_variables, ONLY : pfunc, ptfunc, tpawp, aug, & USE paw_variables, ONLY : pfunc, ptfunc, &
aevloc_at, psvloc_at, ra=>nraug aevloc_at, psvloc_at, ra=>nraug
USE uspp_param, ONLY : nh, nhm, lmaxq USE uspp_param, ONLY : nh, nhm, lmaxq, upf
REAL(DP), TARGET, INTENT(INOUT) :: d_ae( nhm, nhm, nat, nspin) REAL(DP), TARGET, INTENT(INOUT) :: d_ae( nhm, nhm, nat, nspin)
REAL(DP), TARGET, INTENT(INOUT) :: d_ps( nhm, nhm, nat, nspin) REAL(DP), TARGET, INTENT(INOUT) :: d_ps( nhm, nhm, nat, nspin)
@ -393,14 +393,15 @@ SUBROUTINE PAW_newd(d_ae, d_ps)
i%a = na ! the index of the atom i%a = na ! the index of the atom
i%t = ityp(na) ! the type of atom na i%t = ityp(na) ! the type of atom na
! skip non-paw atoms ! skip non-paw atoms
IF ( .not. tpawp(i%t) ) CYCLE IF ( .not. upf(i%t)%tpawp ) CYCLE
i%m = g(i%t)%mesh! radial mesh size for atom na i%m = g(i%t)%mesh! radial mesh size for atom na
i%w = upf(i%t)%nwfc
i%l = upf(i%t)%nqlc-1
! !
! Different atoms may have different mesh sizes: ! Different atoms may have different mesh sizes:
ALLOCATE(pfunc_lm(i%m, lmaxq**2, nspin)) ALLOCATE(pfunc_lm(i%m, lmaxq**2, nspin))
! !
whattodo: DO i_what = AE, PS whattodo: DO i_what = AE, PS
i%w = i_what
! !
spins: DO is = 1, nspin spins: DO is = 1, nspin
nmb = 0 nmb = 0
@ -411,12 +412,12 @@ SUBROUTINE PAW_newd(d_ae, d_ps)
! !
! compute the density from a single pfunc ! compute the density from a single pfunc
becfake(nmb,na,is) = 1._dp becfake(nmb,na,is) = 1._dp
IF (i%w == AE) THEN IF (i_what == AE) THEN
CALL PAW_rho_lm(i, becfake, pfunc, pfunc_lm) CALL PAW_rho_lm(i, becfake, pfunc, pfunc_lm)
v_at => aevloc_at v_at => aevloc_at
d => d_ae d => d_ae
ELSE ELSE
CALL PAW_rho_lm(i, becfake, ptfunc, pfunc_lm, aug) CALL PAW_rho_lm(i, becfake, ptfunc, pfunc_lm, upf(i%t)%paw%aug)
v_at => psvloc_at v_at => psvloc_at
d => d_ps d => d_ps
ENDIF ENDIF
@ -426,10 +427,10 @@ SUBROUTINE PAW_newd(d_ae, d_ps)
DO lm = 1,lmaxq**2 DO lm = 1,lmaxq**2
IF ( lm == 1 ) THEN IF ( lm == 1 ) THEN
pfunc_lm(1:i%m,lm,is) = pfunc_lm(1:i%m,lm,is) * & pfunc_lm(1:i%m,lm,is) = pfunc_lm(1:i%m,lm,is) * &
( saved(i%a)%v(1:i%m,lm,is,i%w) + e2*sqrtpi*v_at(1:i%m,i%t) ) ( saved(i%a)%v(1:i%m,lm,is,i_what) + e2*sqrtpi*v_at(1:i%m,i%t) )
ELSE ELSE
pfunc_lm(1:i%m,lm,is) = pfunc_lm(1:i%m,lm,is) * & pfunc_lm(1:i%m,lm,is) = pfunc_lm(1:i%m,lm,is) * &
saved(i%a)%v(1:i%m,lm,is,i%w) saved(i%a)%v(1:i%m,lm,is,i_what)
ENDIF ENDIF
! !
! Integrate! ! Integrate!
@ -470,8 +471,8 @@ FUNCTION PAW_integrate(becsum)
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE atom, ONLY : g => rgrid USE atom, ONLY : g => rgrid
USE paw_variables, ONLY : pfunc, ptfunc, tpawp, aug USE paw_variables, ONLY : pfunc, ptfunc
USE uspp_param, ONLY : nhm, lmaxq USE uspp_param, ONLY : nhm, lmaxq, upf
REAL(DP) :: PAW_integrate REAL(DP) :: PAW_integrate
@ -499,19 +500,20 @@ FUNCTION PAW_integrate(becsum)
i%a = na ! the index of the atom i%a = na ! the index of the atom
i%t = ityp(na) ! the type of atom na i%t = ityp(na) ! the type of atom na
i%m = g(i%t)%mesh! radial mesh size for atom na i%m = g(i%t)%mesh! radial mesh size for atom na
i%w = upf(i%t)%nwfc
i%l = upf(i%t)%nqlc-1
! !
ifpaw: IF (tpawp(i%t)) THEN ifpaw: IF (upf(i%t)%tpawp) THEN
! !
ALLOCATE(rho_lm(i%m,lmaxq**2,nspin)) ALLOCATE(rho_lm(i%m,lmaxq**2,nspin))
! !
whattodo: DO i_what = AE, PS whattodo: DO i_what = AE, PS
i%w = i_what
! !
IF (i%w == AE) THEN IF (i_what == AE) THEN
CALL PAW_rho_lm(i, becsum, pfunc, rho_lm) CALL PAW_rho_lm(i, becsum, pfunc, rho_lm)
i_sign = +1._dp i_sign = +1._dp
ELSE ELSE
CALL PAW_rho_lm(i, becsum, ptfunc, rho_lm, aug) CALL PAW_rho_lm(i, becsum, ptfunc, rho_lm, upf(i%t)%paw%aug)
i_sign = -1._dp i_sign = -1._dp
ENDIF ENDIF
! !
@ -520,7 +522,7 @@ FUNCTION PAW_integrate(becsum)
DO lm = 1, lmaxq**2 DO lm = 1, lmaxq**2
! I can use rho_lm itself as workspace ! I can use rho_lm itself as workspace
rho_lm(1:i%m,lm,is) = rho_lm(1:i%m,lm,is) & rho_lm(1:i%m,lm,is) = rho_lm(1:i%m,lm,is) &
* saved(i%a)%v(1:i%m,lm,is,i%w) * saved(i%a)%v(1:i%m,lm,is,i_what)
CALL simpson (i%m,rho_lm(:,lm,is),g(i%t)%rab,integral) CALL simpson (i%m,rho_lm(:,lm,is),g(i%t)%rab,integral)
!WRITE(6,"(3i3,10f20.10)") i_what, is, lm, integral !WRITE(6,"(3i3,10f20.10)") i_what, is, lm, integral
PAW_integrate = PAW_integrate + i_sign*integral PAW_integrate = PAW_integrate + i_sign*integral
@ -552,8 +554,8 @@ FUNCTION PAW_ddot(bec1,bec2)
USE lsda_mod, ONLY : nspin USE lsda_mod, ONLY : nspin
USE ions_base, ONLY : nat, ityp USE ions_base, ONLY : nat, ityp
USE atom, ONLY : g => rgrid USE atom, ONLY : g => rgrid
USE paw_variables, ONLY : pfunc, ptfunc, tpawp, aug USE paw_variables, ONLY : pfunc, ptfunc
USE uspp_param, ONLY : nhm, lmaxq USE uspp_param, ONLY : nhm, lmaxq, upf
REAL(DP) :: PAW_ddot REAL(DP) :: PAW_ddot
@ -586,8 +588,10 @@ FUNCTION PAW_ddot(bec1,bec2)
i%a = na ! the index of the atom i%a = na ! the index of the atom
i%t = ityp(na) ! the type of atom na i%t = ityp(na) ! the type of atom na
i%m = g(i%t)%mesh ! radial mesh size for atom na i%m = g(i%t)%mesh ! radial mesh size for atom na
i%w = upf(i%t)%nwfc
i%l = upf(i%t)%nqlc-1
! !
ifpaw: IF (tpawp(i%t)) THEN ifpaw: IF (upf(i%t)%tpawp) THEN
! !
ALLOCATE(rho_lm(i%m,lmaxq**2,nspin)) ALLOCATE(rho_lm(i%m,lmaxq**2,nspin))
ALLOCATE(v_lm(i%m,lmaxq**2)) ALLOCATE(v_lm(i%m,lmaxq**2))
@ -598,7 +602,7 @@ FUNCTION PAW_ddot(bec1,bec2)
CALL PAW_rho_lm(i, bec1, pfunc, rho_lm) CALL PAW_rho_lm(i, bec1, pfunc, rho_lm)
i_sign = +1._dp i_sign = +1._dp
ELSE ELSE
CALL PAW_rho_lm(i, bec1, ptfunc, rho_lm, aug) !fun) CALL PAW_rho_lm(i, bec1, ptfunc, rho_lm, upf(i%t)%paw%aug)
i_sign = -1._dp i_sign = -1._dp
ENDIF ENDIF
! !
@ -609,7 +613,7 @@ FUNCTION PAW_ddot(bec1,bec2)
IF (i_what == AE) THEN IF (i_what == AE) THEN
CALL PAW_rho_lm(i, bec2, pfunc, rho_lm) CALL PAW_rho_lm(i, bec2, pfunc, rho_lm)
ELSE ELSE
CALL PAW_rho_lm(i, bec2, ptfunc, rho_lm, aug) !fun) CALL PAW_rho_lm(i, bec2, ptfunc, rho_lm, upf(i%t)%paw%aug)
ENDIF ENDIF
! !
! Finally compute the integral ! Finally compute the integral
@ -1132,15 +1136,15 @@ SUBROUTINE PAW_rho_lm(i, becsum, pfunc, rho_lm, aug)
USE parameters, ONLY : lqmax USE parameters, ONLY : lqmax
USE constants, ONLY : eps12 USE constants, ONLY : eps12
USE radial_grids, ONLY : ndmx USE radial_grids, ONLY : ndmx
USE paw_variables, ONLY : augfun_t, nbrx USE paw_variables, ONLY : nbrx
USE atom, ONLY : g => rgrid USE atom, ONLY : g => rgrid
TYPE(paw_info) :: i ! atom's minimal info TYPE(paw_info) :: i ! atom's minimal info
REAL(DP), INTENT(IN) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupation REAL(DP), INTENT(IN) :: becsum(nhm*(nhm+1)/2,nat,nspin)! cross band occupation
REAL(DP), INTENT(IN) :: pfunc(ndmx,nbrx,nbrx,ntyp) ! psi_i * psi_j REAL(DP), INTENT(IN) :: pfunc(ndmx,nbrx,nbrx,ntyp) ! psi_i * psi_j
REAL(DP), INTENT(OUT) :: rho_lm(i%m,lmaxq**2,nspin) ! AE charge density on rad. grid REAL(DP), INTENT(OUT) :: rho_lm(i%m,lmaxq**2,nspin) ! AE charge density on rad. grid
TYPE(augfun_t), OPTIONAL,INTENT(IN) :: & REAL(DP), OPTIONAL,INTENT(IN) :: &
aug(ntyp) ! augmentation functions (only for PS part) aug(i%m,i%w,i%w,0:2*i%l) ! augmentation functions (only for PS part)
REAL(DP) :: pref ! workspace (ap*becsum) REAL(DP) :: pref ! workspace (ap*becsum)
@ -1188,7 +1192,7 @@ SUBROUTINE PAW_rho_lm(i, becsum, pfunc, rho_lm, aug)
! if I'm doing the pseudo part I have to add the augmentation charge ! if I'm doing the pseudo part I have to add the augmentation charge
l = INT(SQRT(DBLE(lm-1))) ! l has to start from zero l = INT(SQRT(DBLE(lm-1))) ! l has to start from zero
rho_lm(1:i%m,lm,ispin) = rho_lm(1:i%m,lm,ispin) +& rho_lm(1:i%m,lm,ispin) = rho_lm(1:i%m,lm,ispin) +&
pref * aug(i%t)%fun(1:i%m, indv(nb,i%t), indv(mb,i%t), l) pref * aug(1:i%m, indv(nb,i%t), indv(mb,i%t), l)
ENDIF ! augfun ENDIF ! augfun
ENDDO angular_momentum ENDDO angular_momentum
ENDDO !mb ENDDO !mb

4
PW/read_file.f90
View File

@ -42,7 +42,7 @@ SUBROUTINE read_file()
USE pw_restart, ONLY : pw_readfile USE pw_restart, ONLY : pw_readfile
USE xml_io_base, ONLY : pp_check_file USE xml_io_base, ONLY : pp_check_file
USE uspp, ONLY : okvan USE uspp, ONLY : okvan
USE paw_variables, ONLY : okpaw, tpawp USE paw_variables, ONLY : okpaw
USE paw_init, ONLY : paw_init_onecenter, allocate_paw_internals USE paw_init, ONLY : paw_init_onecenter, allocate_paw_internals
USE ldaU, ONLY : v_hub, eth USE ldaU, ONLY : v_hub, eth
! !
@ -153,7 +153,7 @@ SUBROUTINE read_file()
CALL readpp() CALL readpp()
! !
okvan = ANY ( upf(:)%tvanp ) okvan = ANY ( upf(:)%tvanp )
okpaw = ANY ( tpawp(1:nsp) ) okpaw = ANY ( upf(1:nsp)%tpawp )
! !
! ... check for spin-orbit pseudopotentials ! ... check for spin-orbit pseudopotentials
! !

2
PW/read_pseudo.f90
View File

@ -27,7 +27,6 @@ subroutine readpp
USE io_global, ONLY : stdout USE io_global, ONLY : stdout
USE ions_base, ONLY : zv USE ions_base, ONLY : zv
USE uspp_param, ONLY : upf USE uspp_param, ONLY : upf
USE paw_variables, ONLY : tpawp
USE read_paw_module, ONLY : paw_io, allocate_pseudo_paw, deallocate_pseudo_paw USE read_paw_module, ONLY : paw_io, allocate_pseudo_paw, deallocate_pseudo_paw
USE paw_to_internal, ONLY : set_pseudo_paw USE paw_to_internal, ONLY : set_pseudo_paw
implicit none implicit none
@ -55,7 +54,6 @@ subroutine readpp
END IF END IF
ALLOCATE ( upf(ntyp) ) ALLOCATE ( upf(ntyp) )
do nt = 1, ntyp do nt = 1, ntyp
!tpawp(nt) = .false.
! !
! obsolescent variables, not read from UPF format, no longer used ! obsolescent variables, not read from UPF format, no longer used
! !

4
PW/setup.f90
View File

@ -86,7 +86,7 @@ SUBROUTINE setup()
USE exx, ONLY : exx_grid_init USE exx, ONLY : exx_grid_init
#endif #endif
USE funct, ONLY : dft_is_meta, dft_is_hybrid, dft_is_gradient USE funct, ONLY : dft_is_meta, dft_is_hybrid, dft_is_gradient
USE paw_variables, ONLY : okpaw, tpawp USE paw_variables, ONLY : okpaw
! !
IMPLICIT NONE IMPLICIT NONE
! !
@ -698,7 +698,7 @@ SUBROUTINE setup()
! ... okvan = .TRUE. : at least one pseudopotential is US ! ... okvan = .TRUE. : at least one pseudopotential is US
! !
okvan = ANY( upf(:)%tvanp ) okvan = ANY( upf(:)%tvanp )
okpaw = ANY( tpawp(1:ntyp) ) okpaw = ANY( upf(1:ntyp)%tpawp )
! !
! ... Needed for LDA+U and PAW ! ... Needed for LDA+U and PAW
! !

3
PW/summary.f90
View File

@ -406,7 +406,6 @@ SUBROUTINE print_ps_info
USE ions_base, ONLY : ntyp => nsp USE ions_base, ONLY : ntyp => nsp
USE atom, ONLY : rgrid USE atom, ONLY : rgrid
USE uspp_param, ONLY : upf USE uspp_param, ONLY : upf
USE paw_variables, ONLY: tpawp
! !
INTEGER :: nt INTEGER :: nt
CHARACTER :: ps*35 CHARACTER :: ps*35
@ -415,7 +414,7 @@ SUBROUTINE print_ps_info
! !
IF ( upf(nt)%tvanp ) THEN IF ( upf(nt)%tvanp ) THEN
ps='Ultrasoft' ps='Ultrasoft'
ELSE IF ( tpawp (nt) ) THEN ELSE IF ( upf(nt)%tpawp ) THEN
ps="Projector augmented-wave" ps="Projector augmented-wave"
ELSE ELSE
ps='Norm-conserving' ps='Norm-conserving'