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

Improved generation of FR PAW dataset. 

Modification to the UPF format: the small components of the all-electron
partial waves exported and imported in the FR case.


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@6487 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
dalcorso 2010-03-13 10:45:42 +00:00
parent c816f3462c
commit bac79d8452
12 changed files with 233 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Modules/pseudo_types.f90
View File

@ -22,7 +22,9 @@
REAL(DP),POINTER :: ae_rho_atc(:) ! AE core charge (pseudo ccharge
! is already included in upf)
REAL(DP),POINTER :: pfunc(:,:,:),&! Psi_i(r)*Psi_j(r)
ptfunc(:,:,:) ! as above, but for pseudo
pfunc_rel(:,:,:), & ! Psi_i(r)*Psi_j(r) small component
ptfunc(:,:,:), & ! as above, but for pseudo
aewfc_rel(:,:) ! as above, but for pseudo
REAL(DP),POINTER :: ae_vloc(:) ! AE local potential (pseudo vloc
! is already included in upf)
REAL(DP),POINTER :: oc(:) ! starting occupation used to init becsum
@ -162,7 +164,9 @@
SUBROUTINE nullify_paw_in_upf( paw )
TYPE( paw_in_upf ), INTENT(INOUT) :: paw
NULLIFY( paw%ae_rho_atc )
NULLIFY( paw%aewfc_rel )
NULLIFY( paw%pfunc )
NULLIFY( paw%pfunc_rel )
NULLIFY( paw%ptfunc )
NULLIFY( paw%ae_vloc )
NULLIFY( paw%augmom )
@ -172,7 +176,9 @@
SUBROUTINE deallocate_paw_in_upf( paw )
TYPE( paw_in_upf ), INTENT(INOUT) :: paw
IF( ASSOCIATED( paw%ae_rho_atc ) ) DEALLOCATE ( paw%ae_rho_atc )
IF( ASSOCIATED( paw%aewfc_rel ) ) DEALLOCATE (paw%aewfc_rel )
IF( ASSOCIATED( paw%pfunc ) ) DEALLOCATE ( paw%pfunc )
IF( ASSOCIATED( paw%pfunc_rel ) ) DEALLOCATE ( paw%pfunc_rel )
IF( ASSOCIATED( paw%ptfunc ) ) DEALLOCATE ( paw%ptfunc )
IF( ASSOCIATED( paw%ae_vloc ) ) DEALLOCATE ( paw%ae_vloc )
IF( ASSOCIATED( paw%augmom ) ) DEALLOCATE ( paw%augmom )

32
Modules/read_upf_v2.f90
View File

@ -344,6 +344,7 @@ SUBROUTINE read_upf_v2(u, upf, grid, ierr) !
DO l = abs(ln-lm),ln+lm,2 ! only even terms
CALL iotk_scan_dat(u, 'PP_QIJL'//iotk_index((/nb,mb,l/)),&
upf%qfuncl(:,nmb,l),default=zeros,attr=attr)
IF( upf%tpawp) upf%qfuncl(upf%paw%iraug+1:,nmb,l) = 0._dp
ENDDO
ELSE q_with_l
CALL iotk_scan_dat(u, 'PP_QIJ'//iotk_index((/nb,mb/)),&
@ -412,6 +413,7 @@ SUBROUTINE read_upf_v2(u, upf, grid, ierr) !
TYPE(pseudo_upf),INTENT(INOUT) :: upf ! the pseudo data
INTEGER :: ierr ! /= 0 if something went wrong
CHARACTER(len=iotk_attlenx) :: attr
LOGICAL :: exst
!
INTEGER :: nb
!
@ -425,6 +427,18 @@ SUBROUTINE read_upf_v2(u, upf, grid, ierr) !
upf%aewfc(:,nb), attr=attr)
ENDDO
IF (upf%has_so .and. upf%tpawp) THEN
ALLOCATE( upf%paw%aewfc_rel(upf%mesh, upf%nbeta) )
nb_loop: DO nb = 1,upf%nbeta
CALL iotk_scan_dat(u, 'PP_AEWFC_REL'//iotk_index(nb), &
upf%paw%aewfc_rel(:,nb), attr=attr, found=exst)
IF (.not.exst) THEN
upf%paw%aewfc_rel=0.0_DP
EXIT nb_loop
ENDIF
ENDDO nb_loop
ENDIF
ALLOCATE( upf%pswfc(upf%mesh, upf%nbeta) )
DO nb = 1,upf%nbeta
CALL iotk_scan_dat(u, 'PP_PSWFC'//iotk_index(nb), &
@ -510,13 +524,31 @@ SUBROUTINE read_upf_v2(u, upf, grid, ierr) !
!
ALLOCATE(upf%paw%pfunc(upf%mesh, upf%nbeta,upf%nbeta) )
upf%paw%pfunc(:,:,:) = 0._dp
IF (upf%has_so) THEN
ALLOCATE(upf%paw%pfunc_rel(upf%mesh, upf%nbeta,upf%nbeta) )
upf%paw%pfunc_rel(:,:,:) = 0._dp
ENDIF
DO nb=1,upf%nbeta
DO nb1=1,nb
upf%paw%pfunc (1:upf%mesh, nb, nb1) = &
upf%aewfc(1:upf%mesh, nb) * upf%aewfc(1:upf%mesh, nb1)
IF (upf%has_so) THEN
upf%paw%pfunc_rel (1:upf%paw%iraug, nb, nb1) = &
upf%paw%aewfc_rel(1:upf%paw%iraug, nb) * &
upf%paw%aewfc_rel(1:upf%paw%iraug, nb1)
!
! The small component is added to pfunc. pfunc_rel is useful only
! to add a small magnetic contribution
!
upf%paw%pfunc (1:upf%paw%iraug, nb, nb1) = &
upf%paw%pfunc (1:upf%paw%iraug, nb, nb1) + &
upf%paw%pfunc_rel (1:upf%paw%iraug, nb, nb1)
ENDIF
upf%paw%pfunc(upf%paw%iraug+1:,nb,nb1) = 0._dp
!
upf%paw%pfunc (1:upf%mesh, nb1, nb) = upf%paw%pfunc (1:upf%mesh, nb, nb1)
IF (upf%has_so) upf%paw%pfunc_rel (1:upf%mesh, nb1, nb) = &
upf%paw%pfunc_rel (1:upf%mesh, nb, nb1)
ENDDO
ENDDO
!

10
Modules/write_upf_v2.f90
View File

@ -479,6 +479,16 @@ SUBROUTINE write_upf_v2(u, upf, conf) !
CALL iotk_write_dat(u, 'PP_AEWFC'//iotk_index(nb), &
upf%aewfc(:,nb), columns=4, attr=attr)
ENDDO
IF (upf%has_so.and.upf%tpawp) THEN
DO nb = 1,upf%nbeta
CALL iotk_write_attr(attr, 'index', nb, first=.true.)
CALL iotk_write_attr(attr, 'label', upf%els_beta(nb))
CALL iotk_write_attr(attr, 'l', upf%lll(nb))
CALL iotk_write_attr(attr, 'j', upf%jjj(nb))
CALL iotk_write_dat(u, 'PP_AEWFC_REL'//iotk_index(nb), &
upf%paw%aewfc_rel(:,nb), columns=4, attr=attr)
ENDDO
ENDIF
! Pseudo wavefunctions
DO nb = 1,upf%nbeta
CALL iotk_write_attr(attr, 'index', nb, first=.true.)

41
atomic/atomic_paw.f90
View File

@ -29,6 +29,7 @@ MODULE atomic_paw
!
! Written by Guido Fratesi, february 2005
! Modified by Riccardo Mazzarello, july 2006
! Fully Relativistic generalization by Andrea Dal Corso, november 2009
! Other people involved: Lorenzo Paulatto and Stefano de Gironcoli
!
!============================================================================
@ -104,8 +105,9 @@ CONTAINS
!
do is=1,nspin_
do n=2,pawset_%grid%mesh
!if (chargeps(n,is)<-1.d-12) &
! call errore('new_paw_hamiltonian','negative rho',1)
! write(6,*) n, pawset_%grid%r(n), chargeps(n,is)
if (chargeps(n,is)<-1.d-12) &
call errore('new_paw_hamiltonian','negative rho',1)
enddo
enddo
@ -161,9 +163,9 @@ CONTAINS
!
SUBROUTINE us2paw (pawset_, &
zval, grid, rmatch_augfun, ikk, &
nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, phis, betas, &
qvan, kindiff, &
nlcc, aerhoc, psrhoc, aevtot, psvtot, which_paw_augfun )
nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, psipaw_rel, &
phis, betas, qvan, kindiff, &
nlcc, aerhoc, psrhoc, aevtot, psvtot, which_paw_augfun,rel )
USE funct, ONLY : dft_name, get_iexch, get_icorr, get_igcx, get_igcc
USE ld1inc, ONLY : zed, file_screen
@ -178,7 +180,7 @@ CONTAINS
REAL(dp), INTENT(IN) :: rmatch_augfun
INTEGER, INTENT(IN) :: ikk(nwfsx)
!
INTEGER, INTENT(IN) :: nbeta
INTEGER, INTENT(IN) :: nbeta, rel
INTEGER, INTENT(IN) :: lls(nwfsx)
REAL(dp), INTENT(IN) :: ocs(nwfsx)
CHARACTER(LEN=2), INTENT(IN) :: els(nwfsx)
@ -186,6 +188,7 @@ CONTAINS
REAL(dp), INTENT(IN) :: rcutus(nwfsx)
REAL(dp), INTENT(IN) :: enls(nwfsx)
REAL(dp), INTENT(IN) :: psipaw(ndmx,nwfsx)
REAL(dp), INTENT(IN) :: psipaw_rel(ndmx,nwfsx)
REAL(dp), INTENT(IN) :: phis(ndmx,nwfsx)
REAL(dp), INTENT(IN) :: betas(ndmx,nwfsx)
!
@ -241,6 +244,7 @@ CONTAINS
!
pawset_%rmatch_augfun = rmatch_augfun
if (rmatch_augfun <= 0.0_dp) pawset_%rmatch_augfun = grid%r(irc)
pawset_%rel = rel
pawset_%irc = irc
pawset_%ikk(1:nbeta)=ikk(1:nbeta)
!
@ -254,6 +258,7 @@ CONTAINS
pawset_%els(1:nbeta)= els(1:nbeta)
pawset_%enl(1:nbeta)= enls(1:nbeta)
pawset_%aewfc(1:mesh,1:nbeta) = psipaw(1:mesh,1:nbeta)
pawset_%aewfc_rel(1:mesh,1:nbeta) = psipaw_rel(1:mesh,1:nbeta)
pawset_%pswfc(1:mesh,1:nbeta) = phis (1:mesh,1:nbeta)
pawset_%proj (1:mesh,1:nbeta) = betas (1:mesh,1:nbeta)
!
@ -290,6 +295,9 @@ CONTAINS
pawset_%augfun(1:mesh,ns,ns1,l3) = &
pawset_%aewfc(1:mesh,ns) * pawset_%aewfc(1:mesh,ns1) - &
pawset_%pswfc(1:mesh,ns) * pawset_%pswfc(1:mesh,ns1)
IF (pawset_%rel==2) &
pawset_%augfun(1:irc,ns,ns1,l3) =pawset_%augfun(1:irc,ns,ns1,l3)&
+pawset_%aewfc_rel(1:irc,ns) * pawset_%aewfc_rel(1:irc,ns1)
pawset_%augfun(1:mesh,ns1,ns,l3) = pawset_%augfun(1:mesh,ns,ns1,l3)
aux(1:irc) = pawset_%augfun(1:irc,ns,ns1,l3) * pawset_%grid%r(1:irc)**l3
lll = l1 + l2 + 2 + l3
@ -462,6 +470,9 @@ CONTAINS
aug_real(1:mesh,ns,ns1) = &
pawset_%aewfc(1:mesh,ns) * pawset_%aewfc(1:mesh,ns1) - &
pawset_%pswfc(1:mesh,ns) * pawset_%pswfc(1:mesh,ns1)
IF (pawset_%rel==2) &
aug_real(1:irc,ns,ns1) = aug_real(1:irc,ns,ns1) + &
pawset_%aewfc_rel(1:irc,ns)*pawset_%aewfc_rel(1:irc,ns1)
aug_real(1:mesh,ns1,ns) = aug_real(1:mesh,ns,ns1)
ENDDO
ENDDO
@ -541,7 +552,7 @@ CONTAINS
! ...
!
SUBROUTINE paw2us (pawset_,zval,grid,nbeta,lls,jjs,ikk,betas,qq,qvan,&
vpsloc, bmat, rhos, els, rcutus, pseudotype)
vpsloc, bmat, rhos, els, rcutus, pseudotype,psipaw_rel)
USE funct, ONLY : set_dft_from_name
use radial_grids, only: radial_grid_type, allocate_radial_grid
IMPLICIT NONE
@ -556,6 +567,7 @@ CONTAINS
REAL(dp), INTENT(OUT) :: jjs(nwfsx)
REAL(dp), INTENT(OUT) :: rcutus(nwfsx)
REAL(dp), INTENT(OUT) :: betas(ndmx,nwfsx)
REAL(dp), INTENT(OUT) :: psipaw_rel(ndmx,nwfsx)
REAL(dp), INTENT(OUT) :: qq(nwfsx,nwfsx)
REAL(dp), INTENT(OUT) :: qvan(ndmx,nwfsx,nwfsx)
REAL(dp), INTENT(OUT) :: vpsloc(ndmx) ! the local pseudopotential
@ -602,6 +614,7 @@ CONTAINS
!
rhos(1:mesh)=pawset_%pscharge(1:mesh)
!
psipaw_rel(1:mesh,1:nbeta)=pawset_%aewfc_rel(1:mesh,1:nbeta)
betas(1:mesh,1:nbeta)=pawset_%proj(1:mesh,1:nbeta)
pseudotype=3
!
@ -903,6 +916,11 @@ CONTAINS
pawset_%aewfc(1:pawset_%grid%mesh,ns ) * &
pawset_%aewfc(1:pawset_%grid%mesh,ns1) * &
veff1_(1:pawset_%grid%mesh,is)
IF (pawset_%rel==2) &
aux(1:pawset_%irc) = aux(1:pawset_%irc) + &
pawset_%aewfc_rel(1:pawset_%irc,ns ) * &
pawset_%aewfc_rel(1:pawset_%irc,ns1) * &
veff1_(1:pawset_%irc,is)
dd = dd + &
int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
! dddd(ns,ns1,2) = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
@ -957,6 +975,11 @@ CONTAINS
pawset_%aewfc(1:pawset_%grid%mesh,ns ) * &
pawset_%aewfc(1:pawset_%grid%mesh,ns1) * &
pawset_%aeloc(1:pawset_%grid%mesh)
IF (pawset_%rel==2) &
aux(1:pawset_%irc) = aux(1:pawset_%irc) + &
pawset_%aewfc_rel(1:pawset_%irc,ns ) * &
pawset_%aewfc_rel(1:pawset_%irc,ns1) * &
pawset_%aeloc(1:pawset_%irc)
dd = int_0_inf_dr(aux,pawset_%grid,pawset_%irc,(pawset_%l(ns)+1)*2)
! Int[ps*v1~*ps + Q*v1~]
aux(1:pawset_%grid%mesh) = &
@ -1126,6 +1149,10 @@ CONTAINS
charge1_(1:pawset_%grid%mesh,is) = charge1_(1:pawset_%grid%mesh,is) + factor * &
projsum_(ns,ns1,is) * pawset_%aewfc(1:pawset_%grid%mesh,ns) * &
pawset_%aewfc(1:pawset_%grid%mesh,ns1)
IF (pawset_%rel==2) &
charge1_(1:pawset_%irc,is) = charge1_(1:pawset_%irc,is) + factor * &
projsum_(ns,ns1,is) * pawset_%aewfc_rel(1:pawset_%irc,ns) * &
pawset_%aewfc_rel(1:pawset_%irc,ns1)
CASE ("PS")
charge1_(1:pawset_%grid%mesh,is) = charge1_(1:pawset_%grid%mesh,is) + factor * &
projsum_(ns,ns1,is) * pawset_%pswfc(1:pawset_%grid%mesh,ns) * &

6
atomic/export_upf.f90
View File

@ -286,6 +286,12 @@ SUBROUTINE export_upf(iunps)
= pawsetup%augmom(1:upf%nbeta,1:upf%nbeta,0:2*upf%lmax)
!
upf%kkbeta = max(upf%kkbeta, upf%paw%iraug)
IF (upf%has_so) THEN
ALLOCATE( upf%paw%aewfc_rel(upf%mesh, upf%nbeta) )
upf%paw%aewfc_rel(1:upf%mesh,1:upf%nbeta) = &
pawsetup%aewfc_rel(1:upf%mesh,1:upf%nbeta)
ENDIF
!
!upf%paw%pfunc(:) = not used when writing, reconstructed from upf%aewfc
!upf%paw%ptfunc(:) = not used when writing, reconstructed from upf%pswfc

189
atomic/gener_pseudo.f90
View File

@ -41,25 +41,26 @@ subroutine gener_pseudo
qvan, qvanl, qq, bmat, ddd, betas, nbeta, ikk, pseudotype, &
pawsetup, zval, vpsloc, vpot, vnl, lpaw, rcloc, rcutus, &
enl, enls, rcut, chis, nstoae, rmatch_augfun,&
lnc2paw, rcutnc2paw, rhos, which_augfun
lnc2paw, rcutnc2paw, rhos, which_augfun, psipaw_rel, &
cau_fact, ik, ikus
use atomic_paw, only : us2paw, paw2us
implicit none
integer :: &
ik, & ! the point corresponding to rc
ikus, & ! the point corresponding to rc ultrasoft
ikloc, & ! the point corresponding to rc local
ns, & ! counter on pseudo functions
ns1, & ! counter on pseudo functions
nnode, & ! the number of nodes of phi
lam ! the angular momentum
lam, & ! the angular momentum
irc ! the maximum radius of the sphere
real(DP) :: &
xc(8), & ! parameters of bessel functions
xc(8), & ! parameters of bessel functions
psi_in(ndmx), & ! the all_electron wavefunction
gi(ndmx,2), & ! auxiliary to compute the integrals
occ, norm1, &
db, work(nwfsx) ! work space
gi(ndmx), & ! auxiliary to compute the integrals
occ, norm1, norm2, & !
speed, & ! an estimate of electron speed at rcutus
db, work(nwfsx) ! work space
real(DP), allocatable :: &
b(:,:), binv(:,:) ! the B matrix and its inverse
@ -89,8 +90,6 @@ subroutine gener_pseudo
! generation (normally the AE potential)
integer :: iknc2paw ! point in rgrid closer to rcutnc2paw
real(DP) :: q, fac, pi, wrk(ndmx), jlq(ndmx), norm(nwfsx), normr(nwfsx)
if (lpaw) then
write(stdout, &
'(/,5x,21(''-''),'' Generating PAW atomic setup '',20(''-''),/)')
@ -135,31 +134,79 @@ subroutine gener_pseudo
if (enls(n) == 0.0_dp) enls(n)=enl(nstoae(n))
enddo
!
ik=0
ikus=0
ikloc=0
do ns=1,nbeta
do n=1,grid%mesh
if (grid%r(n).lt.rcut(ns)) ik(ns)=n
if (grid%r(n).lt.rcutus(ns)) ikus(ns)=n
if (grid%r(n).lt.rcloc) ikloc=n
enddo
if (mod(ik(ns),2) == 0) ik(ns)=ik(ns)+1
if (mod(ikus(ns),2) == 0) ikus(ns)=ikus(ns)+1
if (mod(ikloc,2) == 0) ikloc=ikloc+1
if (ik(ns).gt.grid%mesh) call errore('gener_pseudo','ik is wrong ',ns)
if (ikus(ns)+10.gt.grid%mesh) call errore('gener_pseudo','ikus is wrong ',ns)
if (ikloc+5.gt.grid%mesh) call errore('gener_pseudo','ikloc is wrong ',ns)
if (pseudotype == 3) then
ikk(ns)=max(ikus(ns)+10,ikloc+5)
else
ikk(ns)=max(ik(ns)+10,ikloc+5)
endif
if (ikk(ns).gt.grid%mesh) call errore('gener_pseudo','ikk is wrong ',ns)
end do
do ns=1,nbeta
do ns1=1,nbeta
if (lls(ns) == lls(ns1).and.ikk(ns1).gt.ikk(ns)) &
ikk(ns)=ikk(ns1)
enddo
enddo
irc = maxval(ikk(1:nbeta))+8
IF (mod(irc,2) == 0) irc=irc+1
IF (irc>grid%mesh) CALL errore('gener_pseudo','irc is too large',1)
!
! Set the all-electron wavefunctions, calculating those at user supplied
! energies. The wavefunctions are written on file at this point, so
! the user can check them also when the pseudopotential generation is
! unsuccessful
!
IF (rel==2) WRITE(stdout,'(/,5x,"Fully relativistic calculation: |psi|^2")')
do ns=1,nbeta
ik=0
nwf0=nstoae(ns)
do n=1,grid%mesh
if (grid%r(n).lt.rcut(ns)) ik=n
enddo
if (mod(ik,2) == 0) ik=ik+1
if (new(ns)) then
call set_psi_in(ik,lls(ns),jjs(ns),enls(ns),psipaw(1,ns))
call set_psi_in(ik(ns),lls(ns),jjs(ns),enls(ns),psipaw(1,ns),&
psipaw_rel(1,ns))
else
lam=lls(ns)
nst=(lam+1)*2
psipaw(:,ns)=psi(:,1,nwf0)
do n=1,grid%mesh
gi(n,1)=psipaw(n,ns)*psipaw(n,ns)
gi(n)=psipaw(n,ns)*psipaw(n,ns)
enddo
norm1=sqrt(int_0_inf_dr(gi,grid,grid%mesh,nst))
IF (rel==2) THEN
psipaw_rel(:,ns)=psi(:,2,nwf0)
DO n=1,irc
gi(n)=psipaw_rel(n,ns)*psipaw_rel(n,ns)
ENDDO
norm2=sqrt(int_0_inf_dr(gi,grid,irc,nst))
WRITE(stdout,'(/,5x,"Wfc ",a2," LC norm =",f12.8," SC norm ="&
&,f12.8," missing",f12.8)') els(ns), norm1**2, &
norm2**2, 1.0_DP-(norm1**2+norm2**2)
norm1=sqrt(norm1**2+norm2**2)
psipaw_rel(:,ns)=psipaw_rel(:,ns)/norm1
IF (lpaw) THEN
speed=psipaw(ikus(ns),ns)*(enls(ns)-vpotpaw(ikus(ns))) &
*psipaw(ikus(ns),ns)
WRITE(stdout,'(5x,"Speed at rcutus ",f10.3, " a.u., &
&v/c =",f10.6,", (v/c)^2 =",f12.8 )') sqrt(abs(speed)), &
sqrt(abs(speed)/cau_fact**2),speed/cau_fact**2
END IF
END IF
psipaw(:,ns)=psipaw(:,ns)/norm1
endif
enddo
END IF
END DO
call write_wfcfile(file_wfcaegen,psipaw,els,nwfs)
!
@ -175,29 +222,12 @@ subroutine gener_pseudo
!
! compute the ik closer to r_cut, r_cutus, rcloc
!
ik=0
ikus=0
ikloc=0
do n=1,grid%mesh
if (grid%r(n).lt.rcut(ns)) ik=n
if (grid%r(n).lt.rcutus(ns)) ikus=n
if (grid%r(n).lt.rcloc) ikloc=n
enddo
if (mod(ik,2) == 0) ik=ik+1
if (mod(ikus,2) == 0) ikus=ikus+1
if (mod(ikloc,2) == 0) ikloc=ikloc+1
if (lnc2paw) then
do n=1,grid%mesh
if (grid%r(n).lt.rcutnc2paw(ns)) iknc2paw=n
end do
if (mod(iknc2paw,2) == 0) iknc2paw=iknc2paw+1
end if
if (ikus.gt.grid%mesh) call errore('gener_pseudo','ik is wrong ',1)
if (pseudotype == 3) then
ikk(ns)=max(ikus+10,ikloc+5)
else
ikk(ns)=max(ik+10,ikloc+5)
endif
if (new(ns)) then
occ=1.0_DP
@ -219,13 +249,15 @@ subroutine gener_pseudo
psipaw(1:grid%mesh,ns)=phis(1:grid%mesh,ns)
endif
!
IF (which_augfun=='PSQ' .and. .not. lpaw) THEN
! IF (which_augfun=='PSQ' .and. .not. lpaw) THEN
IF ((which_augfun=='PSQ'.AND.ik(ns)/=ikus(ns)).OR.&
(lpaw.AND..NOT.lnc2paw)) THEN
psipsus(:,ns)=psi_in(:)
ELSE
if (tm) then
call compute_phi_tm(lam,ik,psi_in,phis(1,ns),1,xc,enls(ns),els(ns))
call compute_phi_tm(lam,ik(ns),psi_in,phis(1,ns),1,xc,enls(ns),els(ns))
else
call compute_phi(lam,ik,psi_in,phis(1,ns),xc,1,occ,enls(ns),els(ns))
call compute_phi(lam,ik(ns),psi_in,phis(1,ns),xc,1,occ,enls(ns),els(ns))
ecutrho=max(ecutrho,8.0_dp*xc(6)**2)
endif
!
@ -233,16 +265,16 @@ subroutine gener_pseudo
!
psipsus(:,ns)=phis(:,ns)
ENDIF
if (ikus.ne.ik) then
call compute_phius(lam,ikus,psipsus(1,ns),phis(1,ns),xc,1,els(ns))
if (ikus(ns).ne.ik(ns)) then
call compute_phius(lam,ikus(ns),psipsus(1,ns),phis(1,ns),xc,1,els(ns))
ecutwfc=max(ecutwfc,2.0_dp*xc(5)**2)
lbes4=.true.
else
lbes4=.false.
if (.not.tm) ecutwfc=max(ecutwfc,2.0_dp*xc(6)**2)
endif
if (tm.and.ik==ikus) then
call compute_chi_tm(lam,ik,ikk(ns),phis(1,ns),chis(1,ns),xc,enls(ns))
if (tm.and.ik(ns)==ikus(ns)) then
call compute_chi_tm(lam,ik(ns),ikk(ns),phis(1,ns),chis(1,ns),xc,enls(ns))
else
call compute_chi(lam,ikk(ns),phis(1,ns),chis(1,ns),xc,enls(ns),lbes4)
endif
@ -250,12 +282,6 @@ subroutine gener_pseudo
!
! for each angular momentum take the same integration point
!
do ns=1,nbeta
do ns1=1,nbeta
if (lls(ns) == lls(ns1).and.ikk(ns1).gt.ikk(ns)) &
ikk(ns)=ikk(ns1)
enddo
enddo
!
! construct B_{ij}
!
@ -266,7 +292,7 @@ subroutine gener_pseudo
nst=(lls(ns)+1)*2
ikl=ikk(ns1)
do n=1,grid%mesh
gi(n,1)=phis(n,ns)*chis(n,ns1)
gi(n)=phis(n,ns)*chis(n,ns1)
enddo
bmat(ns,ns1)=int_0_inf_dr(gi,grid,ikl,nst)
endif
@ -340,11 +366,21 @@ subroutine gener_pseudo
do ns=1,nbeta
do ns1=1,ns
ikl=max(ikk(ns),ikk(ns1))
do n=1, ikl
qvan(n,ns,ns1) = psipsus(n,ns) * psipsus(n,ns1) &
- phis(n,ns) * phis(n,ns1)
gi(n,1)=qvan(n,ns,ns1)
enddo
if (which_augfun=='PSQ'.and.rel==2) then
do n=1, ikl
qvan(n,ns,ns1) = psipsus(n,ns) * psipsus(n,ns1) &
+ psipaw_rel(n,ns) * psipaw_rel(n,ns1) &
- phis(n,ns) * phis(n,ns1)
gi(n)=qvan(n,ns,ns1)
enddo
else
do n=1, ikl
qvan(n,ns,ns1) = psipsus(n,ns) * psipsus(n,ns1) &
- phis(n,ns) * phis(n,ns1)
gi(n)=qvan(n,ns,ns1)
enddo
end if
do n=ikl+1,grid%mesh
qvan(n,ns,ns1)=0.0_dp
enddo
@ -353,7 +389,7 @@ subroutine gener_pseudo
!
if (lls(ns) == lls(ns1).and.abs(jjs(ns)-jjs(ns1)).lt.1.e-8_dp) then
nst=(lls(ns)+1)*2
qq(ns,ns1)=int_0_inf_dr(gi,grid,ikk(ns),nst)
qq(ns,ns1)=int_0_inf_dr(gi,grid,ikl,nst)
endif
!
! set the bmat with the eigenvalue part
@ -385,14 +421,6 @@ subroutine gener_pseudo
ddd(:,:,is)=bmat(:,:)
enddo
!
! Pseudize the Q functions if required. This might be needed for
! pseudo-potentials with semicore states. In this case the cut-off radius
! for the norm conserving wavefunctions is quite small and without
! the Q pseudization the augmentation charges are very hard making the
! ASR in phonon calculation very difficult to converge.
!
IF (which_augfun=='PSQ'.and..not.lpaw) CALL pseudo_q(qvan,qvanl)
!
! generate a PAW dataset if required
!
if (lpaw) then
@ -409,13 +437,21 @@ subroutine gener_pseudo
ikl=max(ikk(ns),ikk(ns1))
nst=2*(lls(ns)+1)
do n=1,ikl
gi(n,1)=psipaw(n,ns)*(enls(ns1)-vpotpaw(n))*psipaw(n,ns1)
gi(n)=psipaw(n,ns)*(enls(ns1)-vpotpaw(n))*psipaw(n,ns1)
end do
aekin(ns,ns1)=int_0_inf_dr(gi(1:grid%mesh,1),grid,ikl,nst)
aekin(ns,ns1)=int_0_inf_dr(gi(1:grid%mesh),grid,ikl,nst)
IF (rel==2) THEN
do n=1,irc
gi(n)=psipaw_rel(n,ns)*(enls(ns1)-vpotpaw(n))*&
psipaw_rel(n,ns1)
enddo
aekin(ns,ns1)=aekin(ns,ns1)+&
int_0_inf_dr(gi(1:grid%mesh),grid,irc,nst)
ENDIF
do n=1,ikl
gi(n,1)=phis(n,ns)*( (enls(ns1)-vpsloc(n))*phis(n,ns1) - chis(n,ns1) )
gi(n)=phis(n,ns)*( (enls(ns1)-vpsloc(n))*phis(n,ns1) - chis(n,ns1) )
end do
pskin(ns,ns1)=int_0_inf_dr(gi(1:grid%mesh,1),grid,ikl,nst)
pskin(ns,ns1)=int_0_inf_dr(gi(1:grid%mesh),grid,ikl,nst)
else
aekin(ns,ns1)=0._dp
pskin(ns,ns1)=0._dp
@ -428,14 +464,23 @@ subroutine gener_pseudo
! create the 'pawsetup' object containing the atomic setup for PAW
call us2paw ( pawsetup, &
zval, grid, rmatch_augfun, ikk, &
nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, phis, betas, &
qvan, kindiff, nlcc, aeccharge, psccharge, vpotpaw, vpsloc, &
which_augfun)
nbeta, lls, jjs, ocs, enls, els, rcutus, psipaw, psipaw_rel, &
phis, betas, qvan, kindiff, nlcc, aeccharge, psccharge, vpotpaw, &
vpsloc, which_augfun, rel)
!
! reread augmentation functions and descreened potentials from PAW
call paw2us ( pawsetup, zval, grid, nbeta, lls, jjs, ikk, betas, &
qq, qvan, vpsloc, bmat, rhos, els, rcutus, pseudotype )
qq, qvan, vpsloc, bmat, rhos, els, rcutus, pseudotype, &
psipaw_rel )
else
!
! Pseudize the Q functions if required. This might be needed for
! pseudo-potentials with semicore states. In this case the cut-off radius
! for the norm conserving wavefunctions is quite small and without
! the Q pseudization the augmentation charges are very hard making the
! ASR in phonon calculation very difficult to converge.
!
IF (which_augfun=='PSQ') CALL pseudo_q(qvan,qvanl)
!
! unscreen the local potential and the D coefficients
!

4
atomic/import_upf.f90
View File

@ -205,14 +205,18 @@ INTEGER :: mesh, nbeta,ih,jh,ijh
pawset_%enl(:) = 0.0_DP
if (upf_%has_so) then
pawset_%jj(1:nbeta) = upf_%jjj(1:nbeta)
pawset_%rel=2
else
pawset_%jj(:) = 0.0_DP
pawset_%rel=1
endif
pawset_%l(1:nbeta) = upf_%lll(1:nbeta)
pawset_%ikk(1:nbeta) = upf_%kbeta(1:nbeta)
pawset_%oc(1:nbeta) = upf_%paw%oc(1:nbeta)
pawset_%aewfc(1:mesh,1:nbeta) = upf_%aewfc(1:mesh,1:nbeta)
pawset_%pswfc(1:mesh,1:nbeta) = upf_%pswfc(1:mesh,1:nbeta)
IF (upf_%has_so) &
pawset_%aewfc_rel(1:mesh,1:nbeta) = upf_%paw%aewfc_rel(1:mesh,1:nbeta)
pawset_%proj(1:mesh,1:nbeta) = upf_%beta(1:mesh,1:nbeta)
DO ih = 1,nbeta

3
atomic/ld1inc.f90
View File

@ -74,6 +74,8 @@ module ld1inc
lls(nwfsx), & ! the angular momentum of pseudopotential
isws(nwfsx),& ! the spin of each pseudo-wavefunctions (not used)
ikk(nwfsx), & ! the maximum ik of each beta functions
ik(nwfsx), & ! the ik that correspond to rcut
ikus(nwfsx), & ! the ik that corresponds to rcutus
nwfs, & ! the number of pseudo wavefunctions
nbeta, & ! the number of projectors
nsloc, & ! the wavefunction which correspond to the loc pot
@ -266,6 +268,7 @@ module ld1inc
real(DP) :: &
rmatch_augfun, & ! define the matching radius for paw aug.fun.
psipaw(ndmx,nwfsx),& ! the all-electron wavefunctions for any beta
psipaw_rel(ndmx,nwfsx),& ! the all-electron wfc small component
aeccharge(ndmx), & ! true, not smoothened, AE core charge for PAW
psccharge(ndmx), & ! smoothened core charge for PAW
rCutNC2paw(nwfsx), & ! a cut-off radius for NC wavefunctions to be used

5
atomic/paw_type.f90
View File

@ -25,6 +25,7 @@ MODULE paw_type
LOGICAL :: nlcc ! nonlinear core correction
INTEGER :: nwfc ! number of wavefunctions/projectors
INTEGER :: lmax ! maximum angular momentum of projectors
INTEGER :: rel ! the relativistic level
INTEGER, POINTER :: l(:) !l(nwfsx) ! angular momentum of projectors
INTEGER, POINTER :: ikk(:) !ikk(nwfsx) ! cutoff radius for the projectors
INTEGER :: irc ! r(irc) = radius of the augmentation sphere
@ -35,6 +36,7 @@ MODULE paw_type
jj (:), & ! the total angular momentum
rcutus (:), & ! the cutoff
aewfc(:,:), & !(ndmx,nwfsx) all-electron wavefunctions
aewfc_rel(:,:), & !(ndmx,nwfsx) all-electron wavefunctions
pswfc(:,:), & !(ndmx,nwfsx) pseudo wavefunctions
proj(:,:), & !(ndmx,nwfsx) projectors
augfun(:,:,:,:), &!(ndmx,nwfsx,nwfsx,0:2*lmaxx+1),
@ -63,7 +65,7 @@ MODULE paw_type
TYPE( paw_t ), INTENT(INOUT) :: paw
CALL nullify_radial_grid( paw%grid ) ! nullify grid object, here grid is not a POINTER!
NULLIFY( paw%l, paw%ikk )
NULLIFY( paw%oc, paw%enl, paw%aewfc, paw%pswfc, paw%proj )
NULLIFY( paw%oc, paw%enl, paw%aewfc, paw%aewfc_rel, paw%pswfc, paw%proj)
NULLIFY( paw%augfun, paw%augmom, paw%aeccharge, paw%psccharge, paw%pscharge )
NULLIFY( paw%aeloc, paw%psloc, paw%dion )
NULLIFY( paw%kdiff )
@ -82,6 +84,7 @@ MODULE paw_type
ALLOCATE ( paw%els(size_nwfc) )
ALLOCATE ( paw%enl(size_nwfc) )
ALLOCATE ( paw%aewfc(size_mesh,size_nwfc) )
ALLOCATE ( paw%aewfc_rel(size_mesh,size_nwfc) )
ALLOCATE ( paw%pswfc(size_mesh,size_nwfc) )
ALLOCATE ( paw%proj (size_mesh,size_nwfc) )
ALLOCATE ( paw%augfun(size_mesh,size_nwfc,size_nwfc,0:2*size_lmax) )

5
atomic/pseudovloc.f90
View File

@ -19,7 +19,8 @@ subroutine pseudovloc
use io_global, only : stdout
use ld1inc, only : lloc, rcloc, grid, vpot, vpsloc, rel, nsloc, &
phis, els, chis, psipsus, &
jjs, nstoae, enls, new, psi, enl, rcut, psipaw
jjs, nstoae, enls, new, psi, enl, rcut, psipaw, &
psipaw_rel
implicit none
integer :: &
@ -118,7 +119,7 @@ subroutine pseudovloc
!
ns=indns(indi)
if (new(ns)) then
call set_psi_in(ik,lloc,jjs(ns),enls(ns),psi_in)
call set_psi_in(ik,lloc,jjs(ns),enls(ns),psi_in,psipaw_rel)
else
psi_in(:)=psi(:,1,nwf0)
endif

10
atomic/set_psi_in.f90
View File

@ -1,4 +1,4 @@
SUBROUTINE set_psi_in(ik,l,j,e,psi_out)
SUBROUTINE set_psi_in(ik,l,j,e,psi_out,psi_out_rel)
!
! This subroutine calculates the all electron wavefunction psi at the
! input energy e
@ -10,15 +10,18 @@ IMPLICIT NONE
INTEGER :: l, ik ! input: angular momentum and index of the cut-off radius
REAL(DP) :: e, j ! input: energy and total angular momentum
REAL(DP) :: psi_out(ndmx) ! output: the function psi.
REAL(DP) :: psi_out_rel(ndmx) ! output: the function psi (small component).
REAL(DP) :: psi_dir(ndmx,2) ! auxiliary function.
REAL(DP) :: ze2, jnor
integer :: n, nstop
psi_out_rel=0.0_DP
IF (rel == 1) THEN
CALL lschps(3,zed,grid,grid%mesh,grid%mesh,1,l,e,psi_out,vpot,nstop)
ELSEIF (rel == 2) THEN
CALL dir_outward(ndmx,grid%mesh,l,j,e,grid%dx,psi_dir,grid%r,grid%rab,vpot)
psi_out(:)=psi_dir(:,1)
psi_out_rel(:)=psi_dir(:,2)
ELSE
ze2=-zed*2.0_dp
CALL intref(l,e,grid%mesh,grid,vpot,ze2,psi_out)
@ -34,6 +37,11 @@ jnor = sqrt(jnor)
DO n=1,grid%mesh
psi_out(n)=psi_out(n)*0.5_dp/jnor
ENDDO
IF (rel==2) THEN
DO n=1,grid%mesh
psi_out_rel(n)=psi_out_rel(n)*0.5_dp/jnor
ENDDO
ENDIF
RETURN
END SUBROUTINE set_psi_in

5
atomic/start_potps.f90
View File

@ -56,7 +56,10 @@ subroutine guess_initial_wfc()
!
! This subroutine guess some initial wavefunction for the test configuration
!
use ld1inc
use kinds, only : DP
use radial_grids, only : ndmx
use ld1inc, only: nwfts, octs, llts, jjts, nstoaets, rcutts, rcutusts, grid, &
lpaw, psi, lpaw, tm, enlts, elts, phits, pseudotype, iswitch
implicit none
integer :: &