TABLE OF CONTENTS
ABINIT/der_int [ Functions ]
NAME
der_int
FUNCTION
Given input function f(i) on Teter radial grid, and grid spacing dr(i), compute function derivative df/dr on points from 0 to n. Integrate function f(i) on grid r(i) from r(0) to r(nlast). Note that array dimensions start at 0.
INPUTS
f(0 to nlast)=function values on grid r(0 to nlast)=radial grid points dr(0 to nlast)=INVERSE of spacing on grid nlast=radial grid point for upper limit
OUTPUT
df(0 to n)=derivative $ \frac{df}{dr}$ on grid
smf= $ \int_{r(0)}^{r(nlast)} f(r) dr $.
PARENTS
psp1lo,psp1nl
CHILDREN
SOURCE
944 subroutine der_int(ff,df,rr,dr,nlast,smf) 945 946 use defs_basis 947 use m_errors 948 use m_profiling_abi 949 950 !This section has been created automatically by the script Abilint (TD). 951 !Do not modify the following lines by hand. 952 #undef ABI_FUNC 953 #define ABI_FUNC 'der_int' 954 !End of the abilint section 955 956 implicit none 957 958 !Arguments ------------------------------------ 959 !nmax sets standard number of grid points ! SHOULD BE REMOVED 960 !scalars 961 integer,parameter :: nmax=2000 962 integer,intent(in) :: nlast 963 real(dp),intent(out) :: smf 964 !no_abirules 965 !Note that dimension here starts at 0 966 real(dp), intent(in) :: dr(0:nmax),ff(0:nmax),rr(0:nmax) 967 real(dp), intent(out) :: df(0:nmax) 968 969 !Local variables------------------------------- 970 !scalars 971 integer :: jj 972 real(dp),parameter :: div12=1.d0/12.d0 973 real(dp) :: hh 974 character(len=500) :: message 975 976 ! ************************************************************************* 977 978 !Check that nlast lie within 0 to nmax 979 if (nlast<0.or.nlast>nmax) then 980 write(message, '(a,i12,a,i12)' )& 981 & ' nlast=',nlast,' lies outside range [0,nmax] with dimension nmax=',nmax 982 MSG_BUG(message) 983 end if 984 985 !Compute derivatives at lower end, near r=0 986 df(0)=-25.d0/12.d0*ff(0)+4.d0*ff(1)-3.d0*ff(2)+4.d0/3.d0*ff(3)& 987 & -1.d0/4.d0*ff(4) 988 df(1)=-1.d0/4.d0*ff(0)-5.d0/6.d0*ff(1)+3.d0/2.d0*ff(2)& 989 & -1.d0/2.d0*ff(3)+1.d0/12.d0*ff(4) 990 991 !Run over range from just past r=0 to near r(n), using central differences 992 do jj=2,nlast-2 993 df(jj)=(ff(jj-2)-8.d0*(ff(jj-1)-ff(jj+1))-ff(jj+2))*div12 994 end do 995 996 !Compute derivative at upper end of range 997 if (nlast < 4) then 998 message = ' der_int: ff does not have enough elements. nlast is too low' 999 MSG_ERROR(message) 1000 end if 1001 1002 df(nlast-1)=-1.d0/12.d0*ff(nlast-4)& 1003 & +1.d0/2.d0*ff(nlast-3)& 1004 & -3.d0/2.d0*ff(nlast-2)& 1005 & +5.d0/6.d0*ff(nlast-1)& 1006 & +1.d0/4.d0*ff(nlast) 1007 df(nlast)=1.d0/4.d0*ff(nlast-4)& 1008 & -4.d0/3.d0*ff(nlast-3)& 1009 & +3.d0*ff(nlast-2)& 1010 & -4.d0*ff(nlast-1)& 1011 & +25.d0/12.d0*ff(nlast) 1012 1013 !Apply correct normalization over full range 1014 do jj=0,nlast 1015 df(jj)=df(jj)*dr(jj) 1016 end do 1017 1018 smf=0.d0 1019 do jj=0,nlast-1 1020 hh=rr(jj+1)-rr(jj) 1021 smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1))) 1022 end do 1023 smf=smf/12.d0 1024 1025 end subroutine der_int
ABINIT/psp1in [ Functions ]
NAME
psp1in
FUNCTION
Initialize pspcod=1 or 4 pseudopotential (Teter format): continue to read the corresponding file, then compute the local and non-local potentials.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FrD, MT) 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
dq= spacing of the q-grid
lloc=angular momentum choice of local pseudopotential
lmax=value of lmax mentioned at the second line of the psp file
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=maximum number of points in real space grid in the psp file
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mqgrid=dimension of q (or G) grid for arrays.
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
pspcod=pseudopotential type
qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
useylm=governs the way the nonlocal operator is to be applied:
1=using Ylm, 0=using Legendre polynomials
zion=nominal valence of atom as specified in psp file
znucl=atomic number of atom as specified in psp file
OUTPUT
ekb(lnmax)=Kleinman-Bylander energy, {{\ \begin{equation} \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]} {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]} \end{equation} }} for each (l,n) ekb1(mpsang)= Kleinman-Bylander energy from the psp file, for iproj=1 ekb2(mpsang)= Kleinman-Bylander energy from the psp file, for iproj=2 epsatm=$(4\pi) \int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree) epspsp(mpsang)=values of epsatm for different angular momenta, from the psp file e990(mpsang)=ecut at which 0.99 of the kinetic energy is recovered e999(mpsang)=ecut at which 0.999 of the kinetic energy is recovered ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) nproj(mpsang)=number of projection functions for each angular momentum qchrg is the total (integrated) core charge rcpsp(mpsang)=cut-off radius for each angular momentum rms(mpsang)=root mean square of the KB psp vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file xcccrc=XC core correction cutoff radius (bohr)
NOTES
there are only minor differences in the two formats
1) With pspcod=1, even for the LOCAL angular momentum, there is
a block for the wfs (can be set to zero, though)
2) The core charge density differs: for pspcod=1, it is a
revised expression for core density of 5 Nov 1992, while
for pspcod=4, it is an older expression, of 7 May 1992 .
PARENTS
pspatm
CHILDREN
psp1cc,psp1lo,psp1nl,psp4cc,spline,wrtout
SOURCE
78 #if defined HAVE_CONFIG_H 79 #include "config.h" 80 #endif 81 82 #include "abi_common.h" 83 84 subroutine psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,& 85 & e990,e999,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 86 & mmax,mpsang,mqgrid,nproj,n1xccc,pspcod,& 87 & qchrg,qgrid,rcpsp,rms,useylm,vlspl,xcccrc,xccc1d,& 88 & zion,znucl) 89 90 use defs_basis 91 use m_errors 92 use m_profiling_abi 93 use m_splines 94 95 !This section has been created automatically by the script Abilint (TD). 96 !Do not modify the following lines by hand. 97 #undef ABI_FUNC 98 #define ABI_FUNC 'psp1in' 99 use interfaces_14_hidewrite 100 use interfaces_64_psp, except_this_one => psp1in 101 !End of the abilint section 102 103 implicit none 104 105 !Arguments ------------------------------------ 106 !scalars 107 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc,pspcod 108 integer,intent(in) :: useylm 109 real(dp),intent(in) :: dq,zion,znucl 110 real(dp),intent(out) :: epsatm,qchrg,xcccrc 111 !arrays 112 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang) 113 real(dp),intent(in) :: qgrid(mqgrid) 114 real(dp),intent(out) :: e990(mpsang),e999(mpsang),ekb(lnmax),ekb1(mpsang) 115 real(dp),intent(out) :: ekb2(mpsang),epspsp(mpsang) 116 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) 117 real(dp),intent(out) :: rcpsp(mpsang),rms(mpsang),vlspl(mqgrid,2) 118 real(dp),intent(inout) :: xccc1d(n1xccc,6) 119 120 !Local variables------------------------------- 121 !scalars 122 integer :: ii,iln,index,ipsang,kk,lhigh,ll,mm,nlmax 123 real(dp) :: arg,dq2pi,fchrg,rchrg,xx,yp1,ypn 124 character(len=500) :: message,errmsg 125 !arrays 126 real(dp),allocatable :: drad(:),ekb_tmp(:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:) 127 real(dp),allocatable :: vpspll(:,:),wfll(:,:),wksincos(:,:,:),work_space(:) 128 real(dp),allocatable :: work_spl1(:),work_spl2(:) 129 130 ! *************************************************************************** 131 132 !Note: Teter s grid is hard-coded at mmax=2001 133 !mmax was read from the pseudopotential file in the calling routine 134 if (mmax/=2001) then 135 write(message, '(a,i12,a,a,a,a)' )& 136 & 'Using Teter grid (pspcod=1 and 4) but mmax=',mmax,ch10,& 137 & 'mmax must be 2001 for Teter grid.',ch10,& 138 & 'Action : check your pseudopotential input file.' 139 MSG_ERROR(message) 140 end if 141 142 !File format of formatted Teter psp input (the 3 first lines 143 !have already been read in calling -pspatm- routine) : 144 145 !(1) title (character) line 146 !(2) znucl,zion,pspdat 147 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well 148 !For each angular momentum : 149 !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll) 150 !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll) 151 !(6) rchrg,fchrg,qchrg 152 !(7) ll 153 !(8) (vpsp(j,ll),j=0,nmax) 154 !Then for iproj=1 to 2 155 !for ll=0,lmax 156 !(10) ll 157 !(11) ((upsp(j,ll,iproj),j=0,nmax) 158 159 do ipsang=1,lmax+1 160 161 read (tmp_unit,*,err=10,iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang) 162 write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) & 163 & ipsang-1,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),& 164 & 'l,e99.0,e99.9,nproj,rcpsp' 165 call wrtout(ab_out,message,'COLL') 166 call wrtout(std_out, message,'COLL') 167 168 read (tmp_unit,*,err=10,iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang) 169 write(message, '(4f13.8,t55,a)' ) & 170 & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),& 171 & ' rms, ekb1, ekb2, epsatm' 172 call wrtout(ab_out,message,'COLL') 173 call wrtout(std_out, message,'COLL') 174 175 end do 176 177 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn 178 index=0;iln=0;indlmn(:,:)=0 179 do ipsang=1,lmax+1 180 if(nproj(ipsang)>0)then 181 ll=ipsang-1 182 do kk=1,nproj(ipsang) 183 iln=iln+1 184 do mm=1,2*ll*useylm+1 185 index=index+1 186 indlmn(1,index)=ll 187 indlmn(2,index)=mm-ll*useylm-1 188 indlmn(3,index)=kk 189 indlmn(4,index)=ll*ll+(1-useylm)*ll+mm 190 indlmn(5,index)=iln 191 indlmn(6,index)=1 192 end do 193 end do 194 end if 195 end do 196 197 read (tmp_unit,*,err=10,iomsg=errmsg) rchrg,fchrg,qchrg 198 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,& 199 & 'rchrg,fchrg,qchrg' 200 call wrtout(ab_out,message,'COLL') 201 call wrtout(std_out, message,'COLL') 202 203 !Generate core charge function and derivatives, if needed 204 if(fchrg>1.0d-15)then 205 if(pspcod==1)then 206 call psp1cc(fchrg,n1xccc,xccc1d) 207 ! The core charge function for pspcod=1 208 ! becomes zero beyond 3*rchrg only. Thus xcccrc must be set 209 ! equal to 3*rchrg . 210 xcccrc=3*rchrg 211 else if(pspcod==4)then 212 call psp4cc(fchrg,n1xccc,xccc1d) 213 ! For pspcod=4, the core charge cut off exactly beyond rchrg 214 xcccrc=rchrg 215 end if 216 else 217 xcccrc=0.0d0 218 xccc1d(:,:)=0.0d0 219 end if 220 221 !-------------------------------------------------------------------- 222 !Will now proceed at the reading of pots and wfs, as well as their treatment 223 224 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials 225 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg. 226 !rad(:)=radial grid r(i) 227 !drad(:)= inverse of d(r(i))/d(i) for radial grid 228 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions 229 230 ABI_ALLOCATE(vloc,(mmax)) 231 ABI_ALLOCATE(vpspll,(mmax,mpsang)) 232 if(lmax==-1) vpspll(:,:)=zero 233 234 !(1) Read atomic pseudopotential for each l, filling up array vpspll 235 !Note: put each l into vpspll(:,l+1) 236 do ipsang=1,lmax+1 237 read (tmp_unit,*,err=10,iomsg=errmsg) ll 238 read (tmp_unit,*,err=10,iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax) 239 ! DEBUG 240 ! write(std_out,*) 'END OF READING PSP',ll,'OK' 241 ! ENDDEBUG 242 243 end do 244 245 !Copy appropriate nonlocal psp for use as local one 246 vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 ) 247 248 !DEBUG 249 !write(std_out,*) 'VLOC=',vloc(1),vloc(2),vloc(3) 250 !write(std_out,*) 'VLOC=',vloc(4),vloc(5),vloc(6) 251 !ENDDEBUG 252 253 !(2) Create radial grid, and associated quantities 254 255 ABI_ALLOCATE(rad,(mmax)) 256 ABI_ALLOCATE(drad,(mmax)) 257 ABI_ALLOCATE(wksincos,(mmax,2,2)) 258 259 !Teter grid--need both r and dr in this case 260 do ii=0,mmax-1 261 xx=dble(ii)/dble(mmax-1) 262 rad (ii+1)=100.d0*(xx+.01d0)**5-1.d-8 263 drad(ii+1)=500.d0*(xx+.01d0)**4/dble(mmax-1) 264 end do 265 266 !DEBUG 267 !write(std_out,*) 'RADIAL GRID CREATED' 268 !ENDDEBUG 269 270 !here compute sin(r(:)*dq) and cos(r(:)*dq) 271 !NOTE : also invert dr !! 272 dq2pi=2.0d0*pi*dq 273 do ii=1,mmax 274 arg=dq2pi*rad(ii) 275 drad(ii)=1.0d0/drad(ii) 276 wksincos(ii,1,1)=sin(arg) 277 wksincos(ii,2,1)=cos(arg) 278 end do 279 280 !(3)Carry out calculations for local (lloc) pseudopotential. 281 !Obtain Fourier transform (1-d sine transform) 282 !to get q^2 V(q). 283 ABI_ALLOCATE(work_space,(mqgrid)) 284 ABI_ALLOCATE(work_spl1,(mqgrid)) 285 ABI_ALLOCATE(work_spl2,(mqgrid)) 286 call psp1lo(drad,epsatm,mmax,mqgrid,qgrid,& 287 & work_spl1,rad,vloc,wksincos,yp1,ypn,zion) 288 289 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 290 call spline (qgrid,work_spl1,mqgrid,yp1,ypn,work_spl2) 291 vlspl(:,1)=work_spl1(:) 292 vlspl(:,2)=work_spl2(:) 293 294 ABI_DEALLOCATE(work_space) 295 ABI_DEALLOCATE(work_spl1) 296 ABI_DEALLOCATE(work_spl2) 297 298 !(4)Take care of non-local part 299 300 !Zero out all Kleinman-Bylander energies to initialize 301 ekb(:)=0.0d0 302 303 !DEBUG 304 !write(std_out,*)' psp1in : before nonlocal corrections ' 305 !write(std_out,*)' psp1in : lloc, lmax = ',lloc,lmax 306 !if(.true.)stop 307 !ENDDEBUG 308 309 !Allow for option of no nonlocal corrections (lloc=lmax=0) 310 if (lloc==0.and.lmax==0) then 311 312 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 313 call wrtout(ab_out,message,'COLL') 314 call wrtout(std_out, message,'COLL') 315 316 else 317 318 ! Proceed to make Kleinman-Bylander form factors for 319 ! each l up to lmax 320 321 ! Read wavefunctions for each l up to lmax 322 ABI_ALLOCATE(wfll,(mmax,mpsang)) 323 do ipsang=1,lmax+1 324 ! For pspcod==4, wfs for the local angular momentum are not written 325 if (nproj(ipsang)/=0 .or. pspcod==1) then 326 read (tmp_unit,*,err=10,iomsg=errmsg) ll 327 if (ipsang/=ll+1) then 328 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 329 & 'Pseudopotential input file does not have',ch10,& 330 & 'angular momenta in order expected for first projection',& 331 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 332 & 'Action : check your pseudopotential input file.' 333 MSG_ERROR(message) 334 end if 335 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 336 337 else 338 wfll(:,ipsang)=0.0d0 339 end if 340 341 end do 342 ! ---------------------------------------------------------------------- 343 ! Compute KB form factors and fit splines 344 345 ! nlmax is highest l for which a nonlocal correction is being computed 346 nlmax=lmax 347 if (lloc==lmax) nlmax=lmax-1 348 349 ! DEBUG 350 ! write(std_out,*)' psp1in : lmax,lloc=',lmax,lloc 351 ! ENDDEBUG 352 ABI_ALLOCATE(ekb_tmp,(mpsang,2)) 353 ABI_ALLOCATE(ffspl_tmp,(mqgrid,2,nlmax+1,2)) 354 355 call psp1nl(drad,ekb_tmp(:,1),ffspl_tmp(:,:,:,1),lloc,& 356 & nlmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos) 357 358 ! Read second wavefunction for second projection operator 359 ! (only read cases where nproj(ll)=2) 360 ! --also find highest l for which nproj(l)=2 361 362 lhigh=-1 363 do ipsang=1,min(lmax+1,mpsang) 364 if (nproj(ipsang)==2) then 365 lhigh=ipsang-1 366 read (tmp_unit,*,err=10,iomsg=errmsg) ll 367 if (ipsang/=ll+1) then 368 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 369 & 'Pseudopotential input file does not have',ch10,& 370 & 'angular momenta in order expected for second projection',& 371 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 372 & 'Action : check your pseudopotential input file.' 373 MSG_ERROR(message) 374 end if 375 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 376 377 else 378 wfll(:,ipsang)=0.0d0 379 380 end if 381 end do 382 383 ! Compute KB form factors and fit splines for second wf if any 384 385 if (lhigh>-1) then 386 call psp1nl(drad,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lloc,& 387 & lhigh,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos) 388 end if 389 390 ! Convert ekb and ffspl 391 iln=0 392 do ii=1,lmnmax 393 kk=indlmn(5,ii) 394 if (kk>iln) then 395 iln=kk 396 ekb(kk)=ekb_tmp(1+indlmn(1,ii),indlmn(3,ii)) 397 ! DEBUG 398 ! write(std_out,*)' psp1in : lmnmax,ii,indlmn(1,ii)=',lmnmax,ii,indlmn(1,ii) 399 ! ENDDEBUG 400 ffspl(:,:,kk)=ffspl_tmp(:,:,1+indlmn(1,ii),indlmn(3,ii)) 401 end if 402 end do 403 404 ABI_DEALLOCATE(ekb_tmp) 405 ABI_DEALLOCATE(ffspl_tmp) 406 ABI_DEALLOCATE(wfll) 407 408 end if 409 410 ABI_DEALLOCATE(vpspll) 411 ABI_DEALLOCATE(rad) 412 ABI_DEALLOCATE(drad) 413 ABI_DEALLOCATE(vloc) 414 ABI_DEALLOCATE(wksincos) 415 416 return 417 418 ! Handle IO error 419 10 continue 420 MSG_ERROR(errmsg) 421 422 end subroutine psp1in
ABINIT/psp1lo [ Functions ]
NAME
psp1lo
FUNCTION
Compute sine transform to transform from v(r) to q^2 v(q) using subroutines related to Teter atomic structure grid.
INPUTS
drad(mmax)=inverse of r grid spacing at each point
mmax=number of radial r grid points (Teter grid)
mqgrid=number of grid points in q from 0 to qmax.
qgrid(mqgrid)=q grid values (bohr**-1).
rad(mmax)=r grid values (bohr).
vloc(mmax)=v(r) on radial grid.
wksincos(mmax,2,2)=contains sine and cosine of 2*pi*r(:)*dq and 2*pi*r(:)*q
at input : wksincos(:,1,1)=sine of 2*pi*r(:)*dq
wksincos(:,2,1)=cosine of 2*pi*r(:)*dq
wksincos(:,:,2) is not initialized, will be used inside the routine
zion=nominal valence charge of atom.
OUTPUT
epsatm= $4\pi \int[r^2 (v(r)+Zv/r) dr]$
q2vq(mqgrid)=$q^2 v(q)$
=$\displaystyle -Zv/\pi+q^2 4\pi\int(\frac{\sin(2\pi q r)}{2 \pi q r})(r^2 v(r)+r Zv)dr$.
yp1,ypn=derivative of q^2 v(q) wrt q at q=0 and q=qmax
(needed for spline fitter).
PARENTS
psp1in
CHILDREN
der_int,sincos
SOURCE
461 subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,& 462 & vloc,wksincos,yp1,ypn,zion) 463 464 use defs_basis 465 use m_profiling_abi 466 467 !This section has been created automatically by the script Abilint (TD). 468 !Do not modify the following lines by hand. 469 #undef ABI_FUNC 470 #define ABI_FUNC 'psp1lo' 471 use interfaces_64_psp, except_this_one => psp1lo 472 !End of the abilint section 473 474 implicit none 475 476 !Arguments ------------------------------------ 477 !scalars 478 integer,intent(in) :: mmax,mqgrid 479 real(dp),intent(in) :: zion 480 real(dp),intent(out) :: epsatm,yp1,ypn 481 !arrays 482 real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 483 real(dp),intent(inout) :: wksincos(mmax,2,2) 484 real(dp),intent(out) :: q2vq(mqgrid) 485 486 !Local variables------------------------------- 487 !scalars 488 integer,parameter :: mma0=2001 489 integer :: iq,ir,irmax 490 real(dp),parameter :: scale=10.0d0 491 real(dp) :: result,test,tpiq 492 !arrays 493 real(dp) :: wk(mma0),wk1(mma0),wk2(mma0) 494 495 ! ************************************************************************* 496 497 !DEBUG 498 !write(std_out,*)' psp1lo : enter ' 499 !if(.true.)stop 500 !ENDDEBUG 501 502 !Do q=0 separately (compute epsatm) 503 !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr] 504 !Treat r=0 by itself 505 wk(1)=0.0d0 506 507 do ir=2,mmax 508 ! (at large r do not want prefactor of r^2 and should see 509 ! V(r)+Zv/r go to 0 at large r) 510 test=vloc(ir)+zion/rad(ir) 511 ! DEBUG 512 ! write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test 513 ! ENDDEBUG 514 ! In this routine, NO cut-off radius is imposed : the input 515 ! vloc MUST be in real(dp) to obtain numerically 516 ! accurate values. The error can be on the order of 0.001 Ha ! 517 if (abs(test)<1.0d-20) then 518 wk(ir)=0.0d0 519 else 520 wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion) 521 end if 522 end do 523 !Do integral from 0 to r(max) (disregard contrib beyond r(max) 524 !(need numerical derivatives to do integral) 525 !Use mmax-1 to convert to Teter s dimensioning starting at 0 526 call der_int(wk,wk2,rad,drad,mmax-1,result) 527 528 !DEBUG 529 !write(std_out,*)' psp1lo : result ',result 530 !stop 531 !ENDDEBUG 532 533 epsatm=4.d0*pi*(result) 534 !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi 535 q2vq(1)=-zion/pi 536 537 !Prepare loop over q values 538 irmax=mmax+1 539 do ir=mmax,2,-1 540 test=vloc(ir)+zion/rad(ir) 541 wk1(ir)=test*rad(ir) 542 ! Will ignore tail within decade of machine precision 543 if ((scale+abs(test))==scale .and. irmax==ir+1) then 544 irmax=ir 545 end if 546 end do 547 !Increase irmax a bit : this is copied from psp1nl 548 irmax=irmax+4 549 if(irmax>mmax)irmax=mmax 550 551 !Loop over q values 552 do iq=2,mqgrid 553 tpiq=two_pi*qgrid(iq) 554 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 555 ! set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral 556 !$\displaystyle -Zv/\pi + q^2 4\pi \int[\frac{\sin(2\pi q r)}{2\pi q r}(r^2 v(r)+r Zv)dr]$. 557 ! Handle r=0 separately 558 wk(1)=0.0d0 559 do ir=2,irmax 560 wk(ir)=wksincos(ir,1,2)*wk1(ir) 561 end do 562 ! do integral from 0 to r(max) 563 if(irmax>mmax-1)irmax=mmax-1 564 565 call der_int(wk,wk2,rad,drad,irmax,result) 566 ! store q^2 v(q) 567 q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result 568 end do 569 570 !Compute derivatives of q^2 v(q) at ends of interval 571 yp1=0.0d0 572 !ypn=$\displaystyle 2\int_0^\infty (\sin (2\pi qmax r)+(2\pi qmax r)\cos (2\pi qmax r)(r V(r)+Z)dr]$ 573 !integral from r(mmax) to infinity is overkill; ignore 574 !set up integrand 575 !Handle r=0 separately 576 wk(1)=0.0d0 577 tpiq=two_pi*qgrid(mqgrid) 578 do ir=2,mmax 579 test=vloc(ir)+zion/rad(ir) 580 ! Ignore contributions within decade of machine precision 581 if ((scale+abs(test))==scale) then 582 wk(ir)=0.0d0 583 else 584 wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * & 585 & (rad(ir)*vloc(ir)+zion) 586 end if 587 end do 588 call der_int(wk,wk2,rad,drad,mmax-1,result) 589 590 ypn=2.0d0*result 591 592 end subroutine psp1lo
ABINIT/psp1nl [ Functions ]
NAME
psp1nl
FUNCTION
Make Kleinman-Bylander form factors f_l(q) for each l from 0 to lmax; Vloc is assumed local potential.
INPUTS
dr(mmax)=inverse of grid spacing for radial grid
lloc=angular momentum of local channel (avoid doing integrals for this l)
lmax=maximum ang momentum for which nonlocal form factor is desired.
mmax=number of radial grid points for atomic grid
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mqgrid=number of grid points for q grid
qgrid(mqgrid)=values at which form factors are returned
rad(mmax)=radial grid values
vloc(mmax)=local pseudopotential on radial grid
vpspll(mmax,lmax+1)=nonlocal pseudopotentials for each l on radial grid
wfll(mmax,lmax+1)=reference state wavefunctions on radial grid
wksincos(mmax,2,2)=contains sine and cosine of 2*pi*r(:)*dq and 2*pi*r(:)*q
at input : wksincos(:,1,1)=sine of 2*pi*r(:)*dq
wksincos(:,2,1)=cosine of 2*pi*r(:)*dq
wksincos(:,:,2) is not initialized, will be used inside the routine
OUTPUT
ekb(mpsang)=Kleinman-Bylander energy, {{\ \begin{equation} \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]} {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]} \end{equation} }} for each l ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum
NOTES
u_l(r) is reference state wavefunction (input as wfll);
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=$\displaystyle \int_0^{rmax}[(u_l(r) dV_l(r))^2 dr]$ is the mean
square value of the nonlocal correction for angular momentum l.
E_KB = $\displaystyle \frac{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.
Bessel functions replaced by besj, which accomodates args near 0.
PARENTS
psp1in
CHILDREN
besjm,der_int,sincos,spline
SOURCE
650 subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,& 651 & qgrid,rad,vloc,vpspll,wfll,wksincos) 652 653 use defs_basis 654 use m_errors 655 use m_profiling_abi 656 use m_splines 657 658 use m_special_funcs, only : besjm 659 660 !This section has been created automatically by the script Abilint (TD). 661 !Do not modify the following lines by hand. 662 #undef ABI_FUNC 663 #define ABI_FUNC 'psp1nl' 664 use interfaces_64_psp, except_this_one => psp1nl 665 !End of the abilint section 666 667 implicit none 668 669 !Arguments ------------------------------------ 670 !scalars 671 integer,intent(in) :: lloc,lmax,mmax,mpsang,mqgrid 672 !arrays 673 real(dp),intent(in) :: dr(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 674 real(dp),intent(in) :: vpspll(mmax,mpsang),wfll(mmax,mpsang) 675 real(dp),intent(inout) :: wksincos(mmax,2,2) 676 real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang) 677 678 !Local variables------------------------------- 679 !scalars 680 integer,parameter :: dpsang=5 681 integer :: iq,ir,irmax,lp1 682 real(dp) :: dvwf,result,test,tpiq,yp1,ypn 683 character(len=500) :: message 684 !arrays 685 real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang) 686 real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:) 687 real(dp),allocatable :: work_spl(:) 688 689 ! ************************************************************************* 690 691 !DEBUG 692 !write(std_out,*)' psp1nl : enter' 693 !stop 694 !ENDDEBUG 695 696 !Zero out Kleinman-Bylander energies ekb 697 ekb(:)=0.0d0 698 !Zero out eta and other parameters too (so 0 s show up in output later) 699 eta(:)=0.0d0 700 dvms(:)=0.0d0 701 ckb(:)=0.0d0 702 703 !Allow for no nonlocal correction (lmax=-1) 704 if (lmax/=-1) then 705 706 ! Check that lmax is within allowed range 707 if (lmax<0.or.lmax>3) then 708 write(message, '(a,i12,a,a,a,a,a,a,a)' )& 709 & 'lmax=',lmax,' is not an allowed value.',ch10,& 710 & 'Allowed values are -1 for no nonlocal correction or else',ch10,& 711 & '0, 1, 2, or 3 for maximum l nonlocal correction.',ch10,& 712 & 'Action: check the input atomic psp data file for lmax.' 713 MSG_ERROR(message) 714 end if 715 716 ! Compute normalizing integrals eta=<dV> and mean square 717 ! nonlocal psp correction dvms=<dV^2> 718 ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction 719 720 ABI_ALLOCATE(work1,(mmax+1)) 721 ABI_ALLOCATE(work2,(mmax+1)) 722 ABI_ALLOCATE(work_spl,(mqgrid)) 723 ABI_ALLOCATE(work5,(mmax)) 724 ABI_ALLOCATE(besjx,(mmax)) 725 726 do lp1=1,lmax+1 727 728 ! Only do the work if nonlocal correction is nonzero 729 if (lp1 /= lloc+1) then 730 731 ! integrand for 0 to r(mmax) 732 do ir=1,mmax 733 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 734 work1(ir)=wfll(ir,lp1)*dvwf 735 end do 736 737 ! do integral 738 ! first need derivative of function; note use of 739 ! shifted indices to accomodate Mike Teter s choice of 0:mmax-1 740 call der_int(work1,work2,rad,dr,mmax-1,result) 741 eta(lp1)=result 742 743 ! DEBUG 744 ! write(std_out,*)' psp1nl : write eta(lp1)' 745 ! write(std_out,*)result 746 ! do ir=1,mmax,61 747 ! write(std_out,*)vpspll(ir,lp1),vloc(ir),wfll(ir,lp1) 748 ! end do 749 ! write(std_out,*) 750 ! do ir=1,mmax,61 751 ! write(std_out,*)work1(ir),rad(ir),dr(ir) 752 ! end do 753 ! ENDDEBUG 754 755 do ir=1,mmax 756 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 757 work1(ir)=dvwf**2 758 end do 759 call der_int(work1,work2,rad,dr,mmax-1,result) 760 761 dvms(lp1)=result 762 763 ! If dvms is not 0 for any given angular momentum l, 764 ! compute Xavier Gonze s definition of the Kleinman-Bylander 765 ! energy E_KB = dvms/eta. In this case also renormalize 766 ! the projection operator to u_KB(r)=$u_l(r) dV(r)/\sqrt{dvms}$. 767 ! This means dvwf gets multiplied by the normalization factor 768 ! "renorm"=$1/\sqrt{dvms}$ as seen below. 769 ! With dvwf=dV(r)*wf(r) for wf(r)=``radial'' wf, the integrand 770 ! for each angular momentum l is 771 ! Bessel_l(2 $\pi$ q r) * wf(r) * dV(r) * r; 772 ! NOTE presence of extra r in integrand. 773 774 if (dvms(lp1)/=0.0d0) then 775 ekb(lp1)=dvms(lp1)/eta(lp1) 776 renorm(lp1)=1.0d0/sqrt(dvms(lp1)) 777 ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze) 778 ckb(lp1)=eta(lp1)/sqrt(dvms(lp1)) 779 else 780 ekb(lp1)=0.0d0 781 end if 782 end if 783 end do 784 785 ! Loop on angular momenta 786 do lp1=1,lmax+1 787 788 ! Compute form factor if ekb(lp1) not 0 789 if (ekb(lp1)/=0.0d0) then 790 791 ! do q=0 separately, non-zero if l=0 792 if(lp1==1)then 793 do ir=1,mmax 794 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 795 work1(ir)=rad(ir)*dvwf 796 end do 797 call der_int(work1,work2,rad,dr,mmax-1,result) 798 ffspl(1,1,lp1)=result 799 else 800 ! For l non-zero, f(q=0) vanishes ! 801 ffspl(1,1,lp1)=0.0d0 802 end if 803 804 ! Prepare loop over q values 805 irmax=mmax+1 806 do ir=mmax,2,-1 807 test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir) 808 work5(ir)=test 809 work1(ir)=0.0d0 810 ! Will ignore tail within decade of machine precision 811 if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then 812 irmax=ir 813 end if 814 end do 815 ! Increase irmax a bit 816 irmax=irmax+4 817 ! Ask irmax to be lower than mmax 818 if(irmax>mmax-1)irmax=mmax-1 819 820 ABI_ALLOCATE(work3,(irmax-1)) 821 ABI_ALLOCATE(work4,(irmax-1)) 822 823 ! Loop over q values 824 do iq=2,mqgrid 825 tpiq=two_pi*qgrid(iq) 826 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 827 work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility) 828 work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility) 829 830 ! Handle r=0 separately 831 work1(1)=0.0d0 832 call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax)) 833 do ir=2,irmax 834 work1(ir)=besjx(ir)*work5(ir) 835 end do 836 ! do integral 837 call der_int(work1,work2,rad,dr,irmax,result) 838 ffspl(iq,1,lp1)=result 839 end do 840 841 ! Compute yp1=derivative of f(q) at q=0 842 if(lp1/=2)then 843 ! For l/=1, yp1=0 844 yp1=0.0d0 845 else 846 ! For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3 847 do ir=1,irmax 848 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 849 work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0 850 end do 851 call der_int(work1,work2,rad,dr,irmax,result) 852 yp1=result 853 end if 854 855 ! Compute ypn=derivative of f(q) at q=qgrid(mqgrid) 856 tpiq=two_pi*qgrid(mqgrid) 857 ! Treat ir=1, r=0, separately 858 work1(1)=0.0d0 859 ! Here, must distinguish l==0 from others 860 if(lp1==1)then 861 ! l==0 : ypn=$\int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$ 862 ! The sine and cosine of the last point were computed in the previous loop 863 ! So, there is no need to call sincos. Note that the rank of besj is 1. 864 call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax)) 865 do ir=2,irmax 866 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 867 work1(ir)=-besjx(ir)*two_pi*rad(ir)*rad(ir)*dvwf 868 end do 869 else 870 ! l==1 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$ 871 ! l==2 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$ 872 ! l==3 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$ 873 ! The sine and cosine of the last point were computed in the previous loop 874 ! Store first previously computed value with besj of order l, then use 875 ! besj of order l-1 (=lp1-2) 876 work1(2:irmax)=besjx(2:irmax) 877 call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax)) 878 do ir=2,irmax 879 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 880 work1(ir)=(two_pi*rad(ir)**2)*dvwf*& 881 & ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) ) 882 end do 883 end if 884 ! work1 is ready for integration 885 call der_int(work1,work2,rad,dr,irmax,result) 886 ypn=result 887 888 ! Fit spline to get second derivatives by spline fit 889 call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,& 890 & ffspl(:,2,lp1)) 891 892 ABI_DEALLOCATE(work3) 893 ABI_DEALLOCATE(work4) 894 895 else 896 897 ! KB energy is zero, put nonlocal correction at l=0 to 0 898 ffspl(:,:,lp1)=0.0d0 899 900 end if 901 902 ! End loop on angular momenta 903 end do 904 905 ABI_DEALLOCATE(work1) 906 ABI_DEALLOCATE(work2) 907 ABI_DEALLOCATE(work_spl) 908 ABI_DEALLOCATE(work5) 909 ABI_DEALLOCATE(besjx) 910 911 ! End of lmax/=-1 condition 912 end if 913 914 end subroutine psp1nl
ABINIT/sincos [ Functions ]
NAME
sincos
FUNCTION
Update the sine and cosine values, needed inside the pseudopotential routines psp1lo and psp1nl.
INPUTS
iq = number of current wavevector q irmax = number of values of r on the radial grid to be computed mmax = dimension of pspwk and rad pspwk(:,1,1) and pspwk(:,2,1) : sine and cosine of 2$\pi$ dq * rad pspwk(:,1,2) and pspwk(:,2,2) : sine and cosine of 2$\pi$ previous q * rad rad(mmax) radial grid tpiq = 2 $\pi$ * current wavevector q
OUTPUT
pspwk(*,1,2) and pspwk(*,2,2) : sine and cosine of 2$\pi$ current q * rad
NOTES
The speed was a special concern, so iterative computation based on addition formula is possible. Interestingly, this algorithm places strong constraints on accuracy, so this routine is machine-dependent.
PARENTS
psp1lo,psp1nl
CHILDREN
SOURCE
1061 subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq) 1062 1063 use defs_basis 1064 use m_profiling_abi 1065 1066 !This section has been created automatically by the script Abilint (TD). 1067 !Do not modify the following lines by hand. 1068 #undef ABI_FUNC 1069 #define ABI_FUNC 'sincos' 1070 !End of the abilint section 1071 1072 implicit none 1073 1074 !Arguments ------------------------------------ 1075 !scalars 1076 integer,intent(in) :: iq,irmax,mmax 1077 real(dp),intent(in) :: tpiq 1078 !arrays 1079 real(dp),intent(in) :: rad(mmax) 1080 real(dp),intent(inout) :: pspwk(mmax,2,2) 1081 1082 !Local variables------------------------------- 1083 !scalars 1084 integer :: ir,nstep 1085 real(dp) :: prevcos,prevsin 1086 logical :: testmipspro 1087 #if defined HAVE_LINALG_MLIB 1088 real(dp) :: halfpi 1089 #endif 1090 1091 1092 ! ************************************************************************* 1093 1094 #if defined HAVE_LINALG_MLIB 1095 halfpi=asin(1.0d0) 1096 #endif 1097 1098 if(iq==2)then 1099 1100 ! Here set up the sin and cos at iq=2 1101 do ir=2,irmax 1102 pspwk(ir,1,2)=pspwk(ir,1,1) 1103 pspwk(ir,2,2)=pspwk(ir,2,1) 1104 end do 1105 1106 else 1107 ! 1108 ! The sensitivity of the algorithm to changes of nstep 1109 ! has been tested : for all the machines except SGI - R10000 , 1110 ! either using only the hard way, or 1111 ! using up to nstep=40 causes changes at the level 1112 ! of 1.0d-16 in the total energy. Larger values of 1113 ! nstep might be possible, but the associated residual 1114 ! is already very small ! The accelerated computation of 1115 ! sine and cosine is essential for a good speed on IBM, but, 1116 ! fortunately, on the SGI - R10000 the normal computation is fast enough. 1117 1118 testmipspro=.false. 1119 #ifdef FC_MIPSPRO 1120 testmipspro=.true. 1121 #endif 1122 nstep=40 1123 if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then 1124 1125 ! Every nstep steps, uses the hard way 1126 do ir=2,irmax 1127 #if defined HAVE_LINALG_MLIB 1128 ! There is a bug in the hp library !! Sine is slightly inaccurate ! 1129 pspwk(ir,1,2)=cos(tpiq*rad(ir)-halfpi) 1130 #else 1131 pspwk(ir,1,2)=sin(tpiq*rad(ir)) 1132 #endif 1133 pspwk(ir,2,2)=cos(tpiq*rad(ir)) 1134 end do 1135 1136 else 1137 1138 ! Here the fastest way, iteratively 1139 do ir=2,irmax 1140 prevsin=pspwk(ir,1,2) 1141 prevcos=pspwk(ir,2,2) 1142 pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1) 1143 pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1) 1144 end do 1145 1146 end if 1147 1148 end if ! iq==2 1149 1150 end subroutine sincos