!{\src2tex{textfont=tt}}
!!****f* ABINIT/psp8in
!! 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.
!!
!! COPYRIGHT
!! Copyright (C) 1999-2014 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 .
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
!!
!! 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.
!!  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
!!  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
!!
!! PARENTS
!!      pspatm_abinit
!!
!! CHILDREN
!!      psp8cc,psp8lo,psp8nl,spline,wrtout
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

#include "abi_common.h"

subroutine psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
&                  mmax,mpsang,mpssoang,mqgrid,nproj,n1xccc,pspso,qchrg,qgrid,&
&                  useylm,vlspl,xcccrc,xccc1d,zion,znucl)

 use defs_basis
 use m_splines
 use m_errors
 use m_profiling

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'psp8in'
 use interfaces_14_hidewrite
 use interfaces_65_psp, except_this_one => psp8in
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid
 integer,intent(in)  :: pspso,n1xccc,useylm
 real(dp),intent(in) :: zion,znucl
 real(dp),intent(out) :: epsatm,qchrg,xcccrc
!arrays
 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang)
 real(dp),intent(in) :: qgrid(mqgrid)
 real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i

!Local variables-------------------------------
!scalars
 integer :: extension_switch,iln,iln0,index,ipsang,irad,jj,kk,ll,ll_err,llin
 integer :: mm,nn,nso,iln0_so
 real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn
 character(len=500) :: message
!arrays
 integer, allocatable :: nproj_tmp(:)
 real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_space(:),work_spl(:)
 real(dp),allocatable :: ffspl_so(:,:,:,:),ffspl_sr(:,:,:,:)

! ***************************************************************************

!DEBUG
!write(std_out,*)' psp8in : enter and stop '
!stop
!ENDDEBUG


!File format of formatted drh psp input, as adapted for use
!by the ABINIT code (the 3 first lines
!have already been read in calling -pspatm- routine) :

!(1) title (character) line
!(2) znucl,zion,pspdat
!(3) pspcod,pspxc,lmax,lloc,mmax,r2well  (r2well not used)
!(4) rchrg,fchrg,qchrg  (fchrg /=0 if core charge, qchrg not used)
!(5) nproj(0:lmax)  (several projectors allowed for each l)
!Then, for ll=0,lmax :
!if(nproj(ll)>0)
!1/<u1|vbkb1>, 1/<u2|vbkb2>, ...
!for  irad=1,mmax  : irad, r(irad), vbkb1(irad,ll), vbkb2(irad,ll), ...
!elseif ll=lloc
!for  irad=1,mmax  : irad, r(irad), vloc(irad)
!end if
!
!If(lloc>lmax)
!for  irad=1,mmax  : irad, r(irad), vloc(irad)
!end if
!
!vbkb are Bloechl-Kleinman-Bylander projectors,(vpsp(r,ll)-vloc(r))*u(r,ll),
!unnormalized
!Note that an arbitrary local potential is allowed.  Set lloc>lmax, and
!provide projectors for all ll<=lmax
!
!Finally, if fchrg>0,
!for  irad=1,mmax  : irad, r(irad), xccc(irad),
!xccc'(irac), xccc''(irad), xccc'''(irad), xccc''''(irad)
!
!Model core charge for nonlinear core xc correction, and 4 derivatives

 read (tmp_unit,*) rchrg,fchrg,qchrg
 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
 call wrtout(ab_out,message,'COLL')
 call wrtout(std_out,  message,'COLL')

 ABI_ALLOCATE(nproj_tmp,(mpssoang))

 nproj_tmp(:)=0
 read (tmp_unit,*) nproj_tmp(1:lmax+1)
 write(message, '(a,5i6)' ) '     nproj',nproj_tmp(1:lmax+1)
 call wrtout(ab_out,message,'COLL')
 call wrtout(std_out,  message,'COLL')

!place holder for future implementation of additional optional header
!lines without invalidating existing psp files
!Now (12/2014) extended to include spin-orbit projectors
 read (tmp_unit,*) extension_switch
 if(extension_switch==2) then
   read (tmp_unit,*) nproj_tmp(lmax+2:2*lmax+1)
   write(message, '(a,i6)' ) 'spin-orbit psp, extension_switch',extension_switch
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,  message,'COLL')
   write(message, '(a,5i6)' ) '   nprojso',nproj_tmp(lmax+2:2*lmax+1)
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,  message,'COLL')
 elseif(extension_switch==0) then
   write(message, '(a,i6)' ) 'extension_switch',extension_switch
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,  message,'COLL')
 else
   write(message, '(a,i6/a)' ) 'invalid extension_switch',extension_switch,&
&   'Should be 0 for scalar-relativisgtic psp or 2 to include spin-orbit'
   MSG_ERROR(message)
 end if

 if(lloc<4) then
   if (nproj_tmp(lloc+1)/=0) then
     write(message, '(a,i4,a,a,i4,5a)' )&
&     '  Pseudopotential input file has nproj=',nproj_tmp(lloc+1),ch10,&
&     '  for angular momentum',lloc,' which is the local potential.',ch10,&
&     '  Should be 0 for the local potential',ch10,&
&     '  Action : check your pseudopotential input file.'
     MSG_ERROR(message)
   end if
 end if

!--------------------------------------------------------------------

!Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn
 if(pspso==2) then
   nso=2
 else
   nso=1
 end if

 index=0;iln=0;indlmn(:,:)=0
 do nn=1,nso
   do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
     ll=ipsang-(nn-1)*lmax-1
     if (nproj_tmp(ipsang)>0) then
       do kk=1,nproj_tmp(ipsang)
         iln=iln+1
         do mm=1,2*ll*useylm+1
           index=index+1
           indlmn(1,index)=ll
           indlmn(2,index)=mm-ll*useylm-1
           indlmn(3,index)=kk
           indlmn(4,index)=ll*ll+(1-useylm)*ll+mm
           indlmn(5,index)=iln
           indlmn(6,index)=nn
         end do
       end do
     end if
   end do
 end do

! repackage nproj_tmp for proper use by pspatm
 nproj(:)=0
 nproj(1:lmax+1)=nproj_tmp(1:lmax+1)
 if(pspso==2) then
   nproj(mpsang+1:mpsang+lmax)=nproj_tmp(lmax+2:2*lmax+1)
 end if

 ABI_ALLOCATE(rad,(mmax))
 ABI_ALLOCATE(vloc,(mmax))
 ABI_ALLOCATE(vpspll,(mmax,lnmax))

!Will now proceed at the reading of pots and projectors

!rad(:)=radial grid r(i)
!vpspll(:,1),...,vpspll(:,lnmax)=nonlocal projectors
!vloc(:)=local potential

!Read Vanderbilt-Kleinman-Bylander energies and projectors for each l
!or read local potential for l=lloc.
!Also get rad array (actually read more than once)
 ll_err=0
 iln0=0
 do nn=1,nso
   do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
     ll=ipsang-(nn-1)*lmax-1
     if (nproj_tmp(ipsang)>0) then
       read(tmp_unit,*) llin,ekb(iln0+1:iln0+nproj_tmp(ipsang))
       if(llin/=ll) then
         ll_err=ipsang
         exit
       end if
       do irad=1,mmax
         read(tmp_unit,*)jj,rad(irad),vpspll(irad,iln0+1:iln0+nproj_tmp(ipsang))
       end do
       iln0=iln0+nproj_tmp(ipsang)
     elseif(ll==lloc .and. nn==1) then
       read(tmp_unit,*) llin
       if(llin/=ll) then
         ll_err=ipsang
         exit
       end if
       do irad=1,mmax
         read(tmp_unit,*)jj,rad(irad),vloc(irad)
       end do
     end if
   end do !ipsang
!Provision for general local potential /= any angular momentum potential
   if(nn==1 .and. lloc>lmax) then
     read(tmp_unit,*) llin
     if(llin==lloc) then
       do irad=1,mmax
         read(tmp_unit,*)jj,rad(irad),vloc(irad)
       end do
     else
       ll_err=lloc+1
       exit
     end if
   end if
 end do !nn

 if(ll_err>0) then
   write(message, '(5a,i4,a,i4,a,a)' )&
&   'Pseudopotential input file does not have angular momenta in order',ch10,&
&   'or has inconsistent general local potential index',ch10,&
&   'Expected',ll_err-1,' , got',ll,ch10,&
&   'Action : check your pseudopotential input file.'
   MSG_ERROR(message)
 end if

!Check that rad grid is linear starting at zero
 amesh=rad(2)-rad(1)
 damesh=zero
 do irad=2,mmax-1
   damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
 end do
 if(damesh>tol8 .or. rad(1)/=zero) then
   write(message, '(5a)' )&
&   'Pseudopotential input file reqires linear radial mesh',ch10,&
&   'starting at zero.',ch10,&
&   'Action : check your pseudopotential input file.'
   MSG_ERROR(message)
 end if

!Get core charge function and derivatives, if needed
 if(fchrg>1.0d-15)then
   call psp8cc(mmax,n1xccc,rchrg,xccc1d)
!  The core charge function for pspcod=8
!  becomes zero beyond rchrg. Thus xcccrc must be set
!  equal to rchrg .
   xcccrc=rchrg
 else
   xccc1d(:,:)=0.0d0
   xcccrc=0.0d0
 end if

!--------------------------------------------------------------------
!Carry out calculations for local (lloc) pseudopotential.
!Obtain Fourier transform (1-d sine transform)
!to get q^2 V(q).

 call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,&
& vlspl(:,1),rad,vloc,yp1,ypn,zion)


!Fit spline to q^2 V(q) (Numerical Recipes subroutine)
 ABI_ALLOCATE(work_space,(mqgrid))
 ABI_ALLOCATE(work_spl,(mqgrid))
 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
 vlspl(:,2)=work_spl(:)
 ABI_DEALLOCATE(work_space)
 ABI_DEALLOCATE(work_spl)

!--------------------------------------------------------------------
!Take care of non-local part

!Allow for option of no nonlocal corrections (lloc=lmax=0)
 if (lloc==0.and.lmax==0) then
   write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
   call wrtout(ab_out,message,'COLL')
   call wrtout(std_out,  message,'COLL')
 else

!  ----------------------------------------------------------------------
!  Compute Vanderbilt-KB form factors and fit splines

   call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,&
&   mqgrid,qgrid,rad,vpspll)

 end if

 ABI_DEALLOCATE(nproj_tmp)
 ABI_DEALLOCATE(vpspll)
 ABI_DEALLOCATE(rad)
 ABI_DEALLOCATE(vloc)

end subroutine psp8in
!!***
