TABLE OF CONTENTS
ABINIT/m_psp8 [ Modules ]
NAME
m_psp8
FUNCTION
Initialize pspcod=8 (pseudopotentials in the format generated by DRH):
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (DRH, XG, AF) 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_psp8 23 24 use defs_basis 25 use m_errors 26 use m_abicore 27 use m_splines 28 29 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 30 use defs_datatypes, only : nctab_t 31 use m_psps, only : nctab_eval_tvalespl 32 use m_psptk, only : psp8lo, psp8nl 33 34 implicit none 35 36 private
ABINIT/psp8in [ Functions ]
NAME
psp8in
FUNCTION
Initialize pspcod=8 (pseudopotentials in the format generated by DRH): continue to read the corresponding file, then compute local and non-local potentials.
INPUTS
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
angular momentum of nonlocal pseudopotential
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mpssoang= Maximum number of channels, including those for treating the spin-orbit coupling
when mpspso=1, mpssoang=mpsang
when mpspso=2, mpssoang=2*mpsang-1
mqgrid=dimension of q (or G) grid for arrays.
mqgrid_vl=dimension of q (or G) grid for valence charge (array qgrid_vl)
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
qgrid_vl(psps%mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for valence charge
pspso=spin-orbit characteristics, govern the content of ffspl and ekb
if =0: this input requires no spin-orbit characteristics of the psp
if =2: this input requires hgh or psp8/upf2 characteristics of the psp
if =3: this input requires hfn characteristics of the psp
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=nuclear number of atom as specified in psp file
OUTPUT
ekb(lnmax)=Kleinman-Bylander energy, as read from input file epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree) 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(mpssoang)=number of projection functions for each angular momentum qchrg is not used, and could be suppressed later vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit xcccrc=XC core correction cutoff radius (bohr) xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file nctab<nctab_t>=NC tables %has_tvale=True if the pseudo contains the pseudo valence charge %tvalespl(mqgrid_vl,2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
SOURCE
98 subroutine psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 99 mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,& 100 useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad) 101 102 !Arguments ------------------------------------ 103 !scalars 104 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid,mqgrid_vl 105 integer,intent(in) :: pspso,n1xccc,useylm 106 real(dp),intent(in) :: zion,znucl 107 real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad 108 type(nctab_t),intent(inout) :: nctab 109 !arrays 110 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang) 111 real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl) 112 real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2) 113 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 114 115 !Local variables------------------------------- 116 !scalars 117 integer :: extension_switch,iln,iln0,pspindex,ipsang,irad,jj,kk,ll,ll_err,llin 118 integer :: mm,nn,nso,ir 119 real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn 120 logical :: has_tvale, debug 121 character(len=500) :: msg,errmsg 122 type(pawrad_type) :: mesh 123 !arrays 124 integer, allocatable :: nproj_tmp(:) 125 real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_spl(:) 126 127 ! *************************************************************************** 128 129 !File format of formatted drh psp input, as adapted for use 130 !by the ABINIT code (the 3 first lines have already been read in calling -pspatm- routine): 131 132 !(1) title (character) line 133 !(2) znucl,zion,pspdat 134 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well (r2well not used) 135 !(4) rchrg,fchrg,qchrg (fchrg /=0 if core charge, qchrg not used) 136 !(5) nproj(0:lmax) (several projectors allowed for each l) 137 !(6) extension_switch(2) (spin-orbit parameters) 138 !Then, for ll=0,lmax : 139 !if(nproj(ll)>0) 140 !1/<u1|vbkb1>, 1/<u2|vbkb2>, ... 141 !for irad=1,mmax : irad, r(irad), vbkb1(irad,ll), vbkb2(irad,ll), ... 142 !else if ll=lloc 143 !for irad=1,mmax : irad, r(irad), vloc(irad) 144 !end if 145 ! 146 !If(lloc>lmax): 147 ! for irad=1,mmax : irad, r(irad), vloc(irad) 148 !end if 149 ! 150 !vbkb are Bloechl-Kleinman-Bylander projectors,(vpsp(r,ll)-vloc(r))*u(r,ll), unnormalized. 151 !Note that an arbitrary local potential is allowed. 152 !Set lloc>lmax, and provide projectors for all ll<=lmax 153 ! 154 !Finally, if fchrg>0: 155 ! 156 ! for irad=1,mmax : irad, r(irad), xccc(irad), 157 ! xccc'(irac), xccc''(irad), xccc'''(irad), xccc''''(irad) 158 ! 159 160 debug = .False.!; debug = .True. 161 162 ! Model core charge for nonlinear core xc correction, and 4 derivatives 163 164 read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg 165 write(msg, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 166 call wrtout([std_out, ab_out], msg) 167 168 ABI_MALLOC(nproj_tmp, (mpssoang)) 169 nproj_tmp = 0 170 171 read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(1:lmax+1) 172 write(msg, '(a,5i6)' ) ' nproj',nproj_tmp(1:lmax+1) 173 call wrtout([std_out, ab_out], msg) 174 175 !place holder for future implementation of additional optional header 176 !lines without invalidating existing psp files 177 !Now (12/2014) extended to include spin-orbit projectors 178 179 ! The integer labeled "extension switch" on line 6 180 ! of the *.psp8 file will be set to 1 (non- or scalar-relativistic) 181 ! or 3 (relativistic) to signal to Abinit that the file contains the pseudo valence charge. 182 183 has_tvale = .False. 184 read (tmp_unit,*, err=10, iomsg=errmsg) extension_switch 185 if (any(extension_switch==[2, 3])) then 186 read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(lmax+2:2*lmax+1) 187 write(msg, '(5x,a,i6)' ) 'spin-orbit psp, extension_switch',extension_switch 188 call wrtout([std_out, ab_out], msg) 189 write(msg, '(5x,a,5i6)' ) ' nprojso',nproj_tmp(lmax+2:2*lmax+1) 190 call wrtout([std_out, ab_out], msg) 191 has_tvale = (extension_switch == 3) 192 else if (any(extension_switch==[0,1])) then 193 write(msg, '(5x,a,i6)' ) 'extension_switch',extension_switch 194 call wrtout([std_out, ab_out], msg) 195 has_tvale = (extension_switch == 1) 196 else 197 write(msg, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,& 198 'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit' 199 ABI_ERROR(msg) 200 end if 201 202 if(lloc<4) then 203 if (nproj_tmp(lloc+1)/=0) then 204 write(msg, '(a,i4,a,a,i4,5a)' )& 205 'Pseudopotential input file has nproj=',nproj_tmp(lloc+1),ch10,& 206 'for angular momentum',lloc,' which is the local potential.',ch10,& 207 'Should be 0 for the local potential',ch10,& 208 'Action: check your pseudopotential input file.' 209 ABI_ERROR(msg) 210 end if 211 end if 212 213 !-------------------------------------------------------------------- 214 215 !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn 216 ! if(pspso==2) then 217 if (any(extension_switch == [0,1])) then 218 nso=1 219 else if (any(extension_switch == [2,3])) then 220 nso=2 221 if (pspso==0) then 222 write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,& 223 ' but internal variable pspso is 0' 224 ABI_COMMENT(msg) 225 end if 226 else 227 write(msg, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,& 228 'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit' 229 ABI_ERROR(msg) 230 end if 231 232 pspindex=0; iln=0; indlmn=0 233 do nn=1,nso 234 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 235 ll = ipsang-(nn-1)*lmax-1 236 if (nproj_tmp(ipsang)>0) then 237 do kk=1,nproj_tmp(ipsang) 238 iln = iln+1 239 do mm=1,2*ll*useylm+1 240 pspindex = pspindex + 1 241 indlmn(1,pspindex) = ll 242 indlmn(2,pspindex) = mm-ll*useylm-1 243 indlmn(3,pspindex) = kk 244 indlmn(4,pspindex) = ll*ll+(1-useylm)*ll+mm 245 indlmn(5,pspindex) = iln 246 indlmn(6,pspindex) = nn 247 !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex 248 end do 249 end do 250 end if 251 end do 252 end do 253 254 ! repackage nproj_tmp for proper use by pspatm 255 nproj(:)=0 256 nproj(1:lmax+1)=nproj_tmp(1:lmax+1) 257 if(pspso==2) then 258 nproj(mpsang+1:mpsang+lmax)=nproj_tmp(lmax+2:2*lmax+1) 259 end if 260 261 !Can now allocate grids, potentials and projectors 262 ABI_MALLOC(rad,(mmax)) 263 ABI_MALLOC(vloc,(mmax)) 264 ABI_MALLOC(vpspll,(mmax,lnmax)) 265 266 !Will now proceed at the reading of pots and projectors 267 268 !rad(:)=radial grid r(i) 269 !vpspll(:,1),...,vpspll(:,lnmax)=nonlocal projectors 270 !vloc(:)=local potential 271 272 !Read Vanderbilt-Kleinman-Bylander energies and projectors for each l 273 !or read local potential for l=lloc. 274 !Also get rad array (actually read more than once) 275 ll_err=0 276 iln0=0 277 do nn=1,nso 278 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 279 ll=ipsang-(nn-1)*lmax-1 280 if (nproj_tmp(ipsang)>0) then 281 read(tmp_unit,*, err=10, iomsg=errmsg) llin,ekb(iln0+1:iln0+nproj_tmp(ipsang)) 282 if(llin/=ll) then 283 ll_err=ipsang 284 exit 285 end if 286 do irad=1,mmax 287 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vpspll(irad,iln0+1:iln0+nproj_tmp(ipsang)) 288 end do 289 iln0=iln0+nproj_tmp(ipsang) 290 elseif(ll==lloc .and. nn==1) then 291 read(tmp_unit,*, err=10, iomsg=errmsg) llin 292 if(llin/=ll) then 293 ll_err=ipsang 294 exit 295 end if 296 do irad=1,mmax 297 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad) 298 end do 299 end if 300 end do !ipsang 301 302 ! Provision for general local potential /= any angular momentum potential 303 if(nn==1 .and. lloc>lmax) then 304 read(tmp_unit,*, err=10, iomsg=errmsg) llin 305 if(llin==lloc) then 306 do irad=1,mmax 307 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad) 308 end do 309 else 310 ll_err=lloc+1 311 exit 312 end if 313 end if 314 end do !nn 315 316 if(ll_err>0) then 317 write(msg, '(5a,i4,a,i4,a,a)' )& 318 'Pseudopotential input file does not have angular momenta in order',ch10,& 319 'or has inconsistent general local potential index',ch10,& 320 'Expected',ll_err-1,' , got',ll,ch10,& 321 'Action: check your pseudopotential input file.' 322 ABI_ERROR(msg) 323 end if 324 325 ! Check that rad grid is linear starting at zero 326 amesh=rad(2)-rad(1) 327 damesh=zero 328 do irad=2,mmax-1 329 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 330 end do 331 if(damesh>tol8 .or. rad(1)/=zero) then 332 write(msg, '(5a)' )& 333 'Pseudopotential input file requires linear radial mesh',ch10,& 334 'starting at zero.',ch10,& 335 'Action: check your pseudopotential input file.' 336 ABI_ERROR(msg) 337 end if 338 339 !Get core charge function and derivatives, if needed 340 if(fchrg>1.0d-15)then 341 call psp8cc(mmax, n1xccc, rchrg, xccc1d) 342 ! The core charge function for pspcod=8 becomes zero beyond rchrg. 343 ! Thus xcccrc must be set equal to rchrg. 344 xcccrc=rchrg 345 else 346 xccc1d(:,:) = zero 347 xcccrc = zero 348 fchrg = zero 349 qchrg = zero 350 end if 351 352 maxrad = rad(mmax) 353 354 !! DEBUG 355 !write(std_out,*)' xcccrc = ', xcccrc, rchrg 356 !write(std_out,*) 357 !write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc 358 !do ii = 1, n1xccc 359 !write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),& 360 ! xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6) 361 !enddo 362 !write(std_out,*) 363 !stop 364 !! ENDDEBUG 365 366 367 !-------------------------------------------------------------------- 368 !Carry out calculations for local (lloc) pseudopotential. 369 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q). 370 371 call psp8lo(amesh, epsatm, mmax, mqgrid, qgrid, vlspl(:,1), rad, vloc, yp1, ypn, zion) 372 373 ! Fit spline to q^2 V(q) (Numerical Recipes subroutine) 374 ABI_MALLOC(work_spl,(mqgrid)) 375 call spline(qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 376 vlspl(:,2)=work_spl(:) 377 ABI_FREE(work_spl) 378 379 if (debug) then 380 write(std_out,*)'# Vlocal psp8 = ' 381 write(std_out,*)' amesh = ', amesh 382 write(std_out,*)' epsatm = ', epsatm 383 write(std_out,*)' mmax = ', mmax 384 write(std_out,*)' mqgrid = ', mqgrid 385 do ir = 1, mqgrid 386 write(std_out,*)' qgrid = ', ir, qgrid(ir) 387 enddo 388 do ir = 1, mqgrid 389 write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2) 390 enddo 391 write(std_out,*) 392 do ir = 1, mmax 393 write(std_out,*)' rad = ', rad(ir), vloc(ir) 394 enddo 395 write(std_out,*) 396 write(std_out,*)' yp1 = ', yp1 397 write(std_out,*)' ypn = ', ypn 398 write(std_out,*)' zion = ', zion 399 stop 400 end if 401 402 403 !-------------------------------------------------------------------- 404 !Take care of non-local part 405 406 !Allow for option of no nonlocal corrections (lloc=lmax=0) 407 if (lloc == 0 .and. lmax == 0) then 408 write(msg, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 409 call wrtout([std_out, ab_out], msg) 410 else 411 412 ! Compute Vanderbilt-KB form factors and fit splines 413 call psp8nl(amesh, ffspl, indlmn, lmax, lmnmax, lnmax, mmax, mqgrid, qgrid, rad, vpspll) 414 end if 415 416 !! DEBUG 417 ! write(std_out,*)'# KB Projectors = ' 418 ! write(std_out,*)' amesh = ', amesh 419 ! do ir = 1, mqgrid 420 ! do il = 1, lnmax 421 ! write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il) 422 ! enddo 423 ! enddo 424 ! do il = 1, lmnmax 425 ! write(std_out,*)' indlmn = ', il, indlmn(:,il) 426 ! enddo 427 ! write(std_out,*)' lmax = ', lmax 428 ! write(std_out,*)' lmnmax = ', lmnmax 429 ! write(std_out,*)' lnmax = ', lnmax 430 ! write(std_out,*)' mmax = ', mmax 431 ! write(std_out,*)' mqgrid = ', mqgrid 432 ! do ir = 1, mqgrid 433 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 434 ! enddo 435 ! do il = 1, lnmax 436 ! write(std_out,*) 437 ! write(std_out,*)'# il = ', il 438 ! do ir = 1, mmax 439 ! write(std_out,*)' rad = ', rad(ir), vpspll(ir,il) 440 ! enddo 441 ! enddo 442 ! stop 443 !! ENDDEBUG 444 445 ! Read pseudo valence charge in real space on the linear mesh 446 ! and transform it to reciprocal space on a regular grid. Use vloc as workspace. 447 vloc(:) = zero 448 if (has_tvale) then 449 do irad=1,mmax 450 read(tmp_unit,*, err=10, iomsg=errmsg)jj, rad(irad), vloc(irad) 451 vloc(irad) = vloc(irad) / four_pi 452 end do 453 454 ! Check that rad grid is linear starting at zero 455 amesh = rad(2) - rad(1); damesh = zero 456 do irad=2,mmax-1 457 damesh = max(damesh, abs(rad(irad)+amesh-rad(irad+1))) 458 end do 459 460 if (damesh > tol8 .or. abs(rad(1)) > tol16) then 461 write(msg,'(3a)')& 462 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,& 463 'Action: check your pseudopotential file.' 464 ABI_ERROR(msg) 465 end if 466 467 ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space. 468 call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh) 469 call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl) 470 call pawrad_free(mesh) 471 end if 472 473 ABI_FREE(vpspll) 474 ABI_FREE(vloc) 475 ABI_FREE(rad) 476 ABI_FREE(nproj_tmp) 477 478 return 479 480 ! Handle IO error 481 10 continue 482 ABI_ERROR(errmsg) 483 484 end subroutine psp8in
m_psp8/psp8cc [ Functions ]
[ Top ] [ m_psp8 ] [ Functions ]
NAME
psp8cc
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for format 8 of the pseudopotentials.
INPUTS
mmax=maximum number of points in real space grid in the psp file n1xccc=dimension of xccc1d; 0 if no XC core correction is used rchrg=cut-off radius for the core density
OUTPUT
xccc1d(n1xccc,6)= 1D core charge function and its four first derivatives
SOURCE
506 subroutine psp8cc(mmax, n1xccc, rchrg, xccc1d) 507 508 !Arguments ------------------------------------ 509 !scalars 510 integer,intent(in) :: mmax,n1xccc 511 real(dp),intent(in) :: rchrg 512 !arrays 513 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 514 515 !Local variables------------------------------- 516 !scalars 517 integer :: i1xccc,idum,irad,jj 518 real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx 519 character(len=500) :: msg,errmsg 520 !arrays 521 real(dp) :: rscale(5) 522 real(dp),allocatable :: ff(:,:),rad(:) 523 524 !********************************************************************** 525 526 ABI_MALLOC(ff,(mmax,5)) 527 ABI_MALLOC(rad,(mmax)) 528 529 pi4i=quarter/pi 530 ! 531 ! Read from pp file the model core charge and its first 4 derivatives 532 ! assumed to be on a linear grid starting at zero. 533 ! The input functions contain the 4pi factor, and must be rescaled. 534 535 do irad=1,mmax 536 read(tmp_unit,*, err=10, iomsg=errmsg) idum,rad(irad),(ff(irad,jj),jj=1,5) 537 end do 538 539 ! Check that rad grid is linear starting at zero 540 amesh=rad(2)-rad(1) 541 damesh=zero 542 do irad=2,mmax-1 543 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 544 end do 545 546 if(damesh>tol8 .or. rad(1)/=zero) then 547 write(msg, '(5a)' )& 548 'Pseudopotential input file requires linear radial mesh',ch10,& 549 'starting at zero.',ch10,& 550 'Action: check your pseudopotential input file.' 551 ABI_ERROR(msg) 552 end if 553 554 ! Check that input rchrg is consistent with last grid point 555 if(rchrg>rad(mmax)) then 556 write(msg, '(5a)' )& 557 'Pseudopotential input file core charge mesh',ch10,& 558 'is inconsistent with rchrg in header.',ch10,& 559 'Action: check your pseudopotential input file.' 560 ABI_ERROR(msg) 561 end if 562 563 !Factors for unit range scaling 564 do jj = 1, 5 565 rscale(jj)=rchrg**(jj-1) 566 end do 567 568 !Generate uniform mesh xx in the box cut by rchrg 569 !and interpolate the core charge and derivatives 570 !Cubic polynomial interpolation is used which is consistent 571 !with the original interpolation of these functions from 572 !a log grid to the input linear grid. 573 574 dri=1.d0/amesh 575 do i1xccc=1,n1xccc 576 xx=(i1xccc-1)* rchrg/dble(n1xccc-1) 577 578 ! index to find bracketing input mesh points 579 irad = int(dri * xx) + 1 580 irad = max(irad,2) 581 irad = min(irad,mmax-2) 582 ! interpolation coefficients 583 xp = dri * (xx - rad(irad)) 584 xpp1 = xp + one 585 xpm1 = xp - one 586 xpm2 = xp - two 587 c1 = -xp * xpm1 * xpm2 * sixth 588 c2 = xpp1 * xpm1 * xpm2 * half 589 c3 = - xp * xpp1 * xpm2 * half 590 c4 = xp * xpp1 * xpm1 * sixth 591 ! Now do the interpolation on all derivatives for this grid point 592 ! Include 1/4pi normalization and unit range scaling 593 do jj=1,5 594 tff = c1 * ff(irad - 1, jj) & 595 & + c2 * ff(irad , jj) & 596 & + c3 * ff(irad + 1, jj) & 597 & + c4 * ff(irad + 2, jj) 598 xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff 599 end do 600 end do 601 602 !5th derivative is apparently not in use, so set to zero 603 xccc1d(:,6)=zero 604 605 ABI_FREE(ff) 606 ABI_FREE(rad) 607 608 return 609 610 ! Handle IO error 611 10 continue 612 ABI_ERROR(errmsg) 613 614 end subroutine psp8cc