TABLE OF CONTENTS
ABINIT/psp8nl [ Functions ]
NAME
psp8nl
FUNCTION
Make Kleinman-Bylander/Bloechl form factors f_ln(q) for each projector n for each angular momentum l excepting an l corresponding to the local potential. Note that an arbitrary local potential can be used, so all l from 0 to lmax may be represented.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, FrD, GZ) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
amesh=grid spacing for uniform (linear) radial grid
indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0)
or i=lmn (if useylm=1)
lmax=maximum ang momentum for which nonlocal form factor is desired.
lmax <= 2 allowed.
lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
=if useylm=0, max number of (l,n) comp. over all type of psps
lnmax=max. number of (l,n) components over all type of psps
mmax=number of radial grid points for atomic grid
mqgrid=number of grid points for q grid
pspso=spin-orbit characteristics, govern the content of ffspl and ekb
if =0 : this input requires NO spin-orbit characteristics of the psp
if =2 : this input requires HGH or psp8 characteristics of the psp
if =3 : this input requires HFN characteristics of the psp
qgrid(mqgrid)=values at which form factors are returned
rad(mmax)=radial grid values
vpspll(mmax,lnmax)=nonlocal projectors for each (l,n) on linear
radial grid. Here, these are the product of the reference
wave functions and (v(l,n)-vloc), calculated in the psp generation
program and normalized so that integral(0,rc(l)) vpsll^2 dr = 1,
which leads to the the usual convention for the energies ekb(l,n)
also calculated in the psp generation program.
OUTPUT
ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_ln(q) and second derivative from spline fit for each (l,n).
NOTES
u_l(r) is reference state wavefunction (input as wf);
j_l(q) is a spherical Bessel function;
dV_l(r) = vpsp_l(r)-vloc(r) for angular momentum l;
f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) dV_l(r) r dr]/\sqrt{dvms}$
where dvms = $\int_0^{rmax} [(u_l(r) dV_l(r))^2 dr]$ is the mean
square value of the nonlocal correction for angular momentum l.
Xavier Gonze s E_KB = $ dvms/\int_0^{rmax}[(u_l(r))^2 dV_l(r) dr]$.
This is the eigenvalue of the Kleinman-Bylander operator and sets
the energy scale of the nonlocal psp corrections.
PARENTS
psp8in,psp9in
CHILDREN
ctrap,sbf8,spline
SOURCE
67 #if defined HAVE_CONFIG_H 68 #include "config.h" 69 #endif 70 71 #include "abi_common.h" 72 73 subroutine psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,& 74 & mqgrid,qgrid,rad,vpspll) 75 76 use defs_basis 77 use m_splines 78 use m_profiling_abi 79 80 use m_special_funcs, only : sbf8 81 use m_numeric_tools, only : ctrap 82 83 !This section has been created automatically by the script Abilint (TD). 84 !Do not modify the following lines by hand. 85 #undef ABI_FUNC 86 #define ABI_FUNC 'psp8nl' 87 !End of the abilint section 88 89 implicit none 90 91 !Arguments---------------------------------------------------------- 92 !scalars 93 integer,intent(in) :: lmax,lmnmax,lnmax,mmax,mqgrid 94 real(dp),intent(in) :: amesh 95 !arrays 96 integer,intent(in) :: indlmn(6,lmnmax) 97 real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vpspll(mmax,lnmax) 98 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i 99 100 !Local variables------------------------------- 101 !Following parameter controls accuracy of Fourier transform based on qmax 102 !and represents the minimun number of integration points in one period. 103 !scalars 104 integer,parameter :: NPT_IN_2PI=200 105 integer :: iln,iln0,ilmn,iq,ir,irmu,irn,ll,mesh_mult,mmax_new,mvpspll 106 real(dp) :: amesh_new,arg,c1,c2,c3,c4,dri,qmesh,result,tv,xp,xpm1,xpm2,xpp1 107 real(dp) :: yp1,ypn 108 !arrays 109 real(dp) :: sb_out(4) 110 real(dp),allocatable :: rad_new(:),vpspll_new(:,:),work(:,:),work2(:) 111 112 ! ************************************************************************* 113 114 !Find r mesh spacing necessary for accurate integration at qmax 115 amesh_new=2.d0*pi/(NPT_IN_2PI*qgrid(mqgrid)) 116 117 !Choose submultiple of input mesh 118 mesh_mult=int(amesh/amesh_new) + 1 119 mmax_new=mesh_mult*(mmax-1)+1 120 amesh_new=amesh/dble(mesh_mult) 121 122 ABI_ALLOCATE(rad_new,(mmax_new)) 123 ABI_ALLOCATE(vpspll_new,(mmax_new,lnmax)) 124 125 if(mesh_mult==1) then 126 rad_new(:)=rad(:) 127 else 128 ! Set up new radial mesh 129 irn=1 130 do ir=1,mmax-1 131 do irmu=0,mesh_mult-1 132 rad_new(irn)=rad(ir)+dble(irmu)*amesh_new 133 irn=irn+1 134 end do 135 end do 136 rad_new(mmax_new)=rad(mmax) 137 end if 138 139 !Interpolate projectors onto new grid if called for 140 !Cubic polynomial interpolation is used which is consistent 141 !with the original interpolation of these functions from 142 !a log grid to the input linear grid. 143 dri = one/amesh 144 do irn=1,mmax_new 145 ! index to find bracketing input mesh points 146 if(mesh_mult>1) then 147 ir = irn/mesh_mult + 1 148 ir = max(ir,2) 149 ir = min(ir,mmax-2) 150 ! interpolation coefficients 151 xp = dri * (rad_new(irn) - rad(ir)) 152 xpp1 = xp + one 153 xpm1 = xp - one 154 xpm2 = xp - two 155 c1 = -xp * xpm1 * xpm2 * sixth 156 c2 = xpp1 * xpm1 * xpm2 * half 157 c3 = - xp * xpp1 * xpm2 * half 158 c4 = xp * xpp1 * xpm1 * sixth 159 ! Now do the interpolation on all projectors for this grid point 160 161 iln0=0 162 do ilmn=1,lmnmax 163 iln=indlmn(5,ilmn) 164 if (iln>iln0) then 165 iln0=iln 166 tv = c1 * vpspll(ir - 1, iln) & 167 & + c2 * vpspll(ir , iln) & 168 & + c3 * vpspll(ir + 1, iln) & 169 & + c4 * vpspll(ir + 2, iln) 170 if(abs(tv)>tol10) then 171 vpspll_new(irn,iln)=tv 172 mvpspll=irn 173 else 174 vpspll_new(irn,iln)=zero 175 end if 176 end if 177 end do 178 179 else 180 ! With no mesh multiplication, just copy projectors 181 ir=irn 182 iln0=0 183 do ilmn=1,lmnmax 184 iln=indlmn(5,ilmn) 185 if (iln>iln0) then 186 iln0=iln 187 tv = vpspll(ir,iln) 188 if(abs(tv)>tol10) then 189 vpspll_new(irn,iln)=tv 190 mvpspll=irn 191 else 192 vpspll_new(irn,iln)=zero 193 end if 194 end if 195 end do 196 197 end if 198 end do !irn 199 200 ABI_ALLOCATE(work,(mvpspll,lnmax)) 201 202 !Loop over q values 203 do iq=1,mqgrid 204 arg=2.d0*pi*qgrid(iq) 205 206 ! Set up integrands 207 do ir=1,mvpspll 208 call sbf8(lmax+1,arg*rad_new(ir),sb_out) 209 iln0=0 210 do ilmn=1,lmnmax 211 iln=indlmn(5,ilmn) 212 if (iln>iln0) then 213 iln0=iln 214 ll=indlmn(1,ilmn) 215 work(ir,iln)=sb_out(ll+1)*vpspll_new(ir,iln)*rad_new(ir) 216 end if 217 end do 218 end do !ir 219 220 ! Do integral from zero to rad_new(mvpspll) 221 iln0=0 222 do ilmn=1,lmnmax 223 iln=indlmn(5,ilmn) 224 if (iln>iln0) then 225 iln0=iln 226 call ctrap(mvpspll,work(1,iln),amesh_new,result) 227 ffspl(iq,1,iln)=result 228 end if 229 end do 230 231 ! End loop over q mesh 232 end do !iq 233 234 !Fit splines for form factors 235 ABI_ALLOCATE(work2,(mqgrid)) 236 qmesh=qgrid(2)-qgrid(1) 237 238 iln0=0 239 do ilmn=1,lmnmax 240 iln=indlmn(5,ilmn) 241 if (iln>iln0) then 242 iln0=iln 243 ! Compute derivatives of form factors at ends of interval 244 yp1=(-50.d0*ffspl(1,1,iln)+96.d0*ffspl(2,1,iln)-72.d0*ffspl(3,1,iln)& 245 & +32.d0*ffspl(4,1,iln)- 6.d0*ffspl(5,1,iln))/(24.d0*qmesh) 246 ypn=(6.d0*ffspl(mqgrid-4,1,iln)-32.d0*ffspl(mqgrid-3,1,iln)& 247 & +72.d0*ffspl(mqgrid-2,1,iln)-96.d0*ffspl(mqgrid-1,1,iln)& 248 & +50.d0*ffspl(mqgrid,1,iln))/(24.d0*qmesh) 249 250 call spline(qgrid,ffspl(1,1,iln),mqgrid,yp1,ypn,ffspl(1,2,iln)) 251 end if 252 end do 253 254 ABI_DEALLOCATE(rad_new) 255 ABI_DEALLOCATE(vpspll_new) 256 ABI_DEALLOCATE(work) 257 ABI_DEALLOCATE(work2) 258 259 end subroutine psp8nl