TABLE OF CONTENTS
ABINIT/m_psp1 [ Modules ]
NAME
m_psp1
FUNCTION
Initialize pspcod=1 or 4 pseudopotential (Teter format)
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 .
PARENTS
CHILDREN
SOURCE
21 #if defined HAVE_CONFIG_H 22 #include "config.h" 23 #endif 24 25 #include "abi_common.h" 26 27 module m_psp1 28 29 use defs_basis 30 use m_errors 31 use m_abicore 32 use m_splines 33 34 use m_special_funcs, only : besjm 35 use m_psptk, only : psp1cc 36 37 implicit none 38 39 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 $.
PARENTS
psp1lo,psp1nl
CHILDREN
SOURCE
908 subroutine der_int(ff,df,rr,dr,nlast,smf) 909 910 911 !This section has been created automatically by the script Abilint (TD). 912 !Do not modify the following lines by hand. 913 #undef ABI_FUNC 914 #define ABI_FUNC 'der_int' 915 !End of the abilint section 916 917 implicit none 918 919 !Arguments ------------------------------------ 920 !nmax sets standard number of grid points ! SHOULD BE REMOVED 921 !scalars 922 integer,parameter :: nmax=2000 923 integer,intent(in) :: nlast 924 real(dp),intent(out) :: smf 925 !no_abirules 926 !Note that dimension here starts at 0 927 real(dp), intent(in) :: dr(0:nmax),ff(0:nmax),rr(0:nmax) 928 real(dp), intent(out) :: df(0:nmax) 929 930 !Local variables------------------------------- 931 !scalars 932 integer :: jj 933 real(dp),parameter :: div12=1.d0/12.d0 934 real(dp) :: hh 935 character(len=500) :: message 936 937 ! ************************************************************************* 938 939 !Check that nlast lie within 0 to nmax 940 if (nlast<0.or.nlast>nmax) then 941 write(message, '(a,i12,a,i12)' )& 942 & ' nlast=',nlast,' lies outside range [0,nmax] with dimension nmax=',nmax 943 MSG_BUG(message) 944 end if 945 946 !Compute derivatives at lower end, near r=0 947 df(0)=-25.d0/12.d0*ff(0)+4.d0*ff(1)-3.d0*ff(2)+4.d0/3.d0*ff(3)& 948 & -1.d0/4.d0*ff(4) 949 df(1)=-1.d0/4.d0*ff(0)-5.d0/6.d0*ff(1)+3.d0/2.d0*ff(2)& 950 & -1.d0/2.d0*ff(3)+1.d0/12.d0*ff(4) 951 952 !Run over range from just past r=0 to near r(n), using central differences 953 do jj=2,nlast-2 954 df(jj)=(ff(jj-2)-8.d0*(ff(jj-1)-ff(jj+1))-ff(jj+2))*div12 955 end do 956 957 !Compute derivative at upper end of range 958 if (nlast < 4) then 959 message = ' der_int: ff does not have enough elements. nlast is too low' 960 MSG_ERROR(message) 961 end if 962 963 df(nlast-1)=-1.d0/12.d0*ff(nlast-4)& 964 & +1.d0/2.d0*ff(nlast-3)& 965 & -3.d0/2.d0*ff(nlast-2)& 966 & +5.d0/6.d0*ff(nlast-1)& 967 & +1.d0/4.d0*ff(nlast) 968 df(nlast)=1.d0/4.d0*ff(nlast-4)& 969 & -4.d0/3.d0*ff(nlast-3)& 970 & +3.d0*ff(nlast-2)& 971 & -4.d0*ff(nlast-1)& 972 & +25.d0/12.d0*ff(nlast) 973 974 !Apply correct normalization over full range 975 do jj=0,nlast 976 df(jj)=df(jj)*dr(jj) 977 end do 978 979 smf=0.d0 980 do jj=0,nlast-1 981 hh=rr(jj+1)-rr(jj) 982 smf=smf+hh*(6.d0*(ff(jj)+ff(jj+1))+hh*(df(jj)-df(jj+1))) 983 end do 984 smf=smf/12.d0 985 986 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 .
PARENTS
pspatm
CHILDREN
psp1cc,psp1lo,psp1nl,psp4cc,spline,wrtout
SOURCE
117 subroutine psp1in(dq,ekb,ekb1,ekb2,epsatm,epspsp,& 118 & e990,e999,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 119 & mmax,mpsang,mqgrid,nproj,n1xccc,pspcod,& 120 & qchrg,qgrid,rcpsp,rms,useylm,vlspl,xcccrc,xccc1d,& 121 & zion,znucl) 122 123 124 !This section has been created automatically by the script Abilint (TD). 125 !Do not modify the following lines by hand. 126 #undef ABI_FUNC 127 #define ABI_FUNC 'psp1in' 128 !End of the abilint section 129 130 implicit none 131 132 !Arguments ------------------------------------ 133 !scalars 134 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc,pspcod 135 integer,intent(in) :: useylm 136 real(dp),intent(in) :: dq,zion,znucl 137 real(dp),intent(out) :: epsatm,qchrg,xcccrc 138 !arrays 139 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang) 140 real(dp),intent(in) :: qgrid(mqgrid) 141 real(dp),intent(out) :: e990(mpsang),e999(mpsang),ekb(lnmax),ekb1(mpsang) 142 real(dp),intent(out) :: ekb2(mpsang),epspsp(mpsang) 143 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) 144 real(dp),intent(out) :: rcpsp(mpsang),rms(mpsang),vlspl(mqgrid,2) 145 real(dp),intent(inout) :: xccc1d(n1xccc,6) 146 147 !Local variables------------------------------- 148 !scalars 149 integer :: ii,iln,index,ipsang,kk,lhigh,ll,mm,nlmax 150 real(dp) :: arg,dq2pi,fchrg,rchrg,xx,yp1,ypn 151 character(len=500) :: message,errmsg 152 !arrays 153 real(dp),allocatable :: drad(:),ekb_tmp(:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:) 154 real(dp),allocatable :: vpspll(:,:),wfll(:,:),wksincos(:,:,:),work_space(:) 155 real(dp),allocatable :: work_spl1(:),work_spl2(:) 156 157 ! *************************************************************************** 158 159 !Note: Teter s grid is hard-coded at mmax=2001 160 !mmax was read from the pseudopotential file in the calling routine 161 if (mmax/=2001) then 162 write(message, '(a,i12,a,a,a,a)' )& 163 & 'Using Teter grid (pspcod=1 and 4) but mmax=',mmax,ch10,& 164 & 'mmax must be 2001 for Teter grid.',ch10,& 165 & 'Action: check your pseudopotential input file.' 166 MSG_ERROR(message) 167 end if 168 169 !File format of formatted Teter psp input (the 3 first lines 170 !have already been read in calling -pspatm- routine) : 171 172 !(1) title (character) line 173 !(2) znucl,zion,pspdat 174 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well 175 !For each angular momentum : 176 !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll) 177 !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll) 178 !(6) rchrg,fchrg,qchrg 179 !(7) ll 180 !(8) (vpsp(j,ll),j=0,nmax) 181 !Then for iproj=1 to 2 182 !for ll=0,lmax 183 !(10) ll 184 !(11) ((upsp(j,ll,iproj),j=0,nmax) 185 186 do ipsang=1,lmax+1 187 188 read (tmp_unit,*,err=10,iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang) 189 write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) & 190 & ipsang-1,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),& 191 & 'l,e99.0,e99.9,nproj,rcpsp' 192 call wrtout(ab_out,message,'COLL') 193 call wrtout(std_out, message,'COLL') 194 195 read (tmp_unit,*,err=10,iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang) 196 write(message, '(4f13.8,t55,a)' ) & 197 & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),& 198 & ' rms, ekb1, ekb2, epsatm' 199 call wrtout(ab_out,message,'COLL') 200 call wrtout(std_out, message,'COLL') 201 202 end do 203 204 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn 205 index=0;iln=0;indlmn(:,:)=0 206 do ipsang=1,lmax+1 207 if(nproj(ipsang)>0)then 208 ll=ipsang-1 209 do kk=1,nproj(ipsang) 210 iln=iln+1 211 do mm=1,2*ll*useylm+1 212 index=index+1 213 indlmn(1,index)=ll 214 indlmn(2,index)=mm-ll*useylm-1 215 indlmn(3,index)=kk 216 indlmn(4,index)=ll*ll+(1-useylm)*ll+mm 217 indlmn(5,index)=iln 218 indlmn(6,index)=1 219 end do 220 end do 221 end if 222 end do 223 224 read (tmp_unit,*,err=10,iomsg=errmsg) rchrg,fchrg,qchrg 225 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 226 call wrtout(ab_out,message,'COLL') 227 call wrtout(std_out, message,'COLL') 228 229 ! Generate core charge function and derivatives, if needed 230 if(fchrg>1.0d-15)then 231 if(pspcod==1)then 232 call psp1cc(fchrg,n1xccc,xccc1d) 233 ! The core charge function for pspcod=1 becomes zero beyond 3*rchrg only. 234 ! Thus xcccrc must be set equal to 3*rchrg . 235 xcccrc=3*rchrg 236 else if(pspcod==4)then 237 call psp4cc(fchrg,n1xccc,xccc1d) 238 ! For pspcod=4, the core charge cut off exactly beyond rchrg 239 xcccrc=rchrg 240 end if 241 else 242 xcccrc=0.0d0 243 xccc1d(:,:)=0.0d0 244 end if 245 246 !-------------------------------------------------------------------- 247 !Will now proceed at the reading of pots and wfs, as well as their treatment 248 249 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials 250 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg. 251 !rad(:)=radial grid r(i) 252 !drad(:)= inverse of d(r(i))/d(i) for radial grid 253 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions 254 255 ABI_ALLOCATE(vloc,(mmax)) 256 ABI_ALLOCATE(vpspll,(mmax,mpsang)) 257 if(lmax==-1) vpspll(:,:)=zero 258 259 !(1) Read atomic pseudopotential for each l, filling up array vpspll 260 !Note: put each l into vpspll(:,l+1) 261 do ipsang=1,lmax+1 262 read (tmp_unit,*,err=10,iomsg=errmsg) ll 263 read (tmp_unit,*,err=10,iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax) 264 end do 265 266 !Copy appropriate nonlocal psp for use as local one 267 vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 ) 268 269 !(2) Create radial grid, and associated quantities 270 ABI_ALLOCATE(rad,(mmax)) 271 ABI_ALLOCATE(drad,(mmax)) 272 ABI_ALLOCATE(wksincos,(mmax,2,2)) 273 274 !Teter grid--need both r and dr in this case 275 do ii=0,mmax-1 276 xx=dble(ii)/dble(mmax-1) 277 rad (ii+1)=100.d0*(xx+.01d0)**5-1.d-8 278 drad(ii+1)=500.d0*(xx+.01d0)**4/dble(mmax-1) 279 end do 280 281 !here compute sin(r(:)*dq) and cos(r(:)*dq) 282 !NOTE: also invert dr !! 283 dq2pi=2.0d0*pi*dq 284 do ii=1,mmax 285 arg=dq2pi*rad(ii) 286 drad(ii)=1.0d0/drad(ii) 287 wksincos(ii,1,1)=sin(arg) 288 wksincos(ii,2,1)=cos(arg) 289 end do 290 291 !(3)Carry out calculations for local (lloc) pseudopotential. 292 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q). 293 ABI_ALLOCATE(work_space,(mqgrid)) 294 ABI_ALLOCATE(work_spl1,(mqgrid)) 295 ABI_ALLOCATE(work_spl2,(mqgrid)) 296 call psp1lo(drad,epsatm,mmax,mqgrid,qgrid,& 297 & work_spl1,rad,vloc,wksincos,yp1,ypn,zion) 298 299 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 300 call spline (qgrid,work_spl1,mqgrid,yp1,ypn,work_spl2) 301 vlspl(:,1)=work_spl1(:) 302 vlspl(:,2)=work_spl2(:) 303 304 ABI_DEALLOCATE(work_space) 305 ABI_DEALLOCATE(work_spl1) 306 ABI_DEALLOCATE(work_spl2) 307 308 !(4)Take care of non-local part 309 310 !Zero out all Kleinman-Bylander energies to initialize 311 ekb(:)=0.0d0 312 !write(std_out,*)' psp1in : before nonlocal corrections ' 313 !write(std_out,*)' psp1in : lloc, lmax = ',lloc,lmax 314 315 !Allow for option of no nonlocal corrections (lloc=lmax=0) 316 if (lloc==0.and.lmax==0) then 317 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 318 call wrtout(ab_out,message,'COLL') 319 call wrtout(std_out, message,'COLL') 320 321 else 322 323 ! Proceed to make Kleinman-Bylander form factors for each l up to lmax 324 325 ! Read wavefunctions for each l up to lmax 326 ABI_ALLOCATE(wfll,(mmax,mpsang)) 327 do ipsang=1,lmax+1 328 ! For pspcod==4, wfs for the local angular momentum are not written 329 if (nproj(ipsang)/=0 .or. pspcod==1) then 330 read (tmp_unit,*,err=10,iomsg=errmsg) ll 331 if (ipsang/=ll+1) then 332 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 333 & 'Pseudopotential input file does not have',ch10,& 334 & 'angular momenta in order expected for first projection',& 335 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 336 & 'Action: check your pseudopotential input file.' 337 MSG_ERROR(message) 338 end if 339 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 340 341 else 342 wfll(:,ipsang)=0.0d0 343 end if 344 345 end do 346 ! ---------------------------------------------------------------------- 347 ! Compute KB form factors and fit splines 348 349 ! nlmax is highest l for which a nonlocal correction is being computed 350 nlmax=lmax 351 if (lloc==lmax) nlmax=lmax-1 352 ! write(std_out,*)' psp1in : lmax,lloc=',lmax,lloc 353 ABI_ALLOCATE(ekb_tmp,(mpsang,2)) 354 ABI_ALLOCATE(ffspl_tmp,(mqgrid,2,nlmax+1,2)) 355 356 call psp1nl(drad,ekb_tmp(:,1),ffspl_tmp(:,:,:,1),lloc,& 357 & nlmax,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos) 358 359 ! Read second wavefunction for second projection operator 360 ! (only read cases where nproj(ll)=2) --also find highest l for which nproj(l)=2 361 lhigh=-1 362 do ipsang=1,min(lmax+1,mpsang) 363 if (nproj(ipsang)==2) then 364 lhigh=ipsang-1 365 read (tmp_unit,*,err=10,iomsg=errmsg) ll 366 if (ipsang/=ll+1) then 367 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 368 & 'Pseudopotential input file does not have',ch10,& 369 & 'angular momenta in order expected for second projection',& 370 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 371 & 'Action: check your pseudopotential input file.' 372 MSG_ERROR(message) 373 end if 374 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 375 376 else 377 wfll(:,ipsang)=0.0d0 378 379 end if 380 end do 381 382 ! Compute KB form factors and fit splines for second wf if any 383 384 if (lhigh>-1) then 385 call psp1nl(drad,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lloc,& 386 & lhigh,mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll,wksincos) 387 end if 388 389 ! Convert ekb and ffspl 390 iln=0 391 do ii=1,lmnmax 392 kk=indlmn(5,ii) 393 if (kk>iln) then 394 iln=kk 395 ekb(kk)=ekb_tmp(1+indlmn(1,ii),indlmn(3,ii)) 396 ! write(std_out,*)' psp1in : lmnmax,ii,indlmn(1,ii)=',lmnmax,ii,indlmn(1,ii) 397 ffspl(:,:,kk)=ffspl_tmp(:,:,1+indlmn(1,ii),indlmn(3,ii)) 398 end if 399 end do 400 401 ABI_DEALLOCATE(ekb_tmp) 402 ABI_DEALLOCATE(ffspl_tmp) 403 ABI_DEALLOCATE(wfll) 404 end if 405 406 ABI_DEALLOCATE(vpspll) 407 ABI_DEALLOCATE(rad) 408 ABI_DEALLOCATE(drad) 409 ABI_DEALLOCATE(vloc) 410 ABI_DEALLOCATE(wksincos) 411 412 return 413 414 ! Handle IO error 415 10 continue 416 MSG_ERROR(errmsg) 417 418 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).
PARENTS
psp1in
CHILDREN
der_int,sincos
SOURCE
457 subroutine psp1lo(drad,epsatm,mmax,mqgrid,qgrid,q2vq,rad,& 458 & vloc,wksincos,yp1,ypn,zion) 459 460 461 !This section has been created automatically by the script Abilint (TD). 462 !Do not modify the following lines by hand. 463 #undef ABI_FUNC 464 #define ABI_FUNC 'psp1lo' 465 !End of the abilint section 466 467 implicit none 468 469 !Arguments ------------------------------------ 470 !scalars 471 integer,intent(in) :: mmax,mqgrid 472 real(dp),intent(in) :: zion 473 real(dp),intent(out) :: epsatm,yp1,ypn 474 !arrays 475 real(dp),intent(in) :: drad(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 476 real(dp),intent(inout) :: wksincos(mmax,2,2) 477 real(dp),intent(out) :: q2vq(mqgrid) 478 479 !Local variables------------------------------- 480 !scalars 481 integer,parameter :: mma0=2001 482 integer :: iq,ir,irmax 483 real(dp),parameter :: scale=10.0d0 484 real(dp) :: result,test,tpiq 485 !arrays 486 real(dp) :: wk(mma0),wk1(mma0),wk2(mma0) 487 488 ! ************************************************************************* 489 490 !Do q=0 separately (compute epsatm) 491 !Set up integrand for q=0: Int[r^2 (V(r)+Zv/r) dr] 492 !Treat r=0 by itself 493 wk(1)=0.0d0 494 495 do ir=2,mmax 496 ! (at large r do not want prefactor of r^2 and should see 497 ! V(r)+Zv/r go to 0 at large r) 498 test=vloc(ir)+zion/rad(ir) 499 ! write(std_out,'(i4,3es20.10)' )ir,rad(ir),test,rad(ir)*test 500 ! In this routine, NO cut-off radius is imposed : the input 501 ! vloc MUST be in real(dp) to obtain numerically 502 ! accurate values. The error can be on the order of 0.001 Ha ! 503 if (abs(test)<1.0d-20) then 504 wk(ir)=0.0d0 505 else 506 wk(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion) 507 end if 508 end do 509 !Do integral from 0 to r(max) (disregard contrib beyond r(max) 510 !(need numerical derivatives to do integral) 511 !Use mmax-1 to convert to Teter s dimensioning starting at 0 512 call der_int(wk,wk2,rad,drad,mmax-1,result) 513 514 epsatm=4.d0*pi*(result) 515 !q=0 value of integral is -zion/Pi + q^2 * epsatm = -zion/Pi 516 q2vq(1)=-zion/pi 517 518 !Prepare loop over q values 519 irmax=mmax+1 520 do ir=mmax,2,-1 521 test=vloc(ir)+zion/rad(ir) 522 wk1(ir)=test*rad(ir) 523 ! Will ignore tail within decade of machine precision 524 if ((scale+abs(test))==scale .and. irmax==ir+1) then 525 irmax=ir 526 end if 527 end do 528 !Increase irmax a bit : this is copied from psp1nl 529 irmax=irmax+4 530 if(irmax>mmax)irmax=mmax 531 532 !Loop over q values 533 do iq=2,mqgrid 534 tpiq=two_pi*qgrid(iq) 535 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 536 ! set up integrand Sin(2Pi q r)(rV(r)+Zv) for integral 537 !$\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]$. 538 ! Handle r=0 separately 539 wk(1)=0.0d0 540 do ir=2,irmax 541 wk(ir)=wksincos(ir,1,2)*wk1(ir) 542 end do 543 ! do integral from 0 to r(max) 544 if(irmax>mmax-1)irmax=mmax-1 545 546 call der_int(wk,wk2,rad,drad,irmax,result) 547 ! store q^2 v(q) 548 q2vq(iq)=-zion/pi+2.d0*qgrid(iq)*result 549 end do 550 551 !Compute derivatives of q^2 v(q) at ends of interval 552 yp1=0.0d0 553 !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]$ 554 !integral from r(mmax) to infinity is overkill; ignore 555 !set up integrand 556 !Handle r=0 separately 557 wk(1)=0.0d0 558 tpiq=two_pi*qgrid(mqgrid) 559 do ir=2,mmax 560 test=vloc(ir)+zion/rad(ir) 561 ! Ignore contributions within decade of machine precision 562 if ((scale+abs(test))==scale) then 563 wk(ir)=0.0d0 564 else 565 wk(ir)=(sin(tpiq*rad(ir))+tpiq*rad(ir)*cos(tpiq*rad(ir))) * & 566 & (rad(ir)*vloc(ir)+zion) 567 end if 568 end do 569 call der_int(wk,wk2,rad,drad,mmax-1,result) 570 571 ypn=2.0d0*result 572 573 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.
PARENTS
psp1in
CHILDREN
besjm,der_int,sincos,spline
SOURCE
631 subroutine psp1nl(dr,ekb,ffspl,lloc,lmax,mmax,mpsang,mqgrid,& 632 & qgrid,rad,vloc,vpspll,wfll,wksincos) 633 634 635 !This section has been created automatically by the script Abilint (TD). 636 !Do not modify the following lines by hand. 637 #undef ABI_FUNC 638 #define ABI_FUNC 'psp1nl' 639 !End of the abilint section 640 641 implicit none 642 643 !Arguments ------------------------------------ 644 !scalars 645 integer,intent(in) :: lloc,lmax,mmax,mpsang,mqgrid 646 !arrays 647 real(dp),intent(in) :: dr(mmax),qgrid(mqgrid),rad(mmax),vloc(mmax) 648 real(dp),intent(in) :: vpspll(mmax,mpsang),wfll(mmax,mpsang) 649 real(dp),intent(inout) :: wksincos(mmax,2,2) 650 real(dp),intent(out) :: ekb(mpsang),ffspl(mqgrid,2,mpsang) 651 652 !Local variables------------------------------- 653 !scalars 654 integer,parameter :: dpsang=5 655 integer :: iq,ir,irmax,lp1 656 real(dp) :: dvwf,result,test,tpiq,yp1,ypn 657 character(len=500) :: message 658 !arrays 659 real(dp) :: ckb(dpsang),dvms(dpsang),eta(dpsang),renorm(dpsang) 660 real(dp),allocatable :: besjx(:),work1(:),work2(:),work3(:),work4(:),work5(:) 661 real(dp),allocatable :: work_spl(:) 662 663 ! ************************************************************************* 664 665 !Zero out Kleinman-Bylander energies ekb 666 ekb(:)=0.0d0 667 !Zero out eta and other parameters too (so 0 s show up in output later) 668 eta(:)=0.0d0 669 dvms(:)=0.0d0 670 ckb(:)=0.0d0 671 672 !Allow for no nonlocal correction (lmax=-1) 673 if (lmax/=-1) then 674 675 ! Check that lmax is within allowed range 676 if (lmax<0.or.lmax>3) then 677 write(message, '(a,i12,a,a,a,a,a,a,a)' )& 678 & 'lmax=',lmax,' is not an allowed value.',ch10,& 679 & 'Allowed values are -1 for no nonlocal correction or else',ch10,& 680 & '0, 1, 2, or 3 for maximum l nonlocal correction.',ch10,& 681 & 'Action: check the input atomic psp data file for lmax.' 682 MSG_ERROR(message) 683 end if 684 685 ! Compute normalizing integrals eta=<dV> and mean square 686 ! nonlocal psp correction dvms=<dV^2> 687 ! "dvwf" consistently refers to dV(r)*wf(r) where dV=nonlocal correction 688 689 ABI_ALLOCATE(work1,(mmax+1)) 690 ABI_ALLOCATE(work2,(mmax+1)) 691 ABI_ALLOCATE(work_spl,(mqgrid)) 692 ABI_ALLOCATE(work5,(mmax)) 693 ABI_ALLOCATE(besjx,(mmax)) 694 695 do lp1=1,lmax+1 696 697 ! Only do the work if nonlocal correction is nonzero 698 if (lp1 /= lloc+1) then 699 700 ! integrand for 0 to r(mmax) 701 do ir=1,mmax 702 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 703 work1(ir)=wfll(ir,lp1)*dvwf 704 end do 705 706 ! do integral 707 ! first need derivative of function; note use of 708 ! shifted indices to accomodate Mike Teter s choice of 0:mmax-1 709 call der_int(work1,work2,rad,dr,mmax-1,result) 710 eta(lp1)=result 711 712 ! DEBUG 713 ! write(std_out,*)' psp1nl : write eta(lp1)' 714 ! write(std_out,*)result 715 ! do ir=1,mmax,61 716 ! write(std_out,*)vpspll(ir,lp1),vloc(ir),wfll(ir,lp1) 717 ! end do 718 ! write(std_out,*) 719 ! do ir=1,mmax,61 720 ! write(std_out,*)work1(ir),rad(ir),dr(ir) 721 ! end do 722 ! ENDDEBUG 723 724 do ir=1,mmax 725 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1) 726 work1(ir)=dvwf**2 727 end do 728 call der_int(work1,work2,rad,dr,mmax-1,result) 729 730 dvms(lp1)=result 731 732 ! If dvms is not 0 for any given angular momentum l, 733 ! compute Xavier Gonze s definition of the Kleinman-Bylander 734 ! energy E_KB = dvms/eta. In this case also renormalize 735 ! the projection operator to u_KB(r)=$u_l(r) dV(r)/\sqrt{dvms}$. 736 ! This means dvwf gets multiplied by the normalization factor 737 ! "renorm"=$1/\sqrt{dvms}$ as seen below. 738 ! With dvwf=dV(r)*wf(r) for wf(r)=``radial'' wf, the integrand 739 ! for each angular momentum l is 740 ! Bessel_l(2 $\pi$ q r) * wf(r) * dV(r) * r; 741 ! NOTE presence of extra r in integrand. 742 743 if (dvms(lp1)/=0.0d0) then 744 ekb(lp1)=dvms(lp1)/eta(lp1) 745 renorm(lp1)=1.0d0/sqrt(dvms(lp1)) 746 ! ckb is Kleinman-Bylander "cosine" (Xavier Gonze) 747 ckb(lp1)=eta(lp1)/sqrt(dvms(lp1)) 748 else 749 ekb(lp1)=0.0d0 750 end if 751 end if 752 end do 753 754 ! Loop on angular momenta 755 do lp1=1,lmax+1 756 757 ! Compute form factor if ekb(lp1) not 0 758 if (ekb(lp1)/=0.0d0) then 759 760 ! do q=0 separately, non-zero if l=0 761 if(lp1==1)then 762 do ir=1,mmax 763 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 764 work1(ir)=rad(ir)*dvwf 765 end do 766 call der_int(work1,work2,rad,dr,mmax-1,result) 767 ffspl(1,1,lp1)=result 768 else 769 ! For l non-zero, f(q=0) vanishes ! 770 ffspl(1,1,lp1)=0.0d0 771 end if 772 773 ! Prepare loop over q values 774 irmax=mmax+1 775 do ir=mmax,2,-1 776 test=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1)*rad(ir) 777 work5(ir)=test 778 work1(ir)=0.0d0 779 ! Will ignore tail within decade of machine precision 780 if ((10.0d0+abs(test))==10.0d0 .and. irmax==ir+1) then 781 irmax=ir 782 end if 783 end do 784 ! Increase irmax a bit 785 irmax=irmax+4 786 ! Ask irmax to be lower than mmax 787 if(irmax>mmax-1)irmax=mmax-1 788 789 ABI_ALLOCATE(work3,(irmax-1)) 790 ABI_ALLOCATE(work4,(irmax-1)) 791 792 ! Loop over q values 793 do iq=2,mqgrid 794 tpiq=two_pi*qgrid(iq) 795 call sincos(iq,irmax,mmax,wksincos,rad,tpiq) 796 work3(:)=wksincos(2:irmax,2,2) !Temporary array (Intel compiler compatibility) 797 work4(:)=wksincos(2:irmax,1,2) !Temporary array (Intel compiler compatibility) 798 799 ! Handle r=0 separately 800 work1(1)=0.0d0 801 call besjm(tpiq,besjx(2:irmax),work3,(lp1-1),irmax-1,work4,rad(2:irmax)) 802 do ir=2,irmax 803 work1(ir)=besjx(ir)*work5(ir) 804 end do 805 ! do integral 806 call der_int(work1,work2,rad,dr,irmax,result) 807 ffspl(iq,1,lp1)=result 808 end do 809 810 ! Compute yp1=derivative of f(q) at q=0 811 if(lp1/=2)then 812 ! For l/=1, yp1=0 813 yp1=0.0d0 814 else 815 ! For l=1, yp1=Int [2 Pi r^2 wf(r) dV(r)]/3 816 do ir=1,irmax 817 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 818 work1(ir)=(two_pi*rad(ir)**2)*dvwf/3.0d0 819 end do 820 call der_int(work1,work2,rad,dr,irmax,result) 821 yp1=result 822 end if 823 824 ! Compute ypn=derivative of f(q) at q=qgrid(mqgrid) 825 tpiq=two_pi*qgrid(mqgrid) 826 ! Treat ir=1, r=0, separately 827 work1(1)=0.0d0 828 ! Here, must distinguish l==0 from others 829 if(lp1==1)then 830 ! l==0 : ypn=$\int [2\pi r (-bes1(2\pi r q)) wf(r) dV(r) r dr]$ 831 ! The sine and cosine of the last point were computed in the previous loop 832 ! So, there is no need to call sincos. Note that the rank of besj is 1. 833 call besjm(tpiq,besjx(2:irmax),work3,1,irmax-1,work4,rad(2:irmax)) 834 do ir=2,irmax 835 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 836 work1(ir)=-besjx(ir)*two_pi*rad(ir)*rad(ir)*dvwf 837 end do 838 else 839 ! l==1 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_0(x)-(2/x)j_1(x)) dr]$ 840 ! l==2 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_1(x)-(3/x)j_2(x)) dr]$ 841 ! l==3 : ypn=$\int [2\pi r^2 wf(r) dV(r) (j_2(x)-(4/x)j_3(x)) dr]$ 842 ! The sine and cosine of the last point were computed in the previous loop 843 ! Store first previously computed value with besj of order l, then use 844 ! besj of order l-1 (=lp1-2) 845 work1(2:irmax)=besjx(2:irmax) 846 call besjm(tpiq,besjx(2:irmax),work3,(lp1-2),irmax-1,work4,rad(2:irmax)) 847 do ir=2,irmax 848 dvwf=(vpspll(ir,lp1)-vloc(ir))*wfll(ir,lp1)*renorm(lp1) 849 work1(ir)=(two_pi*rad(ir)**2)*dvwf*& 850 & ( besjx(ir) - ( dble(lp1)*work1(ir)/(tpiq*rad(ir)) ) ) 851 end do 852 end if 853 ! work1 is ready for integration 854 call der_int(work1,work2,rad,dr,irmax,result) 855 ypn=result 856 857 ! Fit spline to get second derivatives by spline fit 858 call spline(qgrid,ffspl(:,1,lp1),mqgrid,yp1,ypn,& 859 & ffspl(:,2,lp1)) 860 861 ABI_DEALLOCATE(work3) 862 ABI_DEALLOCATE(work4) 863 864 else 865 ! KB energy is zero, put nonlocal correction at l=0 to 0 866 ffspl(:,:,lp1)=0.0d0 867 end if 868 869 end do ! End loop on angular momenta 870 871 ABI_DEALLOCATE(work1) 872 ABI_DEALLOCATE(work2) 873 ABI_DEALLOCATE(work_spl) 874 ABI_DEALLOCATE(work5) 875 ABI_DEALLOCATE(besjx) 876 end if ! End of lmax/=-1 condition 877 878 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.
PARENTS
psp1in
CHILDREN
spline
SOURCE
1146 subroutine psp4cc(fchrg,n1xccc,xccc1d) 1147 1148 1149 !This section has been created automatically by the script Abilint (TD). 1150 !Do not modify the following lines by hand. 1151 #undef ABI_FUNC 1152 #define ABI_FUNC 'psp4cc' 1153 !End of the abilint section 1154 1155 implicit none 1156 1157 !Arguments ------------------------------------ 1158 !scalars 1159 integer,intent(in) :: n1xccc 1160 real(dp),intent(in) :: fchrg 1161 !arrays 1162 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 1163 1164 !Local variables------------------------------- 1165 !scalars 1166 integer :: i1xccc,ider 1167 real(dp),parameter :: a10=-0.1156854803757563d5,a12=+0.2371534625455588d5 1168 real(dp),parameter :: a14=-0.3138755797827918d5,a16=+0.2582842713241039d5 1169 real(dp),parameter :: a18=-0.1200356429115204d5,a20=+0.2405099057118771d4 1170 real(dp),parameter :: a2=-0.8480751097855989d1,a4=+0.9684600878284791d2 1171 real(dp),parameter :: a6=-0.7490894651588015d3,a8=+0.3670890998130434d4 1172 real(dp) :: der1,dern,factor 1173 character(len=500) :: message 1174 !arrays 1175 real(dp),allocatable :: ff(:),ff2(:),work(:),xx(:) 1176 real(dp) :: x 1177 1178 ! ************************************************************************* 1179 1180 ABI_ALLOCATE(ff,(n1xccc)) 1181 ABI_ALLOCATE(ff2,(n1xccc)) 1182 ABI_ALLOCATE(work,(n1xccc)) 1183 ABI_ALLOCATE(xx,(n1xccc)) 1184 1185 1186 if(n1xccc > 1)then 1187 factor=1.0d0/dble(n1xccc-1) 1188 do i1xccc=1,n1xccc 1189 xx(i1xccc)=(i1xccc-1)*factor 1190 end do 1191 else 1192 write(message, '(a,i0)' )' n1xccc should larger than 1, while it is n1xccc=',n1xccc 1193 MSG_BUG(message) 1194 end if 1195 1196 !Initialization, to avoid some problem with some compilers 1197 xccc1d(1,:)=zero ; xccc1d(n1xccc,:)=zero 1198 1199 !Take care of each derivative separately 1200 do ider=0,2 1201 1202 if(ider==0)then 1203 ! Generate spline fitting for the function gg 1204 do i1xccc=1,n1xccc 1205 ! ff(i1xccc)=fchrg*gg(xx(i1xccc)) 1206 ff(i1xccc)=fchrg*gg_psp4(xx(i1xccc)) 1207 end do 1208 ! Complete with derivatives at end points 1209 der1=0.0d0 1210 ! dern=fchrg*gp(1.0d0) 1211 dern=fchrg*gp_psp4(1.0d0) 1212 else if(ider==1)then 1213 ! Generate spline fitting for the function gp 1214 do i1xccc=1,n1xccc 1215 ! ff(i1xccc)=fchrg*gp(xx(i1xccc)) 1216 ff(i1xccc)=fchrg*gp_psp4(xx(i1xccc)) 1217 end do 1218 ! Complete with derivatives at end points, already estimated 1219 der1=xccc1d(1,ider+2) 1220 dern=xccc1d(n1xccc,ider+2) 1221 else if(ider==2)then 1222 ! Generate spline fitting for the function gpp 1223 ! (note : the function gpp has already been estimated, for the spline 1224 ! fitting of the function gg, but it is replaced here by the more 1225 ! accurate analytic derivative) 1226 do i1xccc=1,n1xccc 1227 x=xx(i1xccc) 1228 ff(i1xccc)=fchrg*(gpp_1_psp4(x)+gpp_2_psp4(x)+gpp_3_psp4(x)) 1229 ! ff(i1xccc)=fchrg*gpp(xx(i1xccc)) 1230 end do 1231 ! Complete with derivatives of end points 1232 der1=xccc1d(1,ider+2) 1233 dern=xccc1d(n1xccc,ider+2) 1234 end if 1235 1236 ! Produce second derivative numerically, for use with splines 1237 call spline(xx,ff,n1xccc,der1,dern,ff2) 1238 xccc1d(:,ider+1)=ff(:) 1239 xccc1d(:,ider+3)=ff2(:) 1240 end do 1241 1242 xccc1d(:,6)=zero 1243 1244 ABI_DEALLOCATE(ff) 1245 ABI_DEALLOCATE(ff2) 1246 ABI_DEALLOCATE(work) 1247 ABI_DEALLOCATE(xx) 1248 1249 !DEBUG 1250 !write(std_out,*)' psp1cc : output of core charge density and derivatives ' 1251 !write(std_out,*)' xx gg gp ' 1252 !do i1xccc=1,n1xccc 1253 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2) 1254 !end do 1255 !write(std_out,*)' xx gpp gg2 ' 1256 !do i1xccc=1,n1xccc 1257 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4) 1258 !end do 1259 !write(std_out,*)' xx gp2 gpp2 ' 1260 !do i1xccc=1,n1xccc 1261 !write(std_out,'(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6) 1262 !end do 1263 !write(std_out,*)' psp1cc : debug done, stop ' 1264 !stop 1265 !ENDDEBUG 1266 1267 contains 1268 1269 function gg_psp4(x) 1270 !Expression of 7 May 1992 1271 1272 !This section has been created automatically by the script Abilint (TD). 1273 !Do not modify the following lines by hand. 1274 #undef ABI_FUNC 1275 #define ABI_FUNC 'gg_psp4' 1276 !End of the abilint section 1277 1278 real(dp) :: gg_psp4 1279 real(dp),intent(in) :: x 1280 gg_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 1281 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 1282 & x**2*(a18+x**2*(a20))))))))))) *(1.0d0-x**2)**2 1283 end function gg_psp4 1284 1285 function gp_psp4(x) 1286 !gp(x) is the derivative of gg(x) wrt x 1287 1288 !This section has been created automatically by the script Abilint (TD). 1289 !Do not modify the following lines by hand. 1290 #undef ABI_FUNC 1291 #define ABI_FUNC 'gp_psp4' 1292 !End of the abilint section 1293 1294 real(dp) :: gp_psp4 1295 real(dp),intent(in) :: x 1296 gp_psp4=2.d0*x*((a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( & 1297 & 4.d0*a8+x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( & 1298 & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*(10.d0*a20))))))))))*& 1299 & (1.d0-x**2)**2 & 1300 & -2.0d0*(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 1301 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 1302 & x**2*(a18+x**2*a20)))))))))) *(1.0d0-x**2) ) 1303 end function gp_psp4 1304 1305 function gpp_1_psp4(x) 1306 !gpp(x) is the second derivative of gg(x) wrt x 1307 1308 !This section has been created automatically by the script Abilint (TD). 1309 !Do not modify the following lines by hand. 1310 #undef ABI_FUNC 1311 #define ABI_FUNC 'gpp_1_psp4' 1312 !End of the abilint section 1313 1314 real(dp) :: gpp_1_psp4 1315 real(dp),intent(in) :: x 1316 gpp_1_psp4= ( 2.d0*a4+ x**2*(3.d0*2.d0*a6 +x**2*( & 1317 & 4.d0*3.d0*a8+ x**2*(5.d0*4.d0*a10+x**2*( & 1318 & 6.d0*5.d0*a12+x**2*(7.d0*6.d0*a14+x**2*( & 1319 & 8.d0*7.d0*a16+x**2*(9.d0*8.d0*a18+x**2*( & 1320 & 10.d0*9.d0*a20) & 1321 & ))))))))*(2.d0*x*(1.d0-x**2))**2 1322 end function gpp_1_psp4 1323 1324 function gpp_2_psp4(x) 1325 1326 1327 !This section has been created automatically by the script Abilint (TD). 1328 !Do not modify the following lines by hand. 1329 #undef ABI_FUNC 1330 #define ABI_FUNC 'gpp_2_psp4' 1331 !End of the abilint section 1332 1333 real(dp) :: gpp_2_psp4 1334 real(dp),intent(in) :: x 1335 gpp_2_psp4=(a2+x**2*(2.d0*a4+x**2*(3.d0*a6+x**2*( & 1336 & 4.d0*a8 +x**2*(5.d0*a10+x**2*(6.d0*a12+x**2*( & 1337 & 7.d0*a14+x**2*(8.d0*a16+x**2*(9.d0*a18+x**2*( & 1338 & 10.d0*a20) & 1339 & )))))))))*(1.d0-x**2)*2*(1.d0-9.d0*x**2) 1340 end function gpp_2_psp4 1341 1342 function gpp_3_psp4(x) 1343 1344 1345 !This section has been created automatically by the script Abilint (TD). 1346 !Do not modify the following lines by hand. 1347 #undef ABI_FUNC 1348 #define ABI_FUNC 'gpp_3_psp4' 1349 !End of the abilint section 1350 1351 real(dp) :: gpp_3_psp4 1352 real(dp),intent(in) :: x 1353 gpp_3_psp4=(1.d0+x**2*(a2 +x**2*(a4 +x**2*(a6 +x**2*(a8 + & 1354 & x**2*(a10+x**2*(a12+x**2*(a14+x**2*(a16+ & 1355 & x**2*(a18+x**2*a20 & 1356 & ))))))))))*(1.0d0-3.d0*x**2)*(-4.d0) 1357 end function gpp_3_psp4 1358 1359 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.
PARENTS
psp1lo,psp1nl
CHILDREN
SOURCE
1022 subroutine sincos(iq,irmax,mmax,pspwk,rad,tpiq) 1023 1024 1025 !This section has been created automatically by the script Abilint (TD). 1026 !Do not modify the following lines by hand. 1027 #undef ABI_FUNC 1028 #define ABI_FUNC 'sincos' 1029 !End of the abilint section 1030 1031 implicit none 1032 1033 !Arguments ------------------------------------ 1034 !scalars 1035 integer,intent(in) :: iq,irmax,mmax 1036 real(dp),intent(in) :: tpiq 1037 !arrays 1038 real(dp),intent(in) :: rad(mmax) 1039 real(dp),intent(inout) :: pspwk(mmax,2,2) 1040 1041 !Local variables------------------------------- 1042 !scalars 1043 integer :: ir,nstep 1044 real(dp) :: prevcos,prevsin 1045 logical :: testmipspro 1046 #if defined HAVE_LINALG_MLIB 1047 real(dp) :: halfpi 1048 #endif 1049 1050 1051 ! ************************************************************************* 1052 1053 #if defined HAVE_LINALG_MLIB 1054 halfpi=asin(1.0d0) 1055 #endif 1056 1057 if(iq==2)then 1058 1059 ! Here set up the sin and cos at iq=2 1060 do ir=2,irmax 1061 pspwk(ir,1,2)=pspwk(ir,1,1) 1062 pspwk(ir,2,2)=pspwk(ir,2,1) 1063 end do 1064 1065 else 1066 ! 1067 ! The sensitivity of the algorithm to changes of nstep 1068 ! has been tested : for all the machines except SGI - R10000 , 1069 ! either using only the hard way, or 1070 ! using up to nstep=40 causes changes at the level 1071 ! of 1.0d-16 in the total energy. Larger values of 1072 ! nstep might be possible, but the associated residual 1073 ! is already very small ! The accelerated computation of 1074 ! sine and cosine is essential for a good speed on IBM, but, 1075 ! fortunately, on the SGI - R10000 the normal computation is fast enough. 1076 1077 testmipspro=.false. 1078 #ifdef FC_MIPSPRO 1079 testmipspro=.true. 1080 #endif 1081 nstep=40 1082 if(iq-(iq/nstep)*nstep == 0 .or. testmipspro)then 1083 1084 ! Every nstep steps, uses the hard way 1085 do ir=2,irmax 1086 #if defined HAVE_LINALG_MLIB 1087 ! There is a bug in the hp library !! Sine is slightly inaccurate ! 1088 pspwk(ir,1,2)=cos(tpiq*rad(ir)-halfpi) 1089 #else 1090 pspwk(ir,1,2)=sin(tpiq*rad(ir)) 1091 #endif 1092 pspwk(ir,2,2)=cos(tpiq*rad(ir)) 1093 end do 1094 1095 else 1096 1097 ! Here the fastest way, iteratively 1098 do ir=2,irmax 1099 prevsin=pspwk(ir,1,2) 1100 prevcos=pspwk(ir,2,2) 1101 pspwk(ir,1,2)=prevsin*pspwk(ir,2,1)+prevcos*pspwk(ir,1,1) 1102 pspwk(ir,2,2)=prevcos*pspwk(ir,2,1)-prevsin*pspwk(ir,1,1) 1103 end do 1104 1105 end if 1106 1107 end if ! iq==2 1108 1109 end subroutine sincos