quantum-espresso/upflib/hgh2qe.f90

855 lines
28 KiB
Fortran

!! PROGRAM WRITTEN TO CONVERT PSEUDOPOTENTIAL FROM HGH/NLCC TO UPF format
!! USED IN Quantum Espresso
!! Written by Santanu Saha, University of Basel, April 2014
!! Under the supervision S.Goedecker
!! Parts of the Program were used from the Conversion Program obtained from QE Group
!! Since the psp file of the elements do not have the rhoatom values
!! you have to provide the values as a separate file from which it is
!! interpolated.Hence,at present the value is provided from the already
!! existing files in the QE website.
!! How to use it:-
!! 1) Compile hgh2qe.f90
!! > gfortran hgh2qe.f90 lininterpol.f90
!!2) Copy the psppar file psppar.<element> to psppar
!!! > cp pspp.file psppar
!! 3) Copy the rhoatom (4 columns) from the website to rhoatom.dat
!! > cp rhoatom_file rhoatom.dat
!! 4) Execute the executible
!! > ./a.out | tee OUT_FILE
!!! Santanu and Stefan
program hgh2qe
implicit real*8 (a-h,o-z)
parameter ( lmx=5, lpmx= 4, nsmx=2 )
allocatable rr(:),rab(:),pot_loc(:),pot_nlcc(:),pot_Dij(:,:),rhoatom(:)
dimension r_l(4),gpot(4),gcore(4),psppar(6,4,2),hsep(6,4,2),lproj(12),&
nproj(4)
character(100)::author,time,functional,forline
character(3)::func,psname
character(14)::filename
character(2)::elname
character(len=2), parameter :: symbol_(94)=(/' H','He', &
'Li','Be',' B',' C',' N',' O',' F','Ne', &
'Na','Mg','Al','Si',' P',' S','Cl','Ar', &
' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni',&
'Cu','Zn','Ga','Ge','As','Se','Br','Kr', &
'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd',&
'Ag','Cd','In','Sn','Sb','Te',' I','Xe', &
'Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd',&
'Tb','Dy','Ho','Er','Tm','Yb',&
'Lu','Hf','Ta',' W','Re','Os','Ir','Pt',&
'Au','Hg','Tl','Pb','Bi','Po','At','Rn', &
'Fr','Ra','Ac','Th','Pa',' U','Np','Pu'/)
!! READ HGH pseudopotential
nproj=0
lproj=-1
call readpsp (lpx,r_l,gpot,gcore,psppar,rcore,nloc,rloc,znuc,zion,ixc, &
ipspcod,lproj,nproj,nkbeta)
if (ixc==-101130.or.ixc==11)func='pbe'
if (ixc==-20.or.ixc==1)func='lda'
!psname='hgh'
!elname=symbol_(int(znuc))
!if ( elname(1).eqv." ") then
! write(filename,'(a1,a1,a3,a1,a3,a4)')elname(2),'.',func,'-',psname,'.UPF'
!else
! write(filename,'(a2,a1,a3,a1,a3,a4)')elname,'.',func,'-',psname,'.UPF'
!end if
!open (unit=6,file=filename)
!! PSP INFO !!
write(6,'(a)')'<UPF version="2.0.1">'
write(6,'(a)')' <PP_INFO>'
write(6,'(a)')' Generated in analytical, separable form'
if (ipspcod==2.or.ipspcod==3) then
write(author,'(a34)')'Goedecker/Hartwigsen/Hutter/Teter"'
write(6,'(a)')' Author: Goedecker/Hartwigsen/Hutter/Teter'
write(time,'(a43)')'Phys.Rev.B58, 3641 (1998);B54, 1703 (1996)"'
write(6,'(a)')' Generation date: Phys.Rev.B58, 3641 (1998); &
B54, 1703 (1996)'
else if (ipspcod==10) then
write(author,'(a40)')'Krack/Goedecker/Hartwigsen/Hutter/Teter"'
write(6,'(a)')' Author: Krack/Goedecker/Hartwigsen/Hutter/Teter'
write(time,'(a35)')'Theor. Chem. Acc(2005) 114:145-152"'
write(6,'(a)')' Generation date: Theor. Chem. Acc(2005) 114:145-152'
else if (ipspcod==12) then
write(author,'(a50)')'Willand/Goedecker/Genovese/Deutsch/Sadeghi/Deb/Mayagoitia/Kvashnin"'
write(6,'(a)')' Author:Willand/Goedecker/Genovese/Deutsch/Sadeghi/Deb/ &
Mayagoitia/Kvashnin'
write(time,'(a34)')'J. Chem. Phys. 138, 104109 (2013)"'
write(6,'(a)')' Generation date: J. Chem. Phys. 138, 104109 (2013)'
end if
write(6,'(a)')' Pseudopotential type: NC'
if (int(znuc) <= 0 .or.int(znuc) > 94) then
stop "Wrong zatom value"
end if
write(6,'(a13,a2)')' Element: ',symbol_(int(znuc))
if (ixc==-20) then
write(6,'(a)')' Functional: SLA-PZ-NOGX-NOGC'
write(functional,'(a17)')'SLA-PZ-NOGX-NOGC"'
else if (ixc==-101130) then
write(functional,'(a15)')'SLA-PW-PBX-PBC"'
write(6,'(a)')' Functional: SLA-PW-PBX-PBC'
end if
write(6,*)
write(6,'(a)')' Suggested minimum cutoff for wavefunctions: 0. Ry'
write(6,'(a)')' Suggested minimum cutoff for charge density: 0. Ry'
write(6,'(a)')' The Pseudo was generated with a Non-Relativistic Calculation'
write(6,'(a)')' Local Potential: unknown format, L component and cutoff &
radius: -3 0.0000'
write(6,*)
write(6,'(a)')' Valence configuration:'
write(6,'(a)')' nl pn l occ Rcut Rcut US E pseu'
write(6,*)
write(6,'(a)')' Generation configuration: not available.'
write(6,'(a)')' Converted from CPMD format using cpmd2upf v.5.0.1 - &
PG 10Jul2012'
write(6,'(a)')' Contains atomic orbitals generated by ld1.x - use with care'
write(6,'(a)')' </PP_INFO>' !! END OF PSP INFO
write(6,'(a)')' <!-- -->'
write(6,'(a)')' <!-- END OF HUMAN READABLE SECTION -->'
write(6,'(a)')' <!-- -->'
write(6,'(a)')' <PP_HEADER generated="Generated in analytical, separable form"'
write(6,'(a21,a)')' author="',author
write(6,'(a19,a)')' date="',time
write(6,'(a)')' comment="Contains atomic orbitals generated by ld1.x - use with care"'
write(6,'(a22,a2,a1)')' element="',symbol_(int(znuc)),'"'
write(6,'(a)')' pseudo_type="NC"'
write(6,'(a)')' relativistic="no"'
write(6,'(a)')' is_ultrasoft="F"'
write(6,'(a)')' is_paw="F"'
write(6,'(a)')' is_coulomb="F"'
write(6,'(a)')' has_so="F"'
write(6,'(a)')' has_wfc="F"'
write(6,'(a)')' has_gipaw="F"'
write(6,'(a)')' paw_as_gipaw="F"'
if (ipspcod==12) then
write(6,'(a)')' core_correction="T"'
else
write(6,'(a)')' core_correction="F"'
end if
write(6,'(a25,a)')' functional="',functional
write(6,'(a24,1pe22.15e3,a)')' z_valence="',zion,'"'
write(6,'(a)')' total_psenergy="0.000000000000000E+000"'
write(6,'(a)')' wfc_cutoff="0.000000000000000E+000"'
write(6,'(a)')' rho_cutoff="0.000000000000000E+000"'
if((lpx-1).ge.0)write(6,'(a20,i1,a)')' l_max="',lpx-1,'"'
if((lpx-1).lt.0)write(6,'(a20,i2,a)')' l_max="',lpx-1,'"'
write(6,'(a)')' l_max_rho="0"'
write(6,'(a)')' l_local="-3"'
!! SETTING UP RADIAL GRID
z=znuc
fact=1.d0
xmin = -7.0d0
amesh=0.0125d0*fact !! Modify this value to increase the no of grid points
rmax =100.0d0
mesh = 1 + (log(z*rmax)-xmin)/amesh
mesh = (mesh/2)*2+1 ! mesh is odd (for historical reasons?)
allocate (rr(mesh),rab(mesh),pot_loc(mesh),pot_nlcc(mesh),pot_Dij(nkbeta,nkbeta),&
rhoatom(mesh))
do i=1, mesh
rr(i) = exp (xmin+(i-1)*amesh)/z
!!r(i)=exp[(i-1)*amesh]*r(1)
rab(i)= rr(i)*amesh
end do
if(mesh.lt.1000)write(6,'(a24,i3,a1)')' mesh_size="',mesh,'"'
if((mesh.ge.1000).and.(mesh.lt.10000))write(6,'(a24,i4,a1)')' mesh_size="',mesh,'"'
if(mesh.ge.10000)write(6,'(a24,i5,a1)')' mesh_size="',mesh,'"'
write(6,'(a)')' number_of_wfc=" 0"'
if(nkbeta.lt.10)write(6,'(a29,i1,a3)')' number_of_proj="',nkbeta,'"/>'
if(nkbeta.ge.10)write(6,'(a29,i2,a3)')' number_of_proj="',nkbeta,'"/>'
rmax=rr(mesh)
xmin=log(z*rr(1))
if(mesh.lt.1000) then
50 format (a15,1pe22.15e3,a8,i3,a8,1pe23.15e3,a8,1pe22.15e3,a1)
write(6,50)' <PP_MESH dx="',amesh,'" mesh="',mesh,'" xmin="',xmin,'" rmax="',rmax,'"'
write(6,'(a7,1pe22.15e3,a2)')'zmesh="',z,'">'
write(6,'(a28,i3,a)')' <PP_R type="real" size="',mesh,'" columns="4">'
write(6,'(4(1x,1pe22.15e3,1x))')rr
write(6,'(a)')' </PP_R>'
write(6,'(a30,i3,a)')' <PP_RAB type="real" size="',mesh,'" columns="4">'
else if((mesh.ge.1000).and.(mesh.lt.10000)) then
51 format (a15,1pe22.15e3,a8,i4,a8,1pe23.15e3,a8,1pe22.15e3,a1)
write(6,51)' <PP_MESH dx="',amesh,'" mesh="',mesh,'" xmin="',xmin,'" rmax="',rmax,'"'
write(6,'(a7,1pe22.15e3,a2)')'zmesh="',z,'">'
write(6,'(a28,i4,a)')' <PP_R type="real" size="',mesh,'" columns="4">'
write(6,'(4(1x,1pe22.15e3,1x))')rr
write(6,'(a)')' </PP_R>'
write(6,'(a30,i4,a)')' <PP_RAB type="real" size="',mesh,'" columns="4">'
else if(mesh.ge.10000) then
52 format (a15,1pe22.15e3,a8,i5,a8,1pe23.15e3,a8,1pe22.15e3,a1)
write(6,52)' <PP_MESH dx="',amesh,'"mesh="',mesh,'" xmin="',xmin,'" rmax="',rmax,'"'
write(6,'(a7,1pe22.15e3,a2)')'zmesh="',z,'">'
write(6,'(a28,i5,a)')' <PP_R type="real" size="',mesh,'" columns="4">'
write(6,'(4(1x,1pe22.15e3,1x))')rr
write(6,'(a)')' </PP_R>'
write(6,'(a30,i5,a)')' <PP_RAB type="real" size="',mesh,'" columns="4">'
end if
write(6,'(4(1x,1pe22.15e3,1x))')rab
write(6,'(a)')' </PP_RAB>'
write(6,'(a)')' </PP_MESH>'
!! END of setting up and writing Radial Grid and Radial Weights
call vloc (mesh,nloc,rloc,zion,rr,gpot,pot_loc) !! Local Potential
if (ipspcod==12)call vnlcc (mesh,znuc,zion,rcore,gcore,rr,pot_nlcc) !! NLCC
!! Non-Local Potential
call vnonloc (mesh,nkbeta,ipspcod,lpx,psppar,r_l,rr,pot_Dij,lproj,nproj)
write(6,'(a)')' <PP_PSWFC>'
write(6,'(a)')' </PP_PSWFC>'
! fourpi=16.d0*datan(1.d0)
! totrho=0.d0
! rcov=1.304d0
! do imesh=1,mesh
! r=rr(imesh)
! rhoatom(imesh)=fourpi*r**2*r**(3.d0/2.d0)*exp(-0.5d0*((r-0.52*rcov)/rcov)**2)
! totrho=totrho+rhoatom(imesh)*rab(imesh)
! end do
! anorm=zion/totrho
! rhoatom=anorm*rhoatom
! !write(21,*)'totrho,anorm,anorm1,zion',totrho,anorm,anorm1,zion
! !write(21,*)'diff',abs(totrho-totrho1)
! !write(21,'(4e25.17)')rhoatom
! !testrho=0.d0
! !do imesh=1,mesh
! ! testrho=testrho+rhoatom(imesh)*rab(imesh)
! !end do
! !write(21,*)'testrho',testrho
call interpolate(znuc,mesh,amesh,rr,rhoatom)
if(mesh.lt.1000)write(6,'(a32,i3,a)')' <PP_RHOATOM type="real" size="',mesh,'" columns="4">'
if((mesh.ge.1000).and.(mesh.lt.10000))write(6,'(a32,i4,a)')' <PP_RHOATOM type="real" size="',mesh,'" columns="4">'
if(mesh.ge.10000)write(6,'(a32,i5,a)')' <PP_RHOATOM type="real" size="',mesh,'" columns="4">'
write(6,'(4(1x,1pe22.15e3,1x))')rhoatom
write(6,'(a)')' </PP_RHOATOM>'
write(6,'(a)')'</UPF>'
close(6)
end program hgh2qe
!! LOCAL PART of PSP
subroutine vloc (mesh,nloc,rloc,zion,rr,gpot,pot_loc)
implicit real*8 (a-h,o-z)
dimension rr(mesh),gpot(4),pot_loc(mesh)
pot_loc=0.d0
do i=1,mesh
a=rr(i)/rloc
pot_loc(i)=-zion*erf(a/dsqrt(2.d0))/rr(i)
factor=0.d0
do j=1,nloc
factor=factor+gpot(j)*a**(2.d0*j-2.d0)
end do
pot_loc(i)=pot_loc(i)+factor*exp(-0.5d0*a*a)
end do
pot_loc=pot_loc*2.d0
if(mesh.lt.1000)write(6,'(a30,i3,a)')' <PP_LOCAL type="real" size="',mesh,'" columns="4">'
if((mesh.ge.1000).and.(mesh.lt.10000))write(6,'(a30,i4,a)')' <PP_LOCAL type="real" size="',mesh,'" columns="4">'
if(mesh.ge.10000)write(6,'(a30,i5,a)')' <PP_LOCAL type="real" size="',mesh,'" columns="4">'
write(6,'(4(1pe23.15e3,2x))')pot_loc
write(6,'(a)')' </PP_LOCAL>'
end subroutine vloc
!! End of setting Vloc
!! NLCC PART OF PSP
subroutine vnlcc (mesh,znuc,zion,rcore,gcore,rr,pot_nlcc)
implicit real*8 (a-h,o-z)
dimension rr(mesh),pot_nlcc(mesh),gcore(4)
pot_nlcc=0.d0
pi=4.d0*datan(1.d0)
sqrt2pi=dsqrt(2.d0*pi)
factor=gcore(1)/(sqrt2pi*rcore)**3.d0
!!factor=gcore(1)
do i=1,mesh
a=(rr(i)/rcore)**2.d0
pot_nlcc(i)=factor*exp(-a*0.5d0)
end do
!!pot_nlcc=2.d0*pot_nlcc
write(6,*)' <PP_NLCC>'
write(6,'(4(1x,1pe22.15e3,1x))')pot_nlcc
write(6,*)' </PP_NLCC>'
end subroutine vnlcc
!! End of setting Vnlcc
!! NONLOCAL PART OF PSP
subroutine vnonloc (mesh,nkbeta,ipspcod,lpx,psppar,r_l,rr,pot_Dij,lproj,nproj)
implicit real*8 (a-h,o-z)
dimension rr(mesh),psppar(6,4,2),pot_Dij(nkbeta,nkbeta),H(3,3),r_l(lpx), &
hsep(6,4,2),ofdcoef(3,4),lproj(nkbeta),proj(mesh,nkbeta),nproj(4), &
rcut(nkbeta),icut(nkbeta),pot_loc(mesh) !! 3 projectors * 4 l values
character(2):: label
hsep=0.d0
!! Setting Up the Projectors
nsproj=0;npproj=0;ndproj=0;nfproj=0;jproj=0
do iproj=1,nkbeta
if (lproj(iproj)==0) then
nsproj=nsproj+1
jproj=nsproj
else if (lproj(iproj)==1) then
npproj=npproj+1
jproj=npproj
else if (lproj(iproj)==2) then
ndproj=ndproj+1
jproj=ndproj
else if (lproj(iproj)==3) then
nfproj=nfproj+1
jproj=nfproj
end if
lq=lproj(iproj)
rnrm=1.d0/dsqrt(.5d0*gamma(lq+1.5d0+2.d0*jproj-2.d0)* &
r_l(lq+1)**(2*lq+3+4*(jproj-1)))
! gamma is standard f2008 - use instead mygamma(i) == gamma(i-1/2)
!rnrm=1.d0/dsqrt(.5d0*mygamma(lq+2*jproj)* &
! r_l(lq+1)**(2*lq+3+4*(jproj-1)))
do imesh=1,mesh
r=rr(imesh)
pr=rnrm*r**(lq+2*(jproj-1))*exp(-.5d0*(r/r_l(lq+1))**2)
proj(imesh,iproj)=2.d0*pr*r
end do
do imesh=mesh,1,-1
if (abs(proj(imesh,iproj)).gt.1d-12)exit
end do
if (imesh<2) then
write(*,*)'Error in Rcut'
else if (mod(imesh,2).ne.0) THEN
! even index
icut(iproj) = imesh
rcut(iproj)=rr(icut(iproj))
else if ((imesh.lt.mesh).and.(mod(imesh,2).eq.0)) then
! odd index
icut(iproj) = imesh+1
rcut(iproj)=rr(icut(iproj))
else
icut(iproj) = mesh
rcut(iproj)=rr(icut(iproj))
end if
end do
!! End of Setting up of Projectors
write(6,'(a)')' <PP_NONLOCAL>'
nsproj=0;npproj=0;ndproj=0;nfproj=0
do ikbeta=1,nkbeta
if(lproj(ikbeta)==0) then
nsproj=nsproj+1
write(label,'(i1,a1)')nsproj,'S'
else if(lproj(ikbeta)==1) then
npproj=npproj+1
write(label,'(i1,a1)')npproj,'P'
else if(lproj(ikbeta)==2) then
ndproj=ndproj+1
write(label,'(i1,a1)')ndproj,'D'
else if(lproj(ikbeta)==3) then
nfproj=nfproj+1
write(label,'(i1,a1)')nfproj,'F'
end if
if (icut(ikbeta).lt.1000) then
if (mesh.lt.1000) then
55 format (a13,i1,a19,i3,a21,i1,a9,a2,a20,i1,a23,i3,a1)
write(6,55)'<PP_BETA.',ikbeta,'type="real" size="',mesh,'" columns="4" index="',ikbeta,'" label="',&
label,'" angular_momentum="',lproj(ikbeta),'" cutoff_radius_index="',icut(ikbeta),'"'
else
56 format (a13,i1,a19,i4,a21,i1,a9,a2,a20,i1,a23,i3,a1)
write(6,56)'<PP_BETA.',ikbeta,'type="real" size="',mesh,'" columns="4" index="',ikbeta,'" label="',&
label,'" angular_momentum="',lproj(ikbeta),'" cutoff_radius_index="',icut(ikbeta),'"'
end if
else
if (mesh.lt.10000) then
57 format (a13,i1,a19,i4,a21,i1,a9,a2,a20,i1,a23,i4,a1)
write(6,57)'<PP_BETA.',ikbeta,'type="real" size="',mesh,'" columns="4" index="',ikbeta,'" label="',&
label,'" angular_momentum="',lproj(ikbeta),'" cutoff_radius_index="',icut(ikbeta),'"'
else
58 format (a13,i1,a19,i5,a21,i1,a9,a2,a20,i1,a23,i4,a1)
write(6,58)'<PP_BETA.',ikbeta,'type="real" size="',mesh,'" columns="4" index="',ikbeta,'" label="',&
label,'" angular_momentum="',lproj(ikbeta),'" cutoff_radius_index="',icut(ikbeta),'"'
end if
end if
write(6,'(a15,1pe22.15e3,a)')'cutoff_radius="',rcut(ikbeta),'" ultrasoft_cutoff_radius="0.000000000000000E+000">'
write(6,'(4(1x,1pe22.15e3,1x))')proj(:,ikbeta)
write(6,'(a14,i1,a)')'</PP_BETA.',ikbeta,'>'
end do
!! Setting Up the hsep elements of HGH psp
if (ipspcod == 2) then !GTH case
! offdiagonal elements are zero per definition.
! simply rearrange hij and fill zero elements
do l=1,lpx
hsep(1,l,1)=psppar(1,l,1) !h11
hsep(2,l,1)=0.0d0 !h12
hsep(3,l,1)=psppar(2,l,1) !h22
hsep(4,l,1)=0.0d0 !h13
hsep(5,l,1)=0.0d0 !h23
hsep(6,l,1)=psppar(3,l,1) !h33
! in the polarized or relativistic case,
! we assume all kij to be zero,
! i.e. same up and down projectors
end do
else if (ipspcod == 3) then !HGH diagonal case
! we need to compute the offdiagonal elements with the following coeffs
ofdcoef(1,1)=-0.5d0*sqrt(3.d0/5.d0) !h2
ofdcoef(2,1)=0.5d0*sqrt(5.d0/21.d0) !h4
ofdcoef(3,1)=-0.5d0*sqrt(100.0d0/63.d0) !h5
ofdcoef(1,2)=-0.5d0*sqrt(5.d0/7.d0) !h2
ofdcoef(2,2)=1.d0/6.d0*sqrt(35.d0/11.d0) !h4
ofdcoef(3,2)=-7.d0/3.d0*sqrt(1.d0/11.d0) !h5
ofdcoef(1,3)=-0.5d0*sqrt(7.d0/9.d0) !h2
ofdcoef(2,3)=0.5d0*sqrt(63.d0/143.d0) !h4
ofdcoef(3,3)=-9.d0*sqrt(1.d0/143.d0) !h5
ofdcoef(1,4)=0.0d0 !h2
ofdcoef(2,4)=0.0d0 !h4
ofdcoef(3,4)=0.0d0 !h5
! this could possibly be done in a simpler way ...
do l=1,lpx
hsep(1,l,1)=psppar(1,l,1)
hsep(2,l,1)=psppar(2,l,1)*ofdcoef(1,l)
hsep(3,l,1)=psppar(2,l,1)
hsep(4,l,1)=psppar(3,l,1)*ofdcoef(2,l)
hsep(5,l,1)=psppar(3,l,1)*ofdcoef(3,l)
hsep(6,l,1)=psppar(3,l,1)
end do
end if
if (ipspcod>9) then !HGH-K or HGH case,
! psppar holds hij and kij in HGHK convention
! fill hsep(up,dn) upper diagonal col by col, as needed for the fit
! for a nonrelativistic calculation, discard the kij elements
do l=1,lpx
hsep(1,l,1)=psppar(1,l,1) !h11
hsep(2,l,1)=psppar(4,l,1) !h12
hsep(3,l,1)=psppar(2,l,1) !h22
hsep(4,l,1)=psppar(5,l,1) !h13
hsep(5,l,1)=psppar(6,l,1) !h23
hsep(6,l,1)=psppar(3,l,1) !h33
! hsep(1,l,1)=psppar(1,l,1) !h11
! hsep(2,l,1)=psppar(4,l,1) !h12
! hsep(3,l,1)=psppar(2,l,1) !h22
! hsep(4,l,1)=psppar(5,l,1) !h13
! hsep(5,l,1)=psppar(6,l,1) !h23
! hsep(6,l,1)=psppar(3,l,1) !h33
end do
end if
H=0.d0
pot_Dij=0.d0
! write(2,*)"l,ix,iy,ip,jp,pot_Dij ,H"
do l=1,lpx
H(1,1)=hsep(1,l,1)
H(1,2)=hsep(2,l,1)
H(2,2)=hsep(3,l,1)
H(1,3)=hsep(4,l,1)
H(2,3)=hsep(5,l,1)
H(3,3)=hsep(6,l,1)
do ip=1,nproj(l)
do jp=ip,nproj(l)
ix=ip;iy=jp
if (l==2)ix=ix+nproj(1)
if (l==2)iy=iy+nproj(1)
if (l==3)ix=ix+nproj(1)+nproj(2)
if (l==3)iy=iy+nproj(1)+nproj(2)
if (l==4)ix=ix+nproj(1)+nproj(2)+nproj(3)
if (l==4)iy=iy+nproj(1)+nproj(2)+nproj(3)
pot_Dij(ix,iy)=H(ip,jp)
! write(2,*)l,ix,iy,ip,jp,pot_Dij(ix,iy),H(ip,jp)
end do
enddo
! write(2,*)" "
end do
pot_Dij=pot_Dij/2.d0
do i=1,nkbeta
do j=1,i-1
pot_Dij(i,j)=pot_Dij(j,i)
end do
end do
if ((nkbeta*nkbeta).lt.10) then
write(6,'(a30,i1,a14)')'<PP_DIJ type="real" size="',nkbeta*nkbeta,'" columns="4">'
else
write(6,'(a30,i2,a14)')'<PP_DIJ type="real" size="',nkbeta*nkbeta,'" columns="4">'
end if
write(6,'(4(1pe23.15e3,2x))')pot_Dij
write(6,'(a)')' </PP_DIJ>'
write(6,'(a)')' </PP_NONLOCAL>'
!
end subroutine vnonloc
!! Reading the psppar file of HGH pseudopotential
subroutine readpsp (lpx,r_l,gpot,gcore,psppar,rcore,nloc,rloc,znuc,zion,ixc, &
ipspcod,lproj,nproj,nkbeta)
implicit real*8 (a-h,o-z)
parameter ( lmx=5, lpmx= 4, nsmx=2 )
dimension r_l2(4),r_l(4),gpot(4),gcore(4),psppar(6,4,2),&
lproj(12),nproj(4)
character (200) :: string
logical exists
inquire(file='psppar',exist=exists)
if(.not.exists)then
write(2,*)'No psppar file to convert'
stop
end if
open(unit=11,file='psppar',form='formatted', &
status='unknown')
!The first line is usually for comments only.
gcore=0d0
rcore=-1d0
!Read 1st line into a string
read(11,'(a)',iostat=ierr) string
if(ierr/=0)then
write(2,*)
write(2,*)' NOTICE'
write(2,*)
write(2,*)'The first line of psppar is just comment'
! stop
ierr=0
end if
read(11,*,iostat=ierr) znuc, zion
! write(2,*)int(znuc),int(zion),"znuc,zion"
if(ierr/=0)then
write(2,*)
write(2,*)' WARNING'
write(2,*)'Could not read nuclear and valence charges'
write(2,*)'on the second line of psppar.'
stop
end if
read(11,*,iostat=ierr) ipspcod,ixc
! write(2,*)ipspcod,ixc,"ipspcod,ixc"
if(ierr/=0)then
! no need to call error handler, shared input file
write(2,*)
write(2,*)' WARNING'
write(2,*)'Could not read PSP format and iXC from'
write(2,*)'the third line of psppar.'
stop
end if
!for convenience: convert LDA and PBE ixc from
!abinit to libxc
if (ixc==1) then
! write(2,*)'LDA pade: ixc = 1 or -20 are equivalent.'
ixc=-20
end if
if (ixc==11) then
! write(2,*)'PBE: ixc = 11 or -101130 are equivalent.'
ixc=-101130
end if
psppar = 0.d0
if (ipspcod==10.or.ipspcod==12)then
! write(2,*)'HGH matrix format'
! ! HGH-K format: all projector elements given.
! dimensions explicitly defined for nonzero output.
! ! local part
gpot=0.d0
read(11,*)rloc,nloc,(gpot(j),j=1,nloc)
!write(2,*)rloc,nloc,(gpot(j),j=1,nloc)
read(11,*)lpx
!write(6,*)lpx,'lpx'
! lpx is here equivalent to nsep. Previous versions used
! it as the max l quantum number, subracting one.
! Here, 0 does not mean s only, but a purely local psppar.
! Negative is no longer used for local, but for r_l2.
if (lpx-1 .gt. lmx ) then
! write(2,*) 'array dimension problem: lpx,lpmx',lpx,lpmx
stop
end if
! ! separable part
! ! relativistic; hij are averages, kij separatons of hsep
do l=1,lpx !l channels
! add some line to read r_l2 if nprl < 0
read(11,'(a)') string
! write(2,*)string
read(string,*) r_l(l),nprl
nproj(l)=nprl
if(nprl>0)then
read(string,*) r_l(l),nprl, &
psppar(1,l,1),(psppar(j+2,l,1), &
j=2,nprl) !h_ij 1st line
do i=2,nprl
! spin up
read(11,*) psppar(i,l,1),(psppar(i+j+1,l,1), &
j=i+1,nprl)!remaining h_ij
! write(2,*) i,psppar(i,l,1),(i+j+1, psppar(i+j+1,l,1), &
! j=i+1,nprl)!remaining h_ij
end do
! do i=1,6
! write(2,*)"l=",l,"i=",i,psppar(i,l,1)
! end do
! disable r_l2, i.e set it equal to r_l
r_l2(l) = r_l(l)
else if (nprl<0) then
! if nprl is negative, read r_l2 from the 2nd line of hij
nprl=-nprl
read(string,*) r_l(l),i, &
psppar(1,l,1),(psppar(j+2,l,1), &
j=2,nprl) !h_ij 1st line
read(11,*)r_l2(l), psppar(2,l,1),(psppar(2+j+1,l,1), &
j=2+1,nprl)!2nd line
if(nprl==3) read(11,*) psppar(3,l,1)! third line
end if
! there are no kij in the s-projector
if (l==1) cycle
do i=1,nprl
read(11,*) psppar(i,l,2),(psppar(i+j+1,l,2), &
j=i+1,nprl)!all k_ij
! write(2,*) psppar(i,l,2),(psppar(i+j+1,l,2), &
! j=i+1,nprl)!all k_ij
end do
end do ! loop over l
do l=1,lpx
k=nproj(l)
do n=k,1,-1
if (psppar(n,l,1).eq.0.d0)nproj(l)=nproj(l)-1
end do
end do
if(ipspcod==12)then
! this psppar uses nlcc
read(11,*,iostat=ierr) rcore, qcore
! write(2,*) rcore, qcore,"rcore,qcore"
if(ierr/=0)then
! write(2,*)' pspcod=12 implies nlcc data on the last &
! line,'
! write(2,*)' but rcore and qcore could not be read!'
! rcore= -1d0
else
!compute gcore(1) from qcore. gcore(2:4) are
!always zero, but we keep them for future testing.
fourpi = 16.d0*datan(1.d0)
sqrt2pi = dsqrt(fourpi*0.5d0)
!gcore(1) = fourpi* qcore * (znuc-zion) / &
! (sqrt2pi*rcore)**3
!qcore=qcore*fourpi/(sqrt2pi*rcore)**3!! qcore
!!stores the value of constant and the qcore itself
gcore(1)=(znuc-zion)*qcore
end if
end if
elseif(ipspcod==3)then
! write(2,*)'HGH diagonal format'
! HGH diagonal part case
! technically, lpx is fixed at the max value of
lpx=4
gpot=0.d0
read(11,*) rloc,(gpot(j),j=1,4)
! write(2,*) rloc,(gpot(j),j=1,4)
read(11,*) r_l(1),psppar(1:3,1,1)
! write(2,*) r_l(1),psppar(1:3,1,1)
nproj(1)=3
do i=3,1,-1
if (psppar(i,1,1)==0.d0)&
nproj(1)=nproj(1)-1
end do
do l=2,4
read(11,*) r_l(l),psppar(1:3,l,1)
! write(2,*) r_l(l),psppar(1:3,l,1)
nproj(l)=3
do i=3,1,-1
if (psppar(i,l,1)==0.d0)&
nproj(l)=nproj(l)-1
end do
read(11,*) psppar(1:3,l,2)
! write(2,*) psppar(1:3,l,2)
end do
elseif(ipspcod==2)then
! write(2,*)'GTH format'
! ! GTH case
! technically, lpx is fixed at s and p
lpx=2
gpot=0.d0
read(11,*) rloc,(gpot(j),j=1,4)
! write(2,*) rloc,(gpot(j),j=1,4)
read(11,*) r_l(1),psppar(1:2,1,1)
nproj(1)=2
if (r_l(1).eq.0.d0) then
nproj(1)=0
else
if (psppar(1,1,1).eq.0.d0)nproj(1)=0
if (psppar(2,1,1).eq.0.d0)nproj(1)=1
end if
! write(2,*) r_l(1),psppar(1:2,1,1)
read(11,*) r_l(2),psppar(1 ,2,1)
nproj(2)=1
if (r_l(2).eq.0.d0) then
nproj(2)=0
else
if (psppar(2,1,1).eq.0.d0)nproj(2)=0
end if
! write(2,*) r_l(2),psppar(1 ,2,1)
! for convenience, if we have no p projector:
if(psppar(1,2,1)<1.d-5) lpx=1
else
! no need to call error handler, shared input file
write(2,*)' WARNING'
write(2,*)'pspcod (1st number of 3rd line) read from'
write(2,*)'psppar is unknown or not supported.'
write(2,*)'supported are 2,3, or 10,12 not ',ipspcod
stop
end if
!done with reading psppar
close(11)
iproj=0
do l=1,lpx
do i=1,nproj(l)
iproj=iproj+1
lproj(iproj)=l-1
end do
end do
nkbeta=sum(nproj)
end subroutine readpsp
subroutine interpolate(znuc,mesh2,amesh2,rr2,rhoatom2)
implicit real*8 (a-h,o-z)
dimension :: rr2(mesh2),rhoatom2(mesh2)
allocatable rr1(:),rhoatom1(:)
z=znuc
xmin = -7.0d0
amesh1=0.0125d0 !! Modify this value to increase the no of grid points
rmax =100.0d0
mesh1 = 1 + (log(z*rmax)-xmin)/amesh1
mesh1 = (mesh1/2)*2+1 ! mesh is odd (for historical reasons?)
allocate(rr1(mesh1),rhoatom1(mesh1))
do i=1,mesh1
rr1(i) = exp (xmin+(i-1)*amesh1)/z
end do
open(12,file='rhoatom.dat')
ir=1
do i=1,(mesh1/4)+1
if ((ir+3).le.mesh1) then
!write(*,*)"ir begin",ir
read(12,*)(rhoatom1(j),j=ir,ir+3)
!write(*,*)i,(rhoatom1(j),j=ir,ir+3)
ir=j
!write(*,*)"ir",ir
else if ((ir+2).le.mesh1) then
!write(*,*)"ir begin",ir
read(12,*)(rhoatom1(j),j=ir,ir+2)
!write(*,*)i,(rhoatom1(j),j=ir,ir+2)
ir=j
!write(*,*)"ir",ir
else if ((ir+1).le.mesh1) then
!write(*,*)"ir begin",ir
read(12,*)(rhoatom1(j),j=ir,ir+1)
!write(*,*)i,(rhoatom1(j),j=ir,ir+1)
ir=j
!write(*,*)"ir",ir
else if (ir.le.mesh1) then
!write(*,*)"ir begin",ir
read(12,*)(rhoatom1(j),j=ir,ir+0)
!write(*,*)i,(rhoatom1(j),j=ir,ir+0)
ir=j
!write(*,*)"ir",ir
end if
end do
close (12)
do i=1,mesh1
write(25,*)rr1(i),rhoatom1(i)
end do
rhoatom2=0.d0
!!do i2=1,mesh2
!! LINEAR INTERPOLATION
! do i1=1,mesh1
! if((rr2(i2).lt.rr1(i1)))rlow=rr1(i1-1)
! if((rr2(i2).lt.rr1(i1)))ilow=i1-1
! if((rr2(i2).lt.rr1(i1)))exit
! end do
! do i1=mesh1,1
! if((rr2(i2).gt.rr1(i1)))rhigh=rr1(i1+1)
! if((rr2(i2).gt.rr1(i1)))ihigh=i1+1
! if((rr2(i2).gt.rr1(i1)))exit
! end do
! dx1=rhigh-rlow
! dx2=rr2(i2)-rlow
! drhoatom=rhoatom1(ihigh)-rhoatom1(ilow)
! rhoatom2(i2)=rhoatom1(ilow)+drhoatom*dx2/dx1
! if (rhoatom2(i2).lt.0.d0)rhoatom2(i2)=0.d0
!! LAGRANGE INTERPOLATION
! do i1=1,mesh1
! do j1=1,mesh1
! if (i1.ne.j1) then
! dr2=rr2(i2)-rr1(j1)
! dr1=rr1(i1)-rr1(j1)
! dfact=dr2/dr1
!!write(*,*)"dr1,dr2",dr1,dr2
! rhoatom2(i2)=rhoatom2(i2)+rhoatom1(i1)*dfact
! end if
! end do
! end do
!! SPLINE INTERPOLATION
do i1=2,mesh1
do i2=1,mesh2
if((rr2(i2).le.rr1(i1)).and.(rr2(i2).ge.rr1(i1-1))) then
r1=rr1(i1-1)
r2=rr1(i1)
if((i1-1).eq.1) then
s1=rhoatom1(i1-1)/rr1(i1-1)
else
dy=rhoatom1(i1-1)-rhoatom1(i1-2)
dx=rr1(i1-1)-rr1(i1-2)
s1=dy/dx
end if
dy=rhoatom1(i1)-rhoatom1(i1-1)
dx=rr1(i1)-rr1(i1-1)
s2=dy/dx
dx1=rr1(i1)-rr1(i1-1)
dy1=rhoatom1(i1)-rhoatom1(i1-1)
t=(rr2(i2)-rr1(i1))/dx1
a=s1*dx1-dy1
b=-s2*dx1+dy1
rhoatom2(i2)=(1.d0-t)*rhoatom1(i1-1)+t*rhoatom1(i1)+t*(1.d0-t)*(a*(1.d0-t)+b*t)
end if
end do
end do
end subroutine interpolate