TABLE OF CONTENTS


ABINIT/m_psp8 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psp8

FUNCTION

 Initialize pspcod=8 (pseudopotentials in the format generated by DRH):

COPYRIGHT

  Copyright (C) 1999-2018 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 .

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_psp8
28 
29  use defs_basis
30  use m_errors
31  use m_abicore
32  use m_splines
33 
34  use defs_datatypes,  only : nctab_t
35  use m_pawrad,        only : pawrad_type, pawrad_init, pawrad_free
36  use m_psps,          only : nctab_eval_tvalespl
37  use m_psptk,  only : psp8lo, psp8nl
38 
39  implicit none
40 
41  private

ABINIT/psp8in [ Functions ]

[ Top ] [ Functions ]

NAME

 psp8in

FUNCTION

 Initialize pspcod=8 (pseudopotentials in the format generated by DRH):
 continue to read the corresponding file, then compute the
 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= 2*maximum angular momentum for nonlocal pseudopotentials - 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 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

PARENTS

      pspatm

CHILDREN

      nctab_eval_tvalespl,pawrad_free,pawrad_init,psp8cc,psp8lo,psp8nl,spline
      wrtout

SOURCE

109 subroutine psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
110 &                  mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,&
111 &                  useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad)
112 
113 
114 !This section has been created automatically by the script Abilint (TD).
115 !Do not modify the following lines by hand.
116 #undef ABI_FUNC
117 #define ABI_FUNC 'psp8in'
118 !End of the abilint section
119 
120  implicit none
121 
122 !Arguments ------------------------------------
123 !scalars
124  integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid,mqgrid_vl
125  integer,intent(in) :: pspso,n1xccc,useylm
126  real(dp),intent(in) :: zion,znucl
127  real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad
128  type(nctab_t),intent(inout) :: nctab
129 !arrays
130  integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang)
131  real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl)
132  real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
133  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
134 
135 !Local variables-------------------------------
136 !scalars
137  integer :: extension_switch,iln,iln0,pspindex,ipsang,irad,jj,kk,ll,ll_err,llin
138  integer :: mm,nn,nso
139  real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn
140  logical :: has_tvale
141  character(len=500) :: message,errmsg
142  type(pawrad_type) :: mesh
143 !arrays
144  integer, allocatable :: nproj_tmp(:)
145  real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_spl(:)
146 
147 ! ***************************************************************************
148 
149 !File format of formatted drh psp input, as adapted for use
150 !by the ABINIT code (the 3 first lines
151 !have already been read in calling -pspatm- routine) :
152 
153 !(1) title (character) line
154 !(2) znucl,zion,pspdat
155 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well  (r2well not used)
156 !(4) rchrg,fchrg,qchrg  (fchrg /=0 if core charge, qchrg not used)
157 !(5) nproj(0:lmax)  (several projectors allowed for each l)
158 !(6) extension_switch(2) (spin-orbit parameters)
159 !Then, for ll=0,lmax :
160 !if(nproj(ll)>0)
161 !1/<u1|vbkb1>, 1/<u2|vbkb2>, ...
162 !for  irad=1,mmax  : irad, r(irad), vbkb1(irad,ll), vbkb2(irad,ll), ...
163 !elseif ll=lloc
164 !for  irad=1,mmax  : irad, r(irad), vloc(irad)
165 !end if
166 !
167 !If(lloc>lmax)
168 !for  irad=1,mmax  : irad, r(irad), vloc(irad)
169 !end if
170 !
171 !vbkb are Bloechl-Kleinman-Bylander projectors,(vpsp(r,ll)-vloc(r))*u(r,ll),
172 !unnormalized
173 !Note that an arbitrary local potential is allowed.  Set lloc>lmax, and
174 !provide projectors for all ll<=lmax
175 !
176 !Finally, if fchrg>0,
177 !for  irad=1,mmax  : irad, r(irad), xccc(irad),
178 !xccc'(irac), xccc''(irad), xccc'''(irad), xccc''''(irad)
179 !
180 !Model core charge for nonlinear core xc correction, and 4 derivatives
181 
182  read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
183  write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
184  call wrtout(ab_out,message,'COLL')
185  call wrtout(std_out,  message,'COLL')
186 
187  ABI_ALLOCATE(nproj_tmp,(mpssoang))
188 
189  nproj_tmp(:)=0
190  read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(1:lmax+1)
191  write(message, '(a,5i6)' ) '     nproj',nproj_tmp(1:lmax+1)
192  call wrtout(ab_out,message,'COLL')
193  call wrtout(std_out,message,'COLL')
194 
195 !place holder for future implementation of additional optional header
196 !lines without invalidating existing psp files
197 !Now (12/2014) extended to include spin-orbit projectors
198 
199 ! The integer labeled "extension switch" on line 6
200 ! of the *.psp8 file will be set to 1 (non- or scalar-relativistic)
201 ! or 3 (relativistic) to signal this to Abinit that the file contains the pseudo valence charge.
202 
203  has_tvale = .False.
204  read (tmp_unit,*, err=10, iomsg=errmsg) extension_switch
205  if (any(extension_switch==[2, 3])) then
206    read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(lmax+2:2*lmax+1)
207    write(message, '(5x,a,i6)' ) 'spin-orbit psp, extension_switch',extension_switch
208    call wrtout(ab_out,message,'COLL')
209    call wrtout(std_out,  message,'COLL')
210    write(message, '(5x,a,5i6)' ) '   nprojso',nproj_tmp(lmax+2:2*lmax+1)
211    call wrtout(ab_out,message,'COLL')
212    call wrtout(std_out,  message,'COLL')
213    has_tvale =  (extension_switch == 3)
214  else if (any(extension_switch==[0,1])) then
215    write(message, '(5x,a,i6)' ) 'extension_switch',extension_switch
216    call wrtout(ab_out,message,'COLL')
217    call wrtout(std_out,  message,'COLL')
218    has_tvale =  (extension_switch == 1)
219  else
220    write(message, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,&
221 &   'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit'
222    MSG_ERROR(message)
223  end if
224 
225  if(lloc<4) then
226    if (nproj_tmp(lloc+1)/=0) then
227      write(message, '(a,i4,a,a,i4,5a)' )&
228 &     'Pseudopotential input file has nproj=',nproj_tmp(lloc+1),ch10,&
229 &     'for angular momentum',lloc,' which is the local potential.',ch10,&
230 &     'Should be 0 for the local potential',ch10,&
231 &     'Action: check your pseudopotential input file.'
232      MSG_ERROR(message)
233    end if
234  end if
235 
236 !--------------------------------------------------------------------
237 
238 !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn
239 ! if(pspso==2) then
240  if (any(extension_switch == [0,1])) then
241    nso=1
242  else if (any(extension_switch == [2,3])) then
243    nso=2
244    if (pspso==0) then
245      write (message, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,&
246 &     ' but internal variable pspso is 0 - this is not possible for oncvpsp. Set nspinor to 2 and so_psp 1'
247      MSG_COMMENT(message)
248    end if
249  else
250    write(message, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,&
251 &   'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit'
252    MSG_ERROR(message)
253  end if
254 
255  pspindex=0;iln=0;indlmn(:,:)=0
256  do nn=1,nso
257    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
258      ll=ipsang-(nn-1)*lmax-1
259      if (nproj_tmp(ipsang)>0) then
260        do kk=1,nproj_tmp(ipsang)
261          iln=iln+1
262          do mm=1,2*ll*useylm+1
263            pspindex=pspindex+1
264            indlmn(1,pspindex)=ll
265            indlmn(2,pspindex)=mm-ll*useylm-1
266            indlmn(3,pspindex)=kk
267            indlmn(4,pspindex)=ll*ll+(1-useylm)*ll+mm
268            indlmn(5,pspindex)=iln
269            indlmn(6,pspindex)=nn
270          end do
271        end do
272      end if
273    end do
274  end do
275 
276 ! repackage nproj_tmp for proper use by pspatm
277  nproj(:)=0
278  nproj(1:lmax+1)=nproj_tmp(1:lmax+1)
279  if(pspso==2) then
280    nproj(mpsang+1:mpsang+lmax)=nproj_tmp(lmax+2:2*lmax+1)
281  end if
282 
283 !Can now allocate grids, potentials and projectors
284  ABI_ALLOCATE(rad,(mmax))
285  ABI_ALLOCATE(vloc,(mmax))
286  ABI_ALLOCATE(vpspll,(mmax,lnmax))
287 
288 !Will now proceed at the reading of pots and projectors
289 
290 !rad(:)=radial grid r(i)
291 !vpspll(:,1),...,vpspll(:,lnmax)=nonlocal projectors
292 !vloc(:)=local potential
293 
294 !Read Vanderbilt-Kleinman-Bylander energies and projectors for each l
295 !or read local potential for l=lloc.
296 !Also get rad array (actually read more than once)
297  ll_err=0
298  iln0=0
299  do nn=1,nso
300    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
301      ll=ipsang-(nn-1)*lmax-1
302      if (nproj_tmp(ipsang)>0) then
303        read(tmp_unit,*, err=10, iomsg=errmsg) llin,ekb(iln0+1:iln0+nproj_tmp(ipsang))
304        if(llin/=ll) then
305          ll_err=ipsang
306          exit
307        end if
308        do irad=1,mmax
309          read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vpspll(irad,iln0+1:iln0+nproj_tmp(ipsang))
310        end do
311        iln0=iln0+nproj_tmp(ipsang)
312      elseif(ll==lloc .and. nn==1) then
313        read(tmp_unit,*, err=10, iomsg=errmsg) llin
314        if(llin/=ll) then
315          ll_err=ipsang
316          exit
317        end if
318        do irad=1,mmax
319          read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad)
320        end do
321      end if
322    end do !ipsang
323 !Provision for general local potential /= any angular momentum potential
324    if(nn==1 .and. lloc>lmax) then
325      read(tmp_unit,*, err=10, iomsg=errmsg) llin
326      if(llin==lloc) then
327        do irad=1,mmax
328          read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad)
329        end do
330      else
331        ll_err=lloc+1
332        exit
333      end if
334    end if
335  end do !nn
336 
337  if(ll_err>0) then
338    write(message, '(5a,i4,a,i4,a,a)' )&
339 &   'Pseudopotential input file does not have angular momenta in order',ch10,&
340 &   'or has inconsistent general local potential index',ch10,&
341 &   'Expected',ll_err-1,' , got',ll,ch10,&
342 &   'Action: check your pseudopotential input file.'
343    MSG_ERROR(message)
344  end if
345 
346 !Check that rad grid is linear starting at zero
347  amesh=rad(2)-rad(1)
348  damesh=zero
349  do irad=2,mmax-1
350    damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
351  end do
352  if(damesh>tol8 .or. rad(1)/=zero) then
353    write(message, '(5a)' )&
354 &   'Pseudopotential input file requires linear radial mesh',ch10,&
355 &   'starting at zero.',ch10,&
356 &   'Action: check your pseudopotential input file.'
357    MSG_ERROR(message)
358  end if
359 
360 !Get core charge function and derivatives, if needed
361  if(fchrg>1.0d-15)then
362    call psp8cc(mmax,n1xccc,rchrg,xccc1d)
363 !  The core charge function for pspcod=8 becomes zero beyond rchrg.
364 !  Thus xcccrc must be set equal to rchrg.
365    xcccrc=rchrg
366  else
367    xccc1d(:,:) = zero
368    xcccrc = zero
369    fchrg = zero
370    qchrg = zero
371  end if
372 
373  maxrad = rad(mmax)
374 
375 !!   DEBUG
376 !    write(std_out,*)' xcccrc = ', xcccrc, rchrg
377 !    write(std_out,*)
378 !    write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc
379 !    do ii = 1, n1xccc
380 !    write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),&
381 ! &         xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6)
382 !    enddo
383 !    write(std_out,*)
384 !    stop
385 !!   ENDDEBUG
386 
387 
388 !--------------------------------------------------------------------
389 !Carry out calculations for local (lloc) pseudopotential.
390 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q).
391 
392  call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,vlspl(:,1),rad,vloc,yp1,ypn,zion)
393 
394 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
395  ABI_ALLOCATE(work_spl,(mqgrid))
396  call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
397  vlspl(:,2)=work_spl(:)
398  ABI_DEALLOCATE(work_spl)
399 
400 !!  DEBUG
401 ! write(std_out,*)'# Vlocal = '
402 ! write(std_out,*)' amesh  = ', amesh
403 ! write(std_out,*)' epsatm = ', epsatm
404 ! write(std_out,*)' mmax   = ', mmax
405 ! write(std_out,*)' mqgrid = ', mqgrid
406 ! do ir = 1, mqgrid
407 !   write(std_out,*)'   qgrid = ', ir, qgrid(ir)
408 ! enddo
409 ! do ir = 1, mqgrid
410 !   write(std_out,'(a,i5,2f20.12)')'   iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2)
411 ! enddo
412 ! write(std_out,*)
413 ! do ir = 1, mmax
414 !   write(std_out,*)'   rad   = ', rad(ir), vloc(ir)
415 ! enddo
416 ! write(std_out,*)
417 ! write(std_out,*)' yp1    = ', yp1
418 ! write(std_out,*)' ypn    = ', ypn
419 ! write(std_out,*)' zion   = ', zion
420 ! stop
421 !!  ENDDEBUG
422 
423 
424 !--------------------------------------------------------------------
425 !Take care of non-local part
426 
427 !Allow for option of no nonlocal corrections (lloc=lmax=0)
428  if (lloc==0.and.lmax==0) then
429    write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
430    call wrtout(ab_out,message,'COLL')
431    call wrtout(std_out,message,'COLL')
432  else
433 
434 !  Compute Vanderbilt-KB form factors and fit splines
435    call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,mqgrid,qgrid,rad,vpspll)
436 
437  end if
438 
439 !!  DEBUG
440 ! write(std_out,*)'# KB Projectors = '
441 ! write(std_out,*)' amesh  = ', amesh
442 ! do ir = 1, mqgrid
443 !   do il = 1, lnmax
444 !     write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il)
445 !   enddo
446 ! enddo
447 ! do il = 1, lmnmax
448 !   write(std_out,*)' indlmn = ', il, indlmn(:,il)
449 ! enddo
450 ! write(std_out,*)' lmax   = ', lmax
451 ! write(std_out,*)' lmnmax = ', lmnmax
452 ! write(std_out,*)' lnmax  = ', lnmax
453 ! write(std_out,*)' mmax   = ', mmax
454 ! write(std_out,*)' mqgrid = ', mqgrid
455 ! do ir = 1, mqgrid
456 !   write(std_out,*)'   qgrid = ', ir, qgrid(ir)
457 ! enddo
458 ! do il = 1, lnmax
459 !   write(std_out,*)
460 !   write(std_out,*)'# il = ', il
461 !   do ir = 1, mmax
462 !     write(std_out,*)'   rad   = ', rad(ir), vpspll(ir,il)
463 !   enddo
464 ! enddo
465 ! stop
466 !!  ENDDEBUG
467 
468 ! Read pseudo valence charge in real space on the linear mesh
469 ! and transform it to reciprocal space on a regular grid. Use vloc as workspace.
470  vloc(:) = zero
471  if (has_tvale) then
472    do irad=1,mmax
473      read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad)
474      vloc(irad) = vloc(irad) / four_pi
475    end do
476 
477 !! DEBUG
478 !  do irad = 1, mmax
479 !    write(std_out,*)' Valence Charge  = ', rad(irad), vloc(irad)
480 !  enddo
481 !  stop
482 !! ENDDEBUG
483 
484    ! Check that rad grid is linear starting at zero
485    amesh=rad(2)-rad(1)
486    damesh=zero
487    do irad=2,mmax-1
488      damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
489    end do
490    if(damesh>tol8 .or. rad(1)/=zero) then
491      write(message, '(5a)' )&
492 &     'Pseudopotential input file requires linear radial mesh',ch10,&
493 &     'starting at zero.',ch10,&
494 &     'Action: check your pseudopotential input file.'
495      MSG_ERROR(message)
496    end if
497 
498    !  Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
499    call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh)
500    call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl)
501    call pawrad_free(mesh)
502  end if
503 
504  ABI_DEALLOCATE(vpspll)
505  ABI_DEALLOCATE(vloc)
506  ABI_DEALLOCATE(rad)
507  ABI_DEALLOCATE(nproj_tmp)
508 
509  return
510 
511  ! Handle IO error
512  10 continue
513  MSG_ERROR(errmsg)
514 
515 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

PARENTS

      psp8in

CHILDREN

SOURCE

542 subroutine psp8cc(mmax,n1xccc,rchrg,xccc1d)
543 
544 
545 !This section has been created automatically by the script Abilint (TD).
546 !Do not modify the following lines by hand.
547 #undef ABI_FUNC
548 #define ABI_FUNC 'psp8cc'
549 !End of the abilint section
550 
551  implicit none
552 
553 !Arguments ------------------------------------
554 !scalars
555  integer,intent(in) :: mmax,n1xccc
556  real(dp),intent(in) :: rchrg
557 !arrays
558  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
559 
560 !Local variables-------------------------------
561 !scalars
562  integer :: i1xccc,idum,irad,jj
563  real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx
564  character(len=500) :: message,errmsg
565 !arrays
566  real(dp) :: rscale(5)
567  real(dp),allocatable :: ff(:,:),rad(:)
568 
569 !**********************************************************************
570 
571  ABI_ALLOCATE(ff,(mmax,5))
572  ABI_ALLOCATE(rad,(mmax))
573 
574  pi4i=quarter/pi
575 !
576 !Read from pp file the model core charge and its first 4 derivatives
577 !assumed to be on a linear grid starting at zero.
578 !The input functions contain the 4pi factor, and must be rescaled.
579 
580  do irad=1,mmax
581    read(tmp_unit,*, err=10, iomsg=errmsg) idum,rad(irad),(ff(irad,jj),jj=1,5)
582  end do
583 
584 !Check that rad grid is linear starting at zero
585  amesh=rad(2)-rad(1)
586  damesh=zero
587  do irad=2,mmax-1
588    damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
589  end do
590 
591  if(damesh>tol8 .or. rad(1)/=zero) then
592    write(message, '(5a)' )&
593 &   'Pseudopotential input file requires linear radial mesh',ch10,&
594 &   'starting at zero.',ch10,&
595 &   'Action: check your pseudopotential input file.'
596    MSG_ERROR(message)
597  end if
598 
599 !Check that input rchrg is consistent with last grid point
600  if(rchrg>rad(mmax)) then
601    write(message, '(5a)' )&
602 &   'Pseudopotential input file core charge mesh',ch10,&
603 &   'is inconsistent with rchrg in header.',ch10,&
604 &   'Action: check your pseudopotential input file.'
605    MSG_ERROR(message)
606  end if
607 
608 !Factors for unit range scaling
609  do jj = 1, 5
610    rscale(jj)=rchrg**(jj-1)
611  end do
612 
613 !Generate uniform mesh xx in the box cut by rchrg
614 !and interpolate the core charge and derivatives
615 !Cubic polynomial interpolation is used which is consistent
616 !with the original interpolation of these functions from
617 !a log grid to the input linear grid.
618 
619  dri=1.d0/amesh
620  do i1xccc=1,n1xccc
621    xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
622 
623 !  index to find bracketing input mesh points
624    irad = int(dri * xx) + 1
625    irad = max(irad,2)
626    irad = min(irad,mmax-2)
627 !  interpolation coefficients
628    xp = dri * (xx - rad(irad))
629    xpp1 = xp + one
630    xpm1 = xp - one
631    xpm2 = xp - two
632    c1 = -xp * xpm1 * xpm2 * sixth
633    c2 = xpp1 * xpm1 * xpm2 * half
634    c3 = - xp * xpp1 * xpm2 * half
635    c4 = xp * xpp1 * xpm1 * sixth
636 !  Now do the interpolation on all derivatives for this grid point
637 !  Include 1/4pi normalization and unit range scaling
638    do jj=1,5
639      tff =  c1 * ff(irad - 1, jj) &
640 &     + c2 * ff(irad    , jj) &
641 &     + c3 * ff(irad + 1, jj) &
642 &     + c4 * ff(irad + 2, jj)
643      xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff
644    end do
645  end do
646 
647 !5th derivative is apparently not in use, so set to zero
648  xccc1d(:,6)=zero
649 
650  ABI_DEALLOCATE(ff)
651  ABI_DEALLOCATE(rad)
652 
653  return
654 
655  ! Handle IO error
656  10 continue
657  MSG_ERROR(errmsg)
658 
659 end subroutine psp8cc