TABLE OF CONTENTS
ABINIT/m_psp1 [ Modules ]
NAME
m_psp1
FUNCTION
Initialize pspcod=1 or 4 pseudopotential (Teter format)
COPYRIGHT
Copyright (C) 1998-2022 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_psp1 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_splines 28 29 use m_special_funcs, only : besjm 30 use m_psptk, only : psp1cc 31 32 implicit none 33 34 private
m_psp1/der_int [ Functions ]
[ Top ] [ m_psp1 ] [ 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 $.
SOURCE
853 subroutine der_int(ff,df,rr,dr,nlast,smf) 854 855 !Arguments ------------------------------------ 856 !nmax sets standard number of grid points ! SHOULD BE REMOVED 857 !scalars 858 integer,parameter :: nmax=2000 859 integer,intent(in) :: nlast 860 real(dp),intent(out) :: smf 861 !no_abirules 862 !Note that dimension here starts at 0 863 real(dp), intent(in) :: dr(0:nmax),ff(0:nmax),rr(0:nmax) 864 real(dp), intent(out) :: df(0:nmax) 865 866 !Local variables------------------------------- 867 !scalars 868 integer :: jj 869 real(dp),parameter :: div12=1.d0/12.d0 870 real(dp) :: hh 871 character(len=500) :: message 872 873 ! ************************************************************************* 874 875 !Check that nlast lie within 0 to nmax 876 if (nlast<0.or.nlast>nmax) then 877 write(message, '(a,i12,a,i12)' )& 878 & ' nlast=',nlast,' lies outside range [0,nmax] with dimension nmax=',nmax 879 ABI_BUG(message) 880 end if 881 882 !Compute derivatives at lower end, near r=0 883 df(0)=-25.d0/12.d0*ff(0)+4.d0*ff(1)-3.d0*ff(2)+4.d0/3.d0*ff(3)& 884 & -1.d0/4.d0*ff(4) 885 df(1)=-1.d0/4.d0*ff(0)-5.d0/6.d0*ff(1)+3.d0/2.d0*ff(2)& 886 & -1.d0/2.d0*ff(3)+1.d0/12.d0*ff(4) 887 888 !Run over range from just past r=0 to near r(n), using central differences 889 do jj=2,nlast-2 890 df(jj)=(ff(jj-2)-8.d0*(ff(jj-1)-ff(jj+1))-ff(jj+2))*div12 891 end do 892 893 !Compute derivative at upper end of range 894 if (nlast < 4) then 895 message = ' der_int: ff does not have enough elements. nlast is too low' 896 ABI_ERROR(message) 897 end if 898 899 df(nlast-1)=-1.d0/12.d0*ff(nlast-4)& 900 & +1.d0/2.d0*ff(nlast-3)& 901 & -3.d0/2.d0*ff(nlast-2)& 902 & +5.d0/6.d0*ff(nlast-1)& 903 & +1.d0/4.d0*ff(nlast) 904 df(nlast)=1.d0/4.d0*ff(nlast-4)& 905 & -4.d0/3.d0*ff(nlast-3)& 906 & +3.d0*ff(nlast-2)& 907 & -4.d0*ff(nlast-1)& 908 & +25.d0/12.d0*ff(nlast) 909 910 !Apply correct normalization over full range 911 do jj=0,nlast 912 df(jj)=df(jj)*dr(jj) 913 end do 914 915 smf=0.d0 916 do jj=0,nlast-1 917 hh=rr(jj+1)-rr(jj) 918 smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1))) 919 end do 920 smf=smf/12.d0 921 922 end subroutine der_int
m_psp1/psp1in [ Functions ]
[ Top ] [ m_psp1 ] [ 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.
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 .
SOURCE
106 subroutine psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,& 107 & e990,e999,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 108 & mmax,mpsang,mqgrid,nproj,n1xccc,pspcod,& 109 & qchrg,qgrid,rcpsp,rms,useylm,vlspl,xcccrc,xccc1d,& 110 & zion,znucl) 111 112 !Arguments ------------------------------------ 113 !scalars 114 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc,pspcod 115 integer,intent(in) :: useylm 116 real(dp),intent(in) :: dq,zion,znucl 117 real(dp),intent(out) :: epsatm,qchrg,xcccrc 118 !arrays 119 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang) 120 real(dp),intent(in) :: qgrid(mqgrid) 121 real(dp),intent(out) :: e990(mpsang),e999(mpsang),ekb(lnmax),ekb1(mpsang) 122 real(dp),intent(out) :: ekb2(mpsang),epspsp(mpsang) 123 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) 124 real(dp),intent(out) :: rcpsp(mpsang),rms(mpsang),vlspl(mqgrid,2) 125 real(dp),intent(inout) :: xccc1d(n1xccc,6) 126 127 !Local variables------------------------------- 128 !scalars 129 integer :: ii,iln,index,ipsang,kk,lhigh,ll,mm,nlmax 130 real(dp) :: arg,dq2pi,fchrg,rchrg,xx,yp1,ypn 131 character(len=500) :: message,errmsg 132 !arrays 133 real(dp),allocatable :: drad(:),ekb_tmp(:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:) 134 real(dp),allocatable :: vpspll(:,:),wfll(:,:),wksincos(:,:,:),work_space(:) 135 real(dp),allocatable :: work_spl1(:),work_spl2(:) 136 137 ! *************************************************************************** 138 139 !Note: Teter s grid is hard-coded at mmax=2001 140 !mmax was read from the pseudopotential file in the calling routine 141 if (mmax/=2001) then 142 write(message, '(a,i12,a,a,a,a)' )& 143 & 'Using Teter grid (pspcod=1 and 4) but mmax=',mmax,ch10,& 144 & 'mmax must be 2001 for Teter grid.',ch10,& 145 & 'Action: check your pseudopotential input file.' 146 ABI_ERROR(message) 147 end if 148 149 !File format of formatted Teter psp input (the 3 first lines 150 !have already been read in calling -pspatm- routine) : 151 152 !(1) title (character) line 153 !(2) znucl,zion,pspdat 154 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well 155 !For each angular momentum : 156 !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll) 157 !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll) 158 !(6) rchrg,fchrg,qchrg 159 !(7) ll 160 !(8) (vpsp(j,ll),j=0,nmax) 161 !Then for iproj=1 to 2 162 !for ll=0,lmax 163 !(10) ll 164 !(11) ((upsp(j,ll,iproj),j=0,nmax) 165 166 do ipsang=1,lmax+1 167 168 read (tmp_unit,*,err=10,iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang) 169 write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) & 170 & ipsang-1,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),& 171 & 'l,e99.0,e99.9,nproj,rcpsp' 172 call wrtout(ab_out,message,'COLL') 173 call wrtout(std_out, message,'COLL') 174 175 read (tmp_unit,*,err=10,iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang) 176 write(message, '(4f13.8,t55,a)' ) & 177 & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),& 178 & ' rms, ekb1, ekb2, epsatm' 179 call wrtout(ab_out,message,'COLL') 180 call wrtout(std_out, message,'COLL') 181 182 end do 183 184 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn 185 index=0;iln=0;indlmn(:,:)=0 186 do ipsang=1,lmax+1 187 if(nproj(ipsang)>0)then 188 ll=ipsang-1 189 do kk=1,nproj(ipsang) 190 iln=iln+1 191 do mm=1,2*ll*useylm+1 192 index=index+1 193 indlmn(1,index)=ll 194 indlmn(2,index)=mm-ll*useylm-1 195 indlmn(3,index)=kk 196 indlmn(4,index)=ll*ll+(1-useylm)*ll+mm 197 indlmn(5,index)=iln 198 indlmn(6,index)=1 199 end do 200 end do 201 end if 202 end do 203 204 read (tmp_unit,*,err=10,iomsg=errmsg) rchrg,fchrg,qchrg 205 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 206 call wrtout(ab_out,message,'COLL') 207 call wrtout(std_out, message,'COLL') 208 209 ! Generate core charge function and derivatives, if needed 210 if(fchrg>1.0d-15)then 211 if(pspcod==1)then 212 call psp1cc(fchrg,n1xccc,xccc1d) 213 ! The core charge function for pspcod=1 becomes zero beyond 3*rchrg only. 214 ! Thus xcccrc must be set equal to 3*rchrg . 215 xcccrc=3*rchrg 216 else if(pspcod==4)then 217 call psp4cc(fchrg,n1xccc,xccc1d) 218 ! For pspcod=4, the core charge cut off exactly beyond rchrg 219 xcccrc=rchrg 220 end if 221 else 222 xcccrc=0.0d0 223 xccc1d(:,:)=0.0d0 224 end if 225 226 !-------------------------------------------------------------------- 227 !Will now proceed at the reading of pots and wfs, as well as their treatment 228 229 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials 230 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg. 231 !rad(:)=radial grid r(i) 232 !drad(:)= inverse of d(r(i))/d(i) for radial grid 233 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions 234 235 ABI_MALLOC(vloc,(mmax)) 236 ABI_MALLOC(vpspll,(mmax,mpsang)) 237 if(lmax==-1) vpspll(:,:)=zero 238 239 !(1) Read atomic pseudopotential for each l, filling up array vpspll 240 !Note: put each l into vpspll(:,l+1) 241 do ipsang=1,lmax+1 242 read (tmp_unit,*,err=10,iomsg=errmsg) ll 243 read (tmp_unit,*,err=10,iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax) 244 end do 245 246 !Copy appropriate nonlocal psp for use as local one 247 vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 ) 248 249 !(2) Create radial grid, and associated quantities 250 ABI_MALLOC(rad,(mmax)) 251 ABI_MALLOC(drad,(mmax)) 252 ABI_MALLOC(wksincos,(mmax,2,2)) 253 254 !Teter grid--need both r and dr in this case 255 do ii=0,mmax-1 256 xx=dble(ii)/dble(mmax-1) 257 rad (ii+1)=100.d0*(xx+.01d0)**5-1.d-8 258 drad(ii+1)=500.d0*(xx+.01d0)**4/dble(mmax-1) 259 end do 260 261 !here compute sin(r(:)*dq) and cos(r(:)*dq) 262 !NOTE: also invert dr !! 263 dq2pi=2.0d0*pi*dq 264 do ii=1,mmax 265 arg=dq2pi*rad(ii) 266 drad(ii)=1.0d0/drad(ii) 267 wksincos(ii,1,1)=sin(arg) 268 wksincos(ii,2,1)=cos(arg) 269 end do 270 271 !(3)Carry out calculations for local (lloc) pseudopotential. 272 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q). 273 ABI_MALLOC(work_space,(mqgrid)) 274 ABI_MALLOC(work_spl1,(mqgrid)) 275 ABI_MALLOC(work_spl2,(mqgrid)) 276 call psp1lo(drad,epsatm,mmax,mqgrid,qgrid,& 277 & work_spl1,rad,vloc,wksincos,yp1,ypn,zion) 278 279 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 280 call spline (qgrid,work_spl1,mqgrid,yp1,ypn,work_spl2) 281 vlspl(:,1)=work_spl1(:) 282 vlspl(:,2)=work_spl2(:) 283 284 ABI_FREE(work_space) 285 ABI_FREE(work_spl1) 286 ABI_FREE(work_spl2) 287 288 !(4)Take care of non-local part 289 290 !Zero out all Kleinman-Bylander energies to initialize 291 ekb(:)=0.0d0 292 !write(std_out,*)' psp1in : before nonlocal corrections ' 293 !write(std_out,*)' psp1in : lloc, lmax = ',lloc,lmax 294 295 !Allow for option of no nonlocal corrections (lloc=lmax=0) 296 if (lloc==0.and.lmax==0) then 297 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 298 call wrtout(ab_out,message,'COLL') 299 call wrtout(std_out, message,'COLL') 300 301 else 302 303 ! Proceed to make Kleinman-Bylander form factors for each l up to lmax 304 305 ! Read wavefunctions for each l up to lmax 306 ABI_MALLOC(wfll,(mmax,mpsang)) 307 do ipsang=1,lmax+1 308 ! For pspcod==4, wfs for the local angular momentum are not written 309 if (nproj(ipsang)/=0 .or. pspcod==1) then 310 read (tmp_unit,*,err=10,iomsg=errmsg) ll 311 if (ipsang/=ll+1) then 312 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 313 & 'Pseudopotential input file does not have',ch10,& 314 & 'angular momenta in order expected for first projection',& 315 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 316 & 'Action: check your pseudopotential input file.' 317 ABI_ERROR(message) 318 end if 319 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 320 321 else 322 wfll(:,ipsang)=0.0d0 323 end if 324 325 end do 326 ! ---------------------------------------------------------------------- 327 ! Compute KB form factors and fit splines 328 329 ! nlmax is highest l for which a nonlocal correction is being computed 330 nlmax=lmax 331 if (lloc==lmax) nlmax=lmax-1 332 ! write(std_out,*)' psp1in : lmax,lloc=',lmax,lloc 333 ABI_MALLOC(ekb_tmp,(mpsang,2)) 334 ABI_MALLOC(ffspl_tmp,(mqgrid,2,nlmax+1,2)) 335 336 call psp1nl(drad,ekb_tmp(:,1),ffspl_tmp(:,:,:,1),lloc,& 337 & nlmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos) 338 339 ! Read second wavefunction for second projection operator 340 ! (only read cases where nproj(ll)=2) --also find highest l for which nproj(l)=2 341 lhigh=-1 342 do ipsang=1,min(lmax+1,mpsang) 343 if (nproj(ipsang)==2) then 344 lhigh=ipsang-1 345 read (tmp_unit,*,err=10,iomsg=errmsg) ll 346 if (ipsang/=ll+1) then 347 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 348 & 'Pseudopotential input file does not have',ch10,& 349 & 'angular momenta in order expected for second projection',& 350 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 351 & 'Action: check your pseudopotential input file.' 352 ABI_ERROR(message) 353 end if 354 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 355 356 else 357 wfll(:,ipsang)=0.0d0 358 359 end if 360 end do 361 362 ! Compute KB form factors and fit splines for second wf if any 363 364 if (lhigh>-1) then 365 call psp1nl(drad,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lloc,& 366 & lhigh,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos) 367 end if 368 369 ! Convert ekb and ffspl 370 iln=0 371 do ii=1,lmnmax 372 kk=indlmn(5,ii) 373 if (kk>iln) then 374 iln=kk 375 ekb(kk)=ekb_tmp(1+indlmn(1,ii),indlmn(3,ii)) 376 ! write(std_out,*)' psp1in : lmnmax,ii,indlmn(1,ii)=',lmnmax,ii,indlmn(1,ii) 377 ffspl(:,:,kk)=ffspl_tmp(:,:,1+indlmn(1,ii),indlmn(3,ii)) 378 end if 379 end do 380 381 ABI_FREE(ekb_tmp) 382 ABI_FREE(ffspl_tmp) 383 ABI_FREE(wfll) 384 end if 385 386 ABI_FREE(vpspll) 387 ABI_FREE(rad) 388 ABI_FREE(drad) 389 ABI_FREE(vloc) 390 ABI_FREE(wksincos) 391 392 return 393 394 ! Handle IO error 395 10 continue 396 ABI_ERROR(errmsg) 397 398 end subroutine psp1in
m_psp1/psp1lo [ Functions ]
[ Top ] [ m_psp1 ] [ 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).
SOURCE
431 subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,& 432 & vloc,wksincos,yp1,ypn,zion) 433 434 !Arguments ------------------------------------ 435 !scalars 436 integer,intent(in) :: mmax,mqgrid 437 real(dp),intent(in) :: zion 438 real(dp),intent(out) :: epsatm,yp1,ypn 439 !arrays 440 real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 441 real(dp),intent(inout) :: wksincos(mmax,2,2) 442 real(dp),intent(out) :: q2vq(mqgrid) 443 444 !Local variables------------------------------- 445 !scalars 446 integer,parameter :: mma0=2001 447 integer :: iq,ir,irmax 448 real(dp),parameter :: scale=10.0d0 449 real(dp) :: result,test,tpiq 450 !arrays 451 real(dp) :: wk(mma0),wk1(mma0),wk2(mma0) 452 453 ! ************************************************************************* 454 455 !Do q=0 separately (compute epsatm) 456 !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr] 457 !Treat r=0 by itself 458 wk(1)=0.0d0 459 460 do ir=2,mmax 461 ! (at large r do not want prefactor of r^2 and should see 462 ! V(r)+Zv/r go to 0 at large r) 463 test=vloc(ir)+zion/rad(ir) 464 ! write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test 465 ! In this routine, NO cut-off radius is imposed : the input 466 ! vloc MUST be in real(dp) to obtain numerically 467 ! accurate values. The error can be on the order of 0.001 Ha ! 468 if (abs(test)<1.0d-20) then 469 wk(ir)=0.0d0 470 else 471 wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion) 472 end if 473 end do 474 !Do integral from 0 to r(max) (disregard contrib beyond r(max) 475 !(need numerical derivatives to do integral) 476 !Use mmax-1 to convert to Teter s dimensioning starting at 0 477 call der_int(wk,wk2,rad,drad,mmax-1,result) 478 479 epsatm=4.d0*pi*(result) 480 !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi 481 q2vq(1)=-zion/pi 482 483 !Prepare loop over q values 484 irmax=mmax+1 485 do ir=mmax,2,-1 486 test=vloc(ir)+zion/rad(ir) 487 wk1(ir)=test*rad(ir) 488 ! Will ignore tail within decade of machine precision 489 if ((scale+abs(test))==scale .and. irmax==ir+1) then 490 irmax=ir 491 end if 492 end do 493 !Increase irmax a bit : this is copied from psp1nl 494 irmax=irmax+4 495 if(irmax>mmax)irmax=mmax 496 497 !Loop over q values 498 do iq=2,mqgrid 499 tpiq=two_pi*qgrid(iq) 500 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 501 ! set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral 502 !$\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]$. 503 ! Handle r=0 separately 504 wk(1)=0.0d0 505 do ir=2,irmax 506 wk(ir)=wksincos(ir,1,2)*wk1(ir) 507 end do 508 ! do integral from 0 to r(max) 509 if(irmax>mmax-1)irmax=mmax-1 510 511 call der_int(wk,wk2,rad,drad,irmax,result) 512 ! store q^2 v(q) 513 q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result 514 end do 515 516 !Compute derivatives of q^2 v(q) at ends of interval 517 yp1=0.0d0 518 !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]$ 519 !integral from r(mmax) to infinity is overkill; ignore 520 !set up integrand 521 !Handle r=0 separately 522 wk(1)=0.0d0 523 tpiq=two_pi*qgrid(mqgrid) 524 do ir=2,mmax 525 test=vloc(ir)+zion/rad(ir) 526 ! Ignore contributions within decade of machine precision 527 if ((scale+abs(test))==scale) then 528 wk(ir)=0.0d0 529 else 530 wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * & 531 & (rad(ir)*vloc(ir)+zion) 532 end if 533 end do 534 call der_int(wk,wk2,rad,drad,mmax-1,result) 535 536 ypn=2.0d0*result 537 538 end subroutine psp1lo
m_psp1/psp1nl [ Functions ]
[ Top ] [ m_psp1 ] [ 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.
SOURCE
590 subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,& 591 & qgrid,rad,vloc,vpspll,wfll,wksincos) 592 593 !Arguments ------------------------------------ 594 !scalars 595 integer,intent(in) :: lloc,lmax,mmax,mpsang,mqgrid 596 !arrays 597 real(dp),intent(in) :: dr(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 598 real(dp),intent(in) :: vpspll(mmax,mpsang),wfll(mmax,mpsang) 599 real(dp),intent(inout) :: wksincos(mmax,2,2) 600 real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang) 601 602 !Local variables------------------------------- 603 !scalars 604 integer,parameter :: dpsang=5 605 integer :: iq,ir,irmax,lp1 606 real(dp) :: dvwf,result,test,tpiq,yp1,ypn 607 character(len=500) :: message 608 !arrays 609 real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang) 610 real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:) 611 real(dp),allocatable :: work_spl(:) 612 613 ! ************************************************************************* 614 615 !Zero out Kleinman-Bylander energies ekb 616 ekb(:)=0.0d0 617 !Zero out eta and other parameters too (so 0 s show up in output later) 618 eta(:)=0.0d0 619 dvms(:)=0.0d0 620 ckb(:)=0.0d0 621 622 !Allow for no nonlocal correction (lmax=-1) 623 if (lmax/=-1) then 624 625 ! Check that lmax is within allowed range 626 if (lmax<0.or.lmax>3) then 627 write(message, '(a,i12,a,a,a,a,a,a,a)' )& 628 & 'lmax=',lmax,' is not an allowed value.',ch10,& 629 & 'Allowed values are -1 for no nonlocal correction or else',ch10,& 630 & '0, 1, 2, or 3 for maximum l nonlocal correction.',ch10,& 631 & 'Action: check the input atomic psp data file for lmax.' 632 ABI_ERROR(message) 633 end if 634 635 ! Compute normalizing integrals eta=<dV> and mean square 636 ! nonlocal psp correction dvms=<dV^2> 637 ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction 638 639 ABI_MALLOC(work1,(mmax+1)) 640 ABI_MALLOC(work2,(mmax+1)) 641 ABI_MALLOC(work_spl,(mqgrid)) 642 ABI_MALLOC(work5,(mmax)) 643 ABI_MALLOC(besjx,(mmax)) 644 645 do lp1=1,lmax+1 646 647 ! Only do the work if nonlocal correction is nonzero 648 if (lp1 /= lloc+1) then 649 650 ! integrand for 0 to r(mmax) 651 do ir=1,mmax 652 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 653 work1(ir)=wfll(ir,lp1)*dvwf 654 end do 655 656 ! do integral 657 ! first need derivative of function; note use of 658 ! shifted indices to accomodate Mike Teter s choice of 0:mmax-1 659 call der_int(work1,work2,rad,dr,mmax-1,result) 660 eta(lp1)=result 661 662 ! DEBUG 663 ! write(std_out,*)' psp1nl : write eta(lp1)' 664 ! write(std_out,*)result 665 ! do ir=1,mmax,61 666 ! write(std_out,*)vpspll(ir,lp1),vloc(ir),wfll(ir,lp1) 667 ! end do 668 ! write(std_out,*) 669 ! do ir=1,mmax,61 670 ! write(std_out,*)work1(ir),rad(ir),dr(ir) 671 ! end do 672 ! ENDDEBUG 673 674 do ir=1,mmax 675 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 676 work1(ir)=dvwf**2 677 end do 678 call der_int(work1,work2,rad,dr,mmax-1,result) 679 680 dvms(lp1)=result 681 682 ! If dvms is not 0 for any given angular momentum l, 683 ! compute Xavier Gonze s definition of the Kleinman-Bylander 684 ! energy E_KB = dvms/eta. In this case also renormalize 685 ! the projection operator to u_KB(r)=$u_l(r) dV(r)/\sqrt{dvms}$. 686 ! This means dvwf gets multiplied by the normalization factor 687 ! "renorm"=$1/\sqrt{dvms}$ as seen below. 688 ! With dvwf=dV(r)*wf(r) for wf(r)=``radial'' wf, the integrand 689 ! for each angular momentum l is 690 ! Bessel_l(2 $\pi$ q r) * wf(r) * dV(r) * r; 691 ! NOTE presence of extra r in integrand. 692 693 if (dvms(lp1)/=0.0d0) then 694 ekb(lp1)=dvms(lp1)/eta(lp1) 695 renorm(lp1)=1.0d0/sqrt(dvms(lp1)) 696 ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze) 697 ckb(lp1)=eta(lp1)/sqrt(dvms(lp1)) 698 else 699 ekb(lp1)=0.0d0 700 end if 701 end if 702 end do 703 704 ! Loop on angular momenta 705 do lp1=1,lmax+1 706 707 ! Compute form factor if ekb(lp1) not 0 708 if (ekb(lp1)/=0.0d0) then 709 710 ! do q=0 separately, non-zero if l=0 711 if(lp1==1)then 712 do ir=1,mmax 713 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 714 work1(ir)=rad(ir)*dvwf 715 end do 716 call der_int(work1,work2,rad,dr,mmax-1,result) 717 ffspl(1,1,lp1)=result 718 else 719 ! For l non-zero, f(q=0) vanishes ! 720 ffspl(1,1,lp1)=0.0d0 721 end if 722 723 ! Prepare loop over q values 724 irmax=mmax+1 725 do ir=mmax,2,-1 726 test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir) 727 work5(ir)=test 728 work1(ir)=0.0d0 729 ! Will ignore tail within decade of machine precision 730 if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then 731 irmax=ir 732 end if 733 end do 734 ! Increase irmax a bit 735 irmax=irmax+4 736 ! Ask irmax to be lower than mmax 737 if(irmax>mmax-1)irmax=mmax-1 738 739 ABI_MALLOC(work3,(irmax-1)) 740 ABI_MALLOC(work4,(irmax-1)) 741 742 ! Loop over q values 743 do iq=2,mqgrid 744 tpiq=two_pi*qgrid(iq) 745 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 746 work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility) 747 work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility) 748 749 ! Handle r=0 separately 750 work1(1)=0.0d0 751 call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax)) 752 do ir=2,irmax 753 work1(ir)=besjx(ir)*work5(ir) 754 end do 755 ! do integral 756 call der_int(work1,work2,rad,dr,irmax,result) 757 ffspl(iq,1,lp1)=result 758 end do 759 760 ! Compute yp1=derivative of f(q) at q=0 761 if(lp1/=2)then 762 ! For l/=1, yp1=0 763 yp1=0.0d0 764 else 765 ! For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3 766 do ir=1,irmax 767 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 768 work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0 769 end do 770 call der_int(work1,work2,rad,dr,irmax,result) 771 yp1=result 772 end if 773 774 ! Compute ypn=derivative of f(q) at q=qgrid(mqgrid) 775 tpiq=two_pi*qgrid(mqgrid) 776 ! Treat ir=1, r=0, separately 777 work1(1)=0.0d0 778 ! Here, must distinguish l==0 from others 779 if(lp1==1)then 780 ! l==0 : ypn=$\int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$ 781 ! The sine and cosine of the last point were computed in the previous loop 782 ! So, there is no need to call sincos. Note that the rank of besj is 1. 783 call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax)) 784 do ir=2,irmax 785 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 786 work1(ir)=-besjx(ir)*two_pi*rad(ir)*rad(ir)*dvwf 787 end do 788 else 789 ! l==1 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$ 790 ! l==2 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$ 791 ! l==3 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$ 792 ! The sine and cosine of the last point were computed in the previous loop 793 ! Store first previously computed value with besj of order l, then use 794 ! besj of order l-1 (=lp1-2) 795 work1(2:irmax)=besjx(2:irmax) 796 call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax)) 797 do ir=2,irmax 798 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 799 work1(ir)=(two_pi*rad(ir)**2)*dvwf*& 800 & ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) ) 801 end do 802 end if 803 ! work1 is ready for integration 804 call der_int(work1,work2,rad,dr,irmax,result) 805 ypn=result 806 807 ! Fit spline to get second derivatives by spline fit 808 call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,& 809 & ffspl(:,2,lp1)) 810 811 ABI_FREE(work3) 812 ABI_FREE(work4) 813 814 else 815 ! KB energy is zero, put nonlocal correction at l=0 to 0 816 ffspl(:,:,lp1)=0.0d0 817 end if 818 819 end do ! End loop on angular momenta 820 821 ABI_FREE(work1) 822 ABI_FREE(work2) 823 ABI_FREE(work_spl) 824 ABI_FREE(work5) 825 ABI_FREE(besjx) 826 end if ! End of lmax/=-1 condition 827 828 end subroutine psp1nl
m_psp1/psp4cc [ Functions ]
[ Top ] [ m_psp1 ] [ Functions ]
NAME
psp4cc
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for the format 4 of pseudopotentials. This is a even polynomial of 24th order for core density, that is cut off exactly beyond rchrg. It has been produced on 7 May 1992 by M. Teter.
INPUTS
fchrg=magnitude of the core charge correction n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
NOTES
The argument of xccc1d is assumed to be normalized, and to vary from xx=0 to 1 (from r=0 to r=xcccrc)
WARNINGS
the fifth derivative is not yet delivered.
SOURCE
1047 subroutine psp4cc(fchrg,n1xccc,xccc1d) 1048 1049 !Arguments ------------------------------------ 1050 !scalars 1051 integer,intent(in) :: n1xccc 1052 real(dp),intent(in) :: fchrg 1053 !arrays 1054 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 1055 1056 !Local variables------------------------------- 1057 !scalars 1058 integer :: i1xccc,ider 1059 real(dp),parameter :: a10=-0.1156854803757563d5,a12=+0.2371534625455588d5 1060 real(dp),parameter :: a14=-0.3138755797827918d5,a16=+0.2582842713241039d5 1061 real(dp),parameter :: a18=-0.1200356429115204d5,a20=+0.2405099057118771d4 1062 real(dp),parameter :: a2=-0.8480751097855989d1,a4=+0.9684600878284791d2 1063 real(dp),parameter :: a6=-0.7490894651588015d3,a8=+0.3670890998130434d4 1064 real(dp) :: der1,dern,factor 1065 character(len=500) :: message 1066 !arrays 1067 real(dp),allocatable :: ff(:),ff2(:),work(:),xx(:) 1068 real(dp) :: x 1069 1070 ! ************************************************************************* 1071 1072 ABI_MALLOC(ff,(n1xccc)) 1073 ABI_MALLOC(ff2,(n1xccc)) 1074 ABI_MALLOC(work,(n1xccc)) 1075 ABI_MALLOC(xx,(n1xccc)) 1076 1077 1078 if(n1xccc > 1)then 1079 factor=1.0d0/dble(n1xccc-1) 1080 do i1xccc=1,n1xccc 1081 xx(i1xccc)=(i1xccc-1)*factor 1082 end do 1083 else 1084 write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc 1085 ABI_BUG(message) 1086 end if 1087 1088 !Initialization, to avoid some problem with some compilers 1089 xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero 1090 1091 !Take care of each derivative separately 1092 do ider=0,2 1093 1094 if(ider==0)then 1095 ! Generate spline fitting for the function gg 1096 do i1xccc=1,n1xccc 1097 ! ff(i1xccc)=fchrg*gg(xx(i1xccc)) 1098 ff(i1xccc)=fchrg*gg_psp4(xx(i1xccc)) 1099 end do 1100 ! Complete with derivatives at end points 1101 der1=0.0d0 1102 ! dern=fchrg*gp(1.0d0) 1103 dern=fchrg*gp_psp4(1.0d0) 1104 else if(ider==1)then 1105 ! Generate spline fitting for the function gp 1106 do i1xccc=1,n1xccc 1107 ! ff(i1xccc)=fchrg*gp(xx(i1xccc)) 1108 ff(i1xccc)=fchrg*gp_psp4(xx(i1xccc)) 1109 end do 1110 ! Complete with derivatives at end points, already estimated 1111 der1=xccc1d(1,ider+2) 1112 dern=xccc1d(n1xccc,ider+2) 1113 else if(ider==2)then 1114 ! Generate spline fitting for the function gpp 1115 ! (note : the function gpp has already been estimated, for the spline 1116 ! fitting of the function gg, but it is replaced here by the more 1117 ! accurate analytic derivative) 1118 do i1xccc=1,n1xccc 1119 x=xx(i1xccc) 1120 ff(i1xccc)=fchrg*(gpp_1_psp4(x)+gpp_2_psp4(x)+gpp_3_psp4(x)) 1121 ! ff(i1xccc)=fchrg*gpp(xx(i1xccc)) 1122 end do 1123 ! Complete with derivatives of end points 1124 der1=xccc1d(1,ider+2) 1125 dern=xccc1d(n1xccc,ider+2) 1126 end if 1127 1128 ! Produce second derivative numerically, for use with splines 1129 call spline(xx,ff,n1xccc,der1,dern,ff2) 1130 xccc1d(:,ider+1)=ff(:) 1131 xccc1d(:,ider+3)=ff2(:) 1132 end do 1133 1134 xccc1d(:,6)=zero 1135 1136 ABI_FREE(ff) 1137 ABI_FREE(ff2) 1138 ABI_FREE(work) 1139 ABI_FREE(xx) 1140 1141 !DEBUG 1142 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 1143 !write(std_out,*)' xx gg gp ' 1144 !do i1xccc=1,n1xccc 1145 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 1146 !end do 1147 !write(std_out,*)' xx gpp gg2 ' 1148 !do i1xccc=1,n1xccc 1149 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 1150 !end do 1151 !write(std_out,*)' xx gp2 gpp2 ' 1152 !do i1xccc=1,n1xccc 1153 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 1154 !end do 1155 !write(std_out,*)' psp1cc : debug done, stop ' 1156 !stop 1157 !ENDDEBUG 1158 1159 contains 1160 1161 function gg_psp4(x) 1162 !Expression of 7 May 1992 1163 real(dp) :: gg_psp4 1164 real(dp),intent(in) :: x 1165 gg_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 1166 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 1167 & x**2*(a18+x**2*(a20))))))))))) *(1.0d0-x**2)**2 1168 end function gg_psp4 1169 1170 function gp_psp4(x) 1171 !gp(x) is the derivative of gg(x) wrt x 1172 real(dp) :: gp_psp4 1173 real(dp),intent(in) :: x 1174 gp_psp4=2.d0*x*((a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( & 1175 & 4.d0*a8+x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( & 1176 & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(10.d0*a20))))))))))*& 1177 & (1.d0-x**2)**2 & 1178 & -2.0d0*(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 1179 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 1180 & x**2*(a18+x**2*a20)))))))))) *(1.0d0-x**2) ) 1181 end function gp_psp4 1182 1183 function gpp_1_psp4(x) 1184 !gpp(x) is the second derivative of gg(x) wrt x 1185 real(dp) :: gpp_1_psp4 1186 real(dp),intent(in) :: x 1187 gpp_1_psp4= ( 2.d0*a4+ x**2*(3.d0*2.d0*a6 +x**2*( & 1188 & 4.d0*3.d0*a8+ x**2*(5.d0*4.d0*a10+x**2*( & 1189 & 6.d0*5.d0*a12+x**2*(7.d0*6.d0*a14+x**2*( & 1190 & 8.d0*7.d0*a16+x**2*(9.d0*8.d0*a18+x**2*( & 1191 & 10.d0*9.d0*a20) & 1192 & ))))))))*(2.d0*x*(1.d0-x**2))**2 1193 end function gpp_1_psp4 1194 1195 function gpp_2_psp4(x) 1196 1197 real(dp) :: gpp_2_psp4 1198 real(dp),intent(in) :: x 1199 gpp_2_psp4=(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( & 1200 & 4.d0*a8 +x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( & 1201 & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*( & 1202 & 10.d0*a20) & 1203 & )))))))))*(1.d0-x**2)*2*(1.d0-9.d0*x**2) 1204 end function gpp_2_psp4 1205 1206 function gpp_3_psp4(x) 1207 1208 real(dp) :: gpp_3_psp4 1209 real(dp),intent(in) :: x 1210 gpp_3_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 1211 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 1212 & x**2*(a18+x**2*a20 & 1213 & ))))))))))*(1.0d0-3.d0*x**2)*(-4.d0) 1214 end function gpp_3_psp4 1215 1216 end subroutine psp4cc
m_psp1/sincos [ Functions ]
[ Top ] [ m_psp1 ] [ 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.
SOURCE
953 subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq) 954 955 !Arguments ------------------------------------ 956 !scalars 957 integer,intent(in) :: iq,irmax,mmax 958 real(dp),intent(in) :: tpiq 959 !arrays 960 real(dp),intent(in) :: rad(mmax) 961 real(dp),intent(inout) :: pspwk(mmax,2,2) 962 963 !Local variables------------------------------- 964 !scalars 965 integer :: ir,nstep 966 real(dp) :: prevcos,prevsin 967 logical :: testmipspro 968 969 970 ! ************************************************************************* 971 972 if(iq==2)then 973 974 ! Here set up the sin and cos at iq=2 975 do ir=2,irmax 976 pspwk(ir,1,2)=pspwk(ir,1,1) 977 pspwk(ir,2,2)=pspwk(ir,2,1) 978 end do 979 980 else 981 ! 982 ! The sensitivity of the algorithm to changes of nstep 983 ! has been tested : for all the machines except SGI - R10000 , 984 ! either using only the hard way, or 985 ! using up to nstep=40 causes changes at the level 986 ! of 1.0d-16 in the total energy. Larger values of 987 ! nstep might be possible, but the associated residual 988 ! is already very small ! The accelerated computation of 989 ! sine and cosine is essential for a good speed on IBM, but, 990 ! fortunately, on the SGI - R10000 the normal computation is fast enough. 991 992 testmipspro=.false. 993 nstep=40 994 if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then 995 996 ! Every nstep steps, uses the hard way 997 do ir=2,irmax 998 pspwk(ir,1,2)=sin(tpiq*rad(ir)) 999 pspwk(ir,2,2)=cos(tpiq*rad(ir)) 1000 end do 1001 1002 else 1003 1004 ! Here the fastest way, iteratively 1005 do ir=2,irmax 1006 prevsin=pspwk(ir,1,2) 1007 prevcos=pspwk(ir,2,2) 1008 pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1) 1009 pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1) 1010 end do 1011 1012 end if 1013 1014 end if ! iq==2 1015 1016 end subroutine sincos