[Wien] ELF
Peter Blaha
peter.blaha at tuwien.ac.at
Fri Oct 28 13:43:54 CEST 2022
Dear Kateryna ,
In fact, I found a big difference between create_elf and
x lapw0 (with VX_ELF); x lapw5 -exchange
I traced it back to normalization errors in tau_w and tau_tf, which
missed a factor of 2.
The attached create_rho.f fixes the problem. It should be copied
into SRC_trig; make
Then you can use create_elf again.
PS: I would always compare the ELF created with both methods as
indicated above. Depending on the numerics, one or the other method may
give smoother plots, but in any case, they should be very similar.
PPS: The agreement to QE-ELF seems reasonable (but not perfect), but
I've not converged my calculations.
Thanks for the report
Peter Blaha
>I attach a pdf showing the differences. Also attached are my wien2k
>struct file and quantum espresso input file.
>Both calculations were done without spin polarization and using PBE.
>To me, the differences are big enough to question whether it is
>meaningful to use ELF at all if it depends on all-electron vs
>pseudopotential so strongly. Unless I am missing something or doing
>something wrong.
>Thank you,
>Kateryna
--
--------------------------------------------------------------------------
Peter BLAHA, Inst.f. Materials Chemistry, TU Vienna, A-1060 Vienna
Phone: +43-1-58801-165300
Email: peter.blaha at tuwien.ac.at WIEN2k: http://www.wien2k.at
WWW: http://www.imc.tuwien.ac.at
-------------------------------------------------------------------------
-------------- next part --------------
program create_rho
!Reads the rho/rho_onedim or xsf files containing the kinetic-energy densities
!tau, tauTF and tauW created by lapw5 or 3ddens (via the script create_elf_lapw),
!and creates new rho/rho_onedim or xsf files containing the function
!alpha=(tau-tauW)/tauTF, z=tauW/tau or ELF=1/(1+alpha**2).
implicit none
integer :: i, indexunit(3), irecl, iunit, ivx, j, k, &
ndim, npx, npy, nunit, x, y, z
real*8 :: relx, rely
real*8, allocatable :: f(:,:,:), f_onedim(:,:), f_xsf(:,:), &
g(:,:), g_onedim(:), g_xsf(:), r_onedim(:)
character*11 :: status, form
character*16 :: dgrid
character*77 :: fname
call getarg(2,fname)
if (fname .eq. ' ') call getarg(1,fname)
open(1,file=fname,status='old',err=8000)
8003 read(1,*,end=8001) iunit, fname, status, form, irecl
open(iunit,file=fname,status=status,form=form,err=8002)
goto 8003
8000 write(*,*) ' ERROR IN OPENING CREATE_RHO.DEF !!!!'
stop 'CREATE_RHO.DEF'
8002 write(*,*) ' ERROR IN OPENING UNIT:', iunit
write(*,*) ' FILENAME: ', fname, ' STATUS: ', status, ' FORM:', form
stop 'OPEN FAILED'
8001 continue
close (1)
read(10,*) ivx
read(10,*) ndim
if (ivx .eq. 1) then !VX_ALPHA
nunit = 3
indexunit(1:3) = (/ 1, 2, 3 /)
elseif (ivx .eq. 2) then !VX_Z
nunit = 2
indexunit(1:2) = (/ 1, 3 /)
elseif (ivx .eq. 3) then !VX_ELF
nunit = 3
indexunit(1:3) = (/ 1, 2, 3 /)
endif
if (ndim .eq. 1) then
do k=1, nunit
read(20+indexunit(k),'(2i5,2f10.5)') npx, npy, relx, rely
if (k .eq. 1) then
allocate(f(npx,npy,3),g(npx,npy))
f = 0d0
g = 0d0
endif
read(20+indexunit(k),'(5es16.8)') ((f(i,j,indexunit(k)),j=1,npy),i=1,npx)
enddo
if (npy .eq. 1) then
allocate(f_onedim(npx,3),g_onedim(npx),r_onedim(npx))
f_onedim = 0d0
g_onedim = 0d0
r_onedim = 0d0
do k=1, nunit
do i=1, npx
read(40+indexunit(k),'(2f20.9)') r_onedim(i), f_onedim(i,indexunit(k))
enddo
enddo
endif
if (ivx .eq. 1) then !VX_ALPHA
do i=1, npx
do j=1, npy
if (abs(f(i,j,2)) .gt. 1d-18) then
g(i,j) = (max(f(i,j,1),f(i,j,3)) - f(i,j,3))/f(i,j,2)
else
g(i,j) = 0d0
endif
enddo
enddo
if (npy .eq. 1) then
do i=1, npx
if (abs(f_onedim(i,2)) .gt. 1d-18) then
g_onedim(i) = (max(f_onedim(i,1),f_onedim(i,3)) - f_onedim(i,3))/f_onedim(i,2)
else
g_onedim(i) = 0d0
endif
enddo
!Fix unphysical values at the nucleus. alpha should be zero
do i=1,npx
if ((g_onedim(i) .gt. 100.d0) .and. (i .eq. 1)) then
if (g_onedim(i+1) .lt. 1.d-2) then
g_onedim(i) = 0.d0
print*, 'alpha set to zero for i=',i
endif
elseif ((g_onedim(i) .gt. 100.d0) .and. (i .eq. npx)) then
if (g_onedim(i-1) .lt. 1.d-2) then
g_onedim(i) = 0.d0
print*, 'alpha set to zero for i=',i
endif
elseif (g_onedim(i) .gt. 100.d0) then
if (g_onedim(i-1) .lt. 1.d-2) then
g_onedim(i) = 0.d0
print*, 'alpha set to zero for i=',i
endif
endif
enddo
endif
elseif (ivx .eq. 2) then !VX_Z
do i=1, npx
do j=1, npy
if (abs(max(f(i,j,1),f(i,j,3))) .gt. 1d-18) then
g(i,j) = f(i,j,3)/max(f(i,j,1),f(i,j,3))
else
g(i,j) = 0d0
endif
enddo
enddo
if (npy .eq. 1) then
do i=1, npx
if (abs(max(f_onedim(i,1),f_onedim(i,3))) .gt. 1d-18) then
g_onedim(i) = f_onedim(i,3)/max(f_onedim(i,1),f_onedim(i,3))
else
g_onedim(i) = 0d0
endif
enddo
endif
elseif (ivx .eq. 3) then !VX_ELF
do i=1, npx
do j=1, npy
if (abs(f(i,j,2)) .gt. 1d-18) then
g(i,j) = 1d0/(1d0 + ((max(f(i,j,1),f(i,j,3)*2.d0)-f(i,j,3)*2.d0)/f(i,j,2)/2.d0)**2)
else
g(i,j) = 0d0
endif
enddo
enddo
if (npy .eq. 1) then
do i=1, npx
if (abs(f_onedim(i,2)) .gt. 1d-18) then
g_onedim(i) = 1d0/(1d0 + ((max(f_onedim(i,1),f_onedim(i,3)*2.d0)-f_onedim(i,3)*2.d0)/f_onedim(i,2)/2.d0)**2)
else
g_onedim(i) = 0d0
endif
enddo
endif
endif
write(20,'(2i5,2f10.5)') npx, npy, relx, rely
write(20,'(5es16.8)') ((g(i,j),j=1,npy),i=1,npx)
if (npy .eq. 1) then
do i=1, npx
write(40,'(2f20.9)') r_onedim(i), g_onedim(i)
enddo
endif
deallocate(f,g)
if (npy .eq. 1) deallocate(f_onedim,g_onedim,r_onedim)
elseif (ndim .eq. 2) then
do
read(60+indexunit(1),'(a16)') dgrid
if (dgrid == 'DATAGRID_3D_GRID') then
read(60+indexunit(1),*) y, x, z !That is how they are given in the xsf file
exit
endif
enddo
rewind(60+indexunit(1))
allocate(f_xsf(x*y*z,3))
f_xsf = 0.0d0
do k=1, nunit
call read_xsf(f_xsf(:,indexunit(k)),x,y,z,60+indexunit(k))
enddo
allocate(g_xsf(x*y*z))
g_xsf = 0.0d0
if (ivx .eq. 1) then !VX_ALPHA
do i=1, x*y*z
if (abs(f_xsf(i,2)) .gt. 1d-18) then
g_xsf(i) = (max(f_xsf(i,1),f_xsf(i,3)) - f_xsf(i,3))/f_xsf(i,2)
else
g_xsf(i) = 0d0
endif
enddo
elseif (ivx .eq. 2) then !VX_Z
do i=1, x*y*z
if (abs(max(f_xsf(i,1),f_xsf(i,3))) .gt. 1d-18) then
g_xsf(i) = f_xsf(i,3)/max(f_xsf(i,1),f_xsf(i,3))
else
g_xsf(i) = 0d0
endif
enddo
elseif (ivx .eq. 3) then !VX_ELF
do i=1, x*y*z
if (abs(f_xsf(i,2)) .gt. 1d-18) then
g_xsf(i) = 1d0/(1d0 + ((max(f_xsf(i,1),f_xsf(i,3)*2.d0)-f_xsf(i,3)*2.d0)/f_xsf(i,2)/2.d0)**2)
else
g_xsf(i) = 0d0
endif
enddo
endif
deallocate(f_xsf)
call write_xsf(g_xsf,x,y,z,60,60+indexunit(1))
deallocate(g_xsf)
endif
end program create_rho
!-----------------------------------------------------------------
subroutine read_xsf(xsf,x,y,z,iunit)
!Reads the xsf file generated by 3ddens.
implicit none
integer :: i, j, l, x, y, z, lines, modxy, xylines, iunit
real*8 :: xsf(x*y*z)
real*8, allocatable :: xsf_tmp(:)
character*16 :: dgrid
do
read(iunit,'(a16)') dgrid
if (dgrid == 'DATAGRID_3D_GRID') then
read(iunit,*) y, x, z !That is how they are given in the xsf file
read(iunit,'(/,/,/)')
exit
endif
enddo
allocate(xsf_tmp(5))
l = 1
xylines = int(x*y/5)
modxy = mod(x*y,5)
lines = xylines
if (modxy /= 0) lines = xylines + 1
do i = 1, z
do j = 1, xylines
read(iunit,'(5(x,es15.8))') xsf_tmp(1:5)
xsf(l:l+4) = xsf_tmp(1:5)
l = l + 5
enddo
if (modxy /= 0) then
read(iunit,'(5(x,es15.8))') xsf_tmp(1:modxy)
xsf(l:l+modxy-1) = xsf_tmp(1:modxy)
l = l + modxy
endif
enddo
deallocate(xsf_tmp)
end subroutine read_xsf
!-----------------------------------------------------------------
subroutine write_xsf(xsf,x,y,z,iunit,iunit2)
!Writes the xsf file
implicit none
integer :: i, j, l, x, y, z, lines, modxy, xylines, iunit, iunit2
real*8 :: xsf(x*y*z)
character*80 :: dummy
rewind(iunit2)
do
read(iunit2,'(a80)') dummy
write(iunit,'(a80)') dummy
if (dummy(1:16) == 'DATAGRID_3D_GRID') then
do i=1, 5
read(iunit2,'(a80)') dummy
write(iunit,'(a80)') dummy
enddo
exit
endif
enddo
l = 1
xylines = int(x*y/5)
modxy = mod(x*y,5)
lines = xylines
if (modxy /= 0) lines = xylines + 1
do i = 1, z
do j = 1, xylines
write(iunit,'(5(x,es15.8))') xsf(l:l+4)
l = l + 5
enddo
if (modxy /= 0) then
write(iunit,'(5(x,es15.8))') xsf(l:l+modxy-1)
l = l + modxy
endif
enddo
write(iunit,'(a)') "END_DATAGRID_3D"
write(iunit,'(a)') "END_BLOCK_DATAGRID3D"
end subroutine write_xsf
More information about the Wien
mailing list