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

Cleanup of the mess with spin-orbit and semilocal PP 

git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@1637 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
giannozz 2005-02-11 11:18:07 +00:00
parent 5d6d65b2d0
commit 552a409910
12 changed files with 141 additions and 227 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
atomic/descreening.f90
View File

@ -12,11 +12,12 @@ subroutine descreening
implicit none
integer :: &
ns, & ! counter on pseudo functions
ns1, & ! counter on pseudo functions
integer :: &
ns, & ! counter on pseudo functions
ns1, & ! counter on pseudo functions
ib,jb, & ! counter on beta functions
lam ! the angular momentum
lam, & ! the angular momentum
ind
real(kind=dp) :: &
vaux(ndm,2),& ! work space
@ -70,9 +71,16 @@ subroutine descreening
do ns=1,nwfts
if (octs(ns).gt.0.0_dp) then
is=iswts(ns)
if (pseudotype.eq.1) then
if (pseudotype ==1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
do n=1,mesh
vaux(n,1)=vpsloc(n)+vnl(n,llts(ns))
vaux(n,1)=vpsloc(n)+vnl(n,llts(ns),ind)
enddo
else
do n=1,mesh

19
atomic/elsdps.f90
View File

@ -34,7 +34,7 @@ implicit none
egcc(:) ! the gga core energy
integer :: &
n,i,ns,is,nst,lam,n1,n2,ikl,ierr
n,i,ns,is,nst,lam,n1,n2,ikl,ierr, ind
logical :: &
gga ! if true it is a gga calculation
@ -120,8 +120,15 @@ implicit none
endif
!
epseu=0.0_DP
if (pseudotype.eq.1) then
if (pseudotype == 1) then
do ns=1,nwfts
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
f1=0.0_DP
lam=llts(ns)
if (octs(ns).gt.0.0_DP) then
@ -130,17 +137,17 @@ implicit none
enddo
endif
do n=1,mesh
f1(n,1) = f1(n,1) * vnl(n,lam)
f1(n,1) = f1(n,1) * vnl(n,lam,ind)
end do
if (ikk(ns) > 0) &
epseu = epseu + int_0_inf_dr(f1,r,r2,dx,ikk(ns),2*(lam+1))
enddo
elseif ((pseudotype.eq.2).or.(pseudotype.eq.3)) then
else
do ns=1,nwfts
if (octs(ns).gt.0.0_DP) then
do n1=1,nbeta
if (llts(ns).eq.lls(n1).and. &
abs(jjts(ns)-jjs(n1)).lt.1.e-7_DP) then
if ( llts(ns).eq.lls(n1).and. &
abs(jjts(ns)-jjs(n1)).lt.1.e-7_DP) then
nst=(llts(ns)+1)*2
ikl=ikk(n1)
do n=1,ikl

27
atomic/gener_pseudo.f90
View File

@ -54,9 +54,9 @@ subroutine gener_pseudo
integer :: &
m, n, l, n1, n2, nwf0, nst, ikl, imax, iwork(nwfsx), &
is, nbf, nc, ios
is, nbf, nc, ios, ind
character(len=5) :: ind
character(len=5) :: indqvan
logical :: &
lbes4 ! use 4 Bessel functions expansion
@ -209,8 +209,15 @@ subroutine gener_pseudo
vnl=0.0_dp
do ns=1,nbeta
lam=lls(ns)
if ( rel < 2 .or. lls(ns) == 0 .or. &
abs(jjs(ns)-lls(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. lls(ns) > 0 .and. &
abs(jjs(ns)-lls(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
do n=1,ikk(ns)
vnl(n,lam) = chis(n,ns)/phis(n,ns)
vnl(n,lam,ind) = chis(n,ns)/phis(n,ns)
enddo
enddo
!
@ -378,15 +385,15 @@ if (file_wavefunctionsps .ne. ' ') then
endif
if (file_qvan .ne. ' ') then
do ns1=1,nbeta
ind=' '
if (ns1.lt.10) then
write(ind,'(".",i1)') ns1
elseif (ns1.lt.100) then
write(ind,'(".",i2)') ns1
indqvan=' '
if (ns1 < 10) then
write(indqvan,'(".",i1)') ns1
elseif (ns1 < 100) then
write(indqvan,'(".",i2)') ns1
else
write(ind,'(".",i3)') ns1
write(indqvan,'(".",i3)') ns1
endif
open(unit=19,file=TRIM(file_qvan)//TRIM(ind), status='unknown', &
open(unit=19,file=TRIM(file_qvan)//TRIM(indqvan), status='unknown', &
iostat=ios, err=700)
700 call errore('gener_pseudo','opening file '//file_qvan,abs(ios))
do n=1,mesh

2
atomic/ld1.f90
View File

@ -15,7 +15,7 @@ program ld1
character :: &
day*9, hour*9
character(len=9), parameter:: version='08-Feb-05'
character(len=9), parameter:: version='11-Feb-05'
!
! write initialization information
!

18
atomic/ld1_readin.f90
View File

@ -399,9 +399,6 @@ subroutine ld1_readin
!
! try old Norm-Conserving format
!
call read_pseudo &
(file_pseudo,zed,xmin,rmax,dx,mesh,ndm,r,r2,sqr, &
dft,lmax,lloc,zval,nlcc,rhoc,vnl,vnlo,vpsloc,rel)
pseudotype = 1
!
else
@ -413,15 +410,24 @@ subroutine ld1_readin
else
!
! Old Norm-Conserving format
!
pseudotype = 1
!
endif
!
if (pseudotype == 1) then
!
call read_pseudo &
(file_pseudo,zed,xmin,rmax,dx,mesh,ndm,r,r2,sqr, &
dft,lmax,lloc,zval,nlcc,rhoc,vnl,vnlo,vpsloc,rel)
pseudotype = 1
dft,lmax,lloc,zval,nlcc,rhoc,vnl,vpsloc,rel)
!
do ns=1,lmax+1
ikk(ns)=mesh
enddo
endif
!
end if
!
if (lpaw) then
if (pseudotype /= 3) call errore('ld1_readin', &
'please start from a US for generating a PAW dataset' ,pseudotype)

16
atomic/ld1_setup.f90
View File

@ -17,11 +17,11 @@ subroutine ld1_setup
! transform dft in a series of codes for the exchange and
! correlation routine
!
if (iswitch.ne.2.or.pseudotype.eq.1) call which_dft(dft)
if (iswitch /= 2 .or. pseudotype == 1) call which_dft(dft)
!
! split the spin-orbit states
!
if (rel.eq.2) then
if (rel == 2) then
call occ_spinorb(nwf,nwfx,el,nn,ll,jj,oc,isw)
do n=1,nwf
oc_old(n)=oc(n)
@ -31,7 +31,7 @@ subroutine ld1_setup
! make the correspondence all-electron pseudopotential
! and split the spin-orbit states
!
if (iswitch.ge.2) then
if (iswitch >= 2) then
do nc=1, nconf
if (rel.eq.2) &
call occ_spinorb(nwftsc(nc),nwfsx,eltsc(1,nc), &
@ -63,7 +63,7 @@ subroutine ld1_setup
!
! divide the core and valence states
!
if (iswitch.eq.3) then
if (iswitch == 3) then
isws=1
nwf0=nwfs
if (rel==2) &
@ -95,7 +95,7 @@ subroutine ld1_setup
new(n)=.false.
endif
enddo
if (lloc.gt.-1) then
if (lloc > -1) then
nsloc=nwfs
nbeta=nwfs-1
if (rel==2.and.lloc.ne.0) then
@ -120,13 +120,7 @@ subroutine ld1_setup
call errore('ld1_setup','strange occupations',1)
endif
enddo
if (pseudotype.eq.1) then
do n=1,nwfs
ikk(n)=mesh
enddo
endif
endif
!
! zero the external potential and total energies
!

93
atomic/ld1_writeout.f90
View File

@ -70,16 +70,6 @@ subroutine ld1_writeout
end if
!
else
!
if (pseudotype == 1) then
!
! prepare for writing UPF file
!
if (rel == 2 ) then
call copy_ncpp_so ()
end if
!
endif
!
call write_upf(iunps)
!
@ -170,86 +160,3 @@ subroutine write_rrkj (iunps)
!
end subroutine write_rrkj
!---------------------------------------------------------------------
subroutine copy_ncpp_so ()
!---------------------------------------------------------------------
!
use ld1inc
implicit none
!
integer :: n, & ! counter on wavefunctions
nch, & ! counter on chi functions
l, & ! counter on angular momentum
ir ! counter on mesh points
real(kind=dp), external :: int_0_inf_dr
real(kind=dp) :: aux(ndm)
!
if ( lloc == 0) then
do ir=1,mesh
vpsloc(ir)=vnlo(ir,lloc,1)
enddo
else if ( lloc > 0 ) then
do ir=1,mesh
vpsloc(ir) = ((lloc+1.0_dp)*vnlo(ir,lloc,2) + &
lloc*vnlo(ir,lloc,1)) / (2.0_dp*lloc + 1.0_dp)
enddo
endif
nbeta=0
do l=0,lmax
if (l /= lloc) then
nbeta=nbeta+1
nch=0
do n=1,nwfts
if (llts(n) == l .and. abs(jjts(n)-l+0.5_dp) < 1e-3_dp) nch=n
enddo
if (l==0) nch=1
if (nch == 0) call errore('copy_ncpp_so','jj not found',1)
do ir=1,mesh
betas(ir,nbeta) = (vnlo(ir,l,1)-vpsloc(ir)) * phis(ir,nch)
enddo
lls(nbeta)=llts(nch)
jjs(nbeta)=jjts(nch)
ikk(nbeta)=mesh
do ir=mesh-1,1,-1
if (abs(betas(ir,nbeta)) < 1.e-11_dp)then
ikk(nbeta)=ir
else
goto 203
endif
enddo
203 continue
do ir = 1, mesh
aux (ir) = phis(ir, nch) * betas (ir, nbeta)
enddo
bmat(nbeta,nbeta)=1.0_dp/int_0_inf_dr(aux,r,r2,dx,mesh,2*(l+1))
if (l /= 0) then
nbeta=nbeta+1
nch=0
do n=1,nwfts
if (llts(n) == l.and.abs(jjts(n)-l-0.5_dp) < 1e-3_dp) nch=n
enddo
if (nch == 0) call errore('convert','jj not found',1)
do ir=1,mesh
betas(ir,nbeta) = (vnlo(ir,l,2)-vpsloc(ir)) * phis(ir,nch)
enddo
lls(nbeta)=llts(nch)
jjs(nbeta)=jjts(nch)
ikk(nbeta)=mesh
do ir=mesh-1,1,-1
if (abs(betas(ir,nbeta)) < 1.e-11_dp) then
ikk(nbeta)=ir
else
goto 204
endif
enddo
204 continue
do ir = 1, mesh
aux(ir) = phis(ir,nch)*betas(ir,nbeta)
enddo
bmat(nbeta,nbeta)=1.0_dp/int_0_inf_dr(aux,r,r2,dx,mesh,2*(l+1))
endif
endif
end do
return
end subroutine copy_ncpp_so

3
atomic/ld1inc.f90
View File

@ -87,8 +87,7 @@ module ld1inc
psipsus(ndm,nwfx),& ! the all-electron wavefunctions for us pseudo
rhos(ndm,2), & ! the pseudo density
zetas(ndm), & ! the pseudo magnetization
vnl(ndm,0:3), & ! the nonlocal pseudopotential
vnlo(ndm,0:3,2), & ! the nonlocal pseudopotential (spin-orbit)
vnl(ndm,0:3,2), & ! the pseudopotential in semilocal form
betas(ndm,nwfsx), & ! the projector functions
chis(ndm,nwfsx), & ! auxiliary functions
rho0, & ! value of the charge at the origin

39
atomic/lderivps.f90
View File

@ -29,7 +29,7 @@ subroutine lderivps
real(kind=dp),allocatable :: &
dlchis(:,:), & ! the logarithmic derivatives
vaux(:), & !
vaux(:), & ! auxiliary: the potential
aux(:), & ! the square of the wavefunction
al(:) ! the known part of the differential equation
@ -80,27 +80,23 @@ subroutine lderivps
nst=(lam+1)**2
nbf=nbeta
if (pseudotype == 1) then
if (rel == 2) then
if (abs(jam-lam+0.5_dp) < 1.e-2_dp .or. lam == 0 ) then
ind=1
else
ind=2
endif
do n=1,mesh
vpstot(n,is)=vpstot(n,is)+vnlo(n,lam,ind)
vaux(n)=vnlo(n,lam,ind)
enddo
else
do n=1,mesh
vpstot(n,is)=vpstot(n,is)+vnl(n,lam)
vaux(n)=vnl(n,lam)
enddo
if (rel < 2 .or. lam == 0 .or. abs(jam-lam+0.5_dp) < 0.001_dp) then
ind=1
else if (rel==2 .and. lam>0 .and. abs(jam-lam-0.5_dp)<0.001_dp) then
ind=2
endif
nbf=0.0_dp
do n=1,mesh
vaux(n) = vpstot(n,is) + vnl(n,lam,ind)
enddo
nbf=0.0
else
do n=1,mesh
vaux(n) = vpstot(n,is)
enddo
endif
do n=1,4
al(n)=vpstot(n,is)-ze2/r(n)
al(n)=vaux(n)-ze2/r(n)
enddo
call series(al,r,r2,b)
@ -123,7 +119,7 @@ subroutine lderivps
aux(2)=rr2/sqr(2)
do n=1,mesh
al(n)=( (vpstot(n,is)-e)*r2(n) + lamsq )*ddx12
al(n)=( (vaux(n)-e)*r2(n) + lamsq )*ddx12
al(n)=1.0_dp-al(n)
enddo
@ -139,11 +135,6 @@ subroutine lderivps
dlchis(ie,nc)=compute_log(aux(ikrld-3),r(ikrld),dx)
enddo
if (pseudotype == 1) then
do n=1,mesh
vpstot(n,is)=vpstot(n,is)-vaux(n)
enddo
endif
enddo
if (nspin == 2 .and. is == 1) then

2
atomic/read_newpseudo.f90
View File

@ -30,7 +30,7 @@ use funct
read( iunps, '(a75)', err=100, iostat=ios ) title
read( iunps, '(i5)',err=100, iostat=ios ) pseudotype
if (pseudotype < 1 .or. pseudotype > 3) &
if (pseudotype /= 2 .and. pseudotype /= 3) &
call errore('read_newpseudo','pseudotype is wrong',1)
read( iunps, '(2l5)',err=100, iostat=ios ) reldum, nlcc

79
atomic/read_pseudo.f90
View File

@ -1,7 +1,7 @@
!
!-----------------------------------------------------------------------
subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
r,r2,sqr,dft,lmax,lloc,zval,nlcc,rhoc,vnl,vnlo,vpsloc,rel)
r,r2,sqr,dft,lmax,lloc,zval,nlcc,rhoc,vnl,vpsloc,rel)
!-----------------------------------------------------------------------
!
use kinds, only : DP
@ -20,8 +20,7 @@ subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
rmax, & ! output: the maximum mesh value
r(ndm),r2(ndm), & ! output: the mesh
sqr(ndm), & ! output: the square root of the mesh
vnl(ndm,0:3), & ! output: the potential in numerical form
vnlo(ndm,0:3,2),& ! output: the spin-orbit potential
vnl(ndm,0:3,2), & ! output: the potential in numerical form
vpsloc(ndm), & ! output: the local pseudopotential
rhoc(ndm) ! output: the core charge
@ -103,7 +102,7 @@ subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
xmax=(mesh-1)*dx+xmin
rmax=exp(xmax)/zed
!
! and generate the mesh: this overwrite the mesh defined in the
! and generate the mesh: this overwrites the mesh defined in the
! input parameters
!
call do_mesh(rmax,zed,xmin,dx,0,ndm,mesh1,r,r2,sqr)
@ -114,6 +113,10 @@ subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
! core charge
!
if (.not.numeric) then
!
! obsolescent format with analytical coefficients
!
read( iunps, '(a)', err=300, iostat=ios ) cdum
if (nlcc) then
do ir=1, mesh
rhoc(ir)=(a_core+b_core*r2(ir))*exp(-alfa_core*r2(ir)) &
@ -130,13 +133,17 @@ subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
( alc(k,l) + r2(ir)*alc(k+3,l) )
enddo
!
! NB: the factor 2 converts from hartree to rydbergs
! NB: the factor 2 converts from hartree to rydberg
! spin-orbit not implemented for analytical PP
!
vnl(ir,l) = 2.0_dp*vnloc-2.0_dp*zval/r(ir)*(cc(1)* &
erf(r(ir)*sqrt(alpc(1))) &
+ cc(2)*erf(r(ir)*sqrt(alpc(2))))
vnl(ir,l,1) = 2.0_dp*vnloc
enddo
enddo
do ir=1,mesh
vpsloc(ir) = -2.0_dp*zval/r(ir)*(cc(1)* &
erf(r(ir)*sqrt(alpc(1))) &
+ cc(2)*erf(r(ir)*sqrt(alpc(2))))
end do
endif
if (numeric) then
@ -144,32 +151,35 @@ subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
! pseudopotentials in numerical form
!
do l = 0, lmax
if (rel.lt.2) then
read( iunps, '(a)', err=300, iostat=ios )
read( iunps, *, err=300, iostat=ios ) &
(vnl(ir,l),ir=1,mesh)
else
!
! this case is not actually implemented
read( iunps, '(a)', err=300, iostat=ios ) cdum
read( iunps, *, err=300, iostat=ios ) &
(vnl(ir,l,1),ir=1,mesh)
if (rel ==2 .and. l > 0) then
!
read( iunps, '(a)', err=300, iostat=ios ) cdum
read( iunps, *, err=300, iostat=ios ) &
(vnlo(ir,l,1),ir=1,mesh)
if (l.gt.0) then
read( iunps, '(a)', err=300, iostat=ios ) cdum
read( iunps, *, err=300, iostat=ios ) &
(vnlo(ir,l,2),ir=1,mesh)
endif
(vnl(ir,l,2),ir=1,mesh)
endif
enddo
if (lloc.eq.-1) then
!
! the local part is subtracted from the non-local part
! and stored in vpsloc - just for consistency
!
if (lloc == -1) then
read( iunps, '(a)', err=300, iostat=ios )
read( iunps, *, err=300, iostat=ios ) &
(vpsloc(ir),ir=1,mesh)
else if (rel < 2) then
vpsloc(1:mesh) = vnl(1:mesh,lloc,1)
else
vpsloc=0.0_dp
endif
vpsloc(1:mesh) = (l*vnl(1:mesh,lloc,1)+(l+1)*vnl(1:mesh,lloc,2)) / &
(2*l+1)
end if
do l=0,lmax
vnl(1:mesh,l,1) = vnl(1:mesh,l,1) - vpsloc(1:mesh)
if (rel == 2) vnl(1:mesh,l,2) = vnl(1:mesh,l,2) - vpsloc(1:mesh)
end do
!
if(nlcc) then
read( iunps, *, err=300, iostat=ios ) &
( rhoc(ir), ir=1,mesh )
@ -184,19 +194,12 @@ subroutine read_pseudo (file_pseudo,zed,xmin,rmax,dx,mesh,ndm, &
300 call errore('read_pseudo','reading pseudofile',abs(ios))
!
! all the components of the nonlocal potential beyond lmax are taken
! equal to vnl of lloc
! equal to the local part
!
if (lloc >= 0 .and. lloc <=3) then
do l=lmax+1,3
vnl(:,l)=vnl(:,lloc)
enddo
else
do l=lmax+1,3
! uncomment once the case lloc=-1 is implemented
! vnl(:,l)=vpsloc(:)
vnl(:,l)=vnl(:,lmax)
enddo
end if
do l=lmax+1,3
vnl(1:mesh,l,1)=vpsloc(1:mesh)
if (rel ==2 .and. l > 0) vnl(1:mesh,l,2)=vpsloc(1:mesh)
enddo
close(iunps)
return
end subroutine read_pseudo

50
atomic/run_pseudo.f90
View File

@ -83,22 +83,19 @@ subroutine run_pseudo
do ns=1,nwfts
if (octs(ns).gt.0.0_dp) then
is=iswts(ns)
if (pseudotype.eq.1) then
if (rel.eq.2) then
if (abs(jjts(ns)-llts(ns)+0.5_dp).lt.1.e-2_dp &
& .or. llts(ns)==0 ) then
ind=1
else
ind=2
endif
do n=1,mesh
vaux(n)=vpstot(n,is)+vnlo(n,llts(ns),ind)
enddo
if (pseudotype == 1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel==2 .and. llts(ns)>0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
else
do n=1,mesh
vaux(n)=vpstot(n,is)+vnl(n,llts(ns))
enddo
call errore('run-pseudo','spin-orbit?!?',3)
endif
do n=1,mesh
vaux(n)=vpstot(n,is)+vnl (n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n)=vpstot(n,is)
@ -168,22 +165,17 @@ subroutine run_pseudo
do ns=1,nwfts
is=iswts(ns)
if (octs(ns).gt.-1.0_dp) then
if (pseudotype.eq.1) then
if (rel.eq.2) then
if (abs(jjts(ns)-llts(ns)+0.5_dp).lt.1.e-2_dp &
& .or. llts(ns)==0 ) then
ind=1
else
ind=2
endif
do n=1,mesh
vaux(n)=vpstot(n,is)+vnlo(n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n)=vpstot(n,is)+vnl(n,llts(ns))
enddo
if (pseudotype == 1) then
if ( rel < 2 .or. llts(ns) == 0 .or. &
abs(jjts(ns)-llts(ns)+0.5_dp) < 0.001_dp) then
ind=1
else if ( rel == 2 .and. llts(ns) > 0 .and. &
abs(jjts(ns)-llts(ns)-0.5_dp) < 0.001_dp) then
ind=2
endif
do n=1,mesh
vaux(n)=vpstot(n,is)+vnl (n,llts(ns),ind)
enddo
else
do n=1,mesh
vaux(n)=vpstot(n,is)