TABLE OF CONTENTS


ABINIT/m_psp9 [ Modules ]

[ Top ] [ Modules ]

NAME

 m_psp9

FUNCTION

 Initialize pspcod=9 (pseudopotentials from the PSML XML format):

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (JJ, MVer, YP)
  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_psp9
28 
29  use defs_basis
30  use m_splines
31  use m_errors
32  use m_abicore
33 #if defined HAVE_PSML
34  use m_psml
35 #endif
36 
37  use defs_datatypes,  only : nctab_t
38  use m_pawrad,        only : pawrad_type, pawrad_init, pawrad_free
39  use m_psps,          only : nctab_eval_tvalespl
40  use m_psptk,         only : psp8lo, psp8nl
41 
42  implicit none
43 
44  private

ABINIT/psp9cc [ Functions ]

[ Top ] [ Functions ]

NAME

 psp9cc

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for format 9 of the pseudopotentials (PSML).

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

OUTPUT

  rchrg=cut-off radius for the core density
  xccc1d(n1xccc,6)= 1D core charge function and its four first derivatives

NOTES

  This routine will be built only if PSML support is enabled.

PARENTS

      psp9in

CHILDREN

      dgesv

SOURCE

709 #if defined HAVE_PSML
710 
711 subroutine psp9cc(psxml,mmax,n1xccc,rad,rchrg,xccc1d)
712 
713 
714 !This section has been created automatically by the script Abilint (TD).
715 !Do not modify the following lines by hand.
716 #undef ABI_FUNC
717 #define ABI_FUNC 'psp9cc'
718 !End of the abilint section
719 
720  implicit none
721 
722 !Arguments ------------------------------------
723 !scalars
724  integer,intent(in) :: mmax,n1xccc
725  real(dp),intent(out) :: rchrg
726  type(ps_t),intent(in) :: psxml
727 !arrays
728  real(dp),intent(in) :: rad(mmax)
729  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
730 
731 !Local variables-------------------------------
732 !scalars
733  integer :: i1xccc,idum,irad,jj
734  real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx,twelvth
735  character(len=500) :: message
736 !arrays
737  integer :: iwork(8)
738  real(dp) :: rscale(5),dpoly(6,6),vpoly(6)
739  real(dp),allocatable :: ff(:,:)
740 
741 !**********************************************************************
742 
743 !Check that rad grid is linear starting at zero
744  amesh=rad(2)-rad(1)
745  damesh=zero
746  do irad=2,mmax-1
747    damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
748  end do
749 
750  if(damesh>tol8 .or. rad(1)/=zero) then
751    write(message, '(5a)' )&
752 &   'Pseudopotential input file requires linear radial mesh',ch10,&
753 &   'starting at zero.',ch10,&
754 &   'Action: check your pseudopotential input file.'
755    MSG_ERROR(message)
756  end if
757 
758  ABI_ALLOCATE(ff,(mmax,5))
759 
760  dri = one / amesh
761  pi4i = quarter / pi
762  twelvth = one / 12.0_dp
763 
764 !Read from pp file the model core charge and calculate its first 4 derivatives
765 !assumed to be on a linear grid starting at zero.
766 !The input functions contain the 4pi factor, and must be rescaled.
767 
768 !Store the value of the pseudo-core charge.
769  ff(:,:) = zero
770  do jj=1,mmax
771    ff(jj,1) = ps_CoreCharge_Value(psxml,rad(jj))
772  end do
773 
774 !Calculate 4 first derivatives with 5-point stencil, except borders
775  do irad=3,mmax-2
776    ff(irad,2) = (-ff(irad+2,1) + 8.0d0*ff(irad+1,1) - &
777 &   8.0d0*ff(irad-1,1) + ff(irad-2,1)) * twelvth * dri
778    ff(irad,3) = (-ff(irad+2,1) + 16.0d0*ff(irad+1,1) - 30.0d0*ff(irad,1) + &
779 &   16.0d0*ff(irad-1,1) - ff(irad-2,1)) * twelvth * dri * dri
780    ff(irad,4) = (ff(irad+2,1) - 2.0d0*ff(irad+1,1) + &
781 &   2.0d0*ff(irad-1,1) - ff(irad-2,1)) * half * dri * dri * dri
782    ff(irad,5) = (ff(irad+2,1) - 4.0d0*ff(irad+1,1) + 6.0d0*ff(irad,1) - &
783 &   4.0d0*ff(irad-1,1) + ff(irad-2,1)) * dri * dri * dri * dri
784  end do
785 
786 !Add border near zero using polynomial fit
787  dpoly(:,:) = zero
788  dpoly(:,1) = one
789  vpoly(:) = zero
790  vpoly(1) = ff(1,1)
791  do irad=2,6
792    do jj=1,6
793      dpoly(irad,jj) = rad(irad)**(jj-1)
794    end do
795    vpoly(irad) = ff(irad,1)
796  end do
797  call dgesv(6,1,dpoly,6,iwork,vpoly,6,idum)
798 
799  do irad=1,2
800    ff(irad,2) = &
801 &   vpoly(2) + 2.0d0*vpoly(3)*rad(irad) + &
802 &   3.0d0*vpoly(4)*rad(irad)*rad(irad) + &
803 &   4.0d0*vpoly(5)*rad(irad)*rad(irad)*rad(irad) + &
804 &   5.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)*rad(irad)
805    ff(irad,3) = &
806 &   2.0d0*vpoly(3)*rad(irad) + &
807 &   6.0d0*vpoly(4)*rad(irad) + &
808 &   12.0d0*vpoly(5)*rad(irad)*rad(irad) + &
809 &   20.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)
810    ff(irad,4) = &
811 &   6.0d0*vpoly(4) + &
812 &   24.0d0*vpoly(5)*rad(irad) + &
813 &   60.0d0*vpoly(6)*rad(irad)*rad(irad)
814    ff(irad,5) = 24.0d0*vpoly(5) + &
815 &   120.0d0*vpoly(6)*rad(irad)
816  end do
817 
818 !Make linear approximation for the tail near mmax
819  do irad=1,2
820    ff(mmax-2+irad,2) = ff(mmax-2,2) + irad * (ff(mmax-2,2) - ff(mmax-3,2))
821    ff(mmax-2+irad,3) = ff(mmax-2,3) + irad * (ff(mmax-2,3) - ff(mmax-3,3))
822    ff(mmax-2+irad,4) = ff(mmax-2,4) + irad * (ff(mmax-2,4) - ff(mmax-3,4))
823    ff(mmax-2+irad,5) = ff(mmax-2,5) + irad * (ff(mmax-2,5) - ff(mmax-3,5))
824  end do
825 
826 !Renormalize core charge
827 ! ff(:,:) = ff(:,:) * pi4i
828 
829 !determine xcccrc where the pseudocore becomes 0
830 !This is a difference with respect the Hamann's treatment of the core
831 !charge when reading PSP8.
832 !In Hamann's case (PSP8), xcccrc = rchrg, and this value is
833 !introduced in the pseudopotential input file.
834 !rchrg is not included in the PSML format
835  rchrg = zero
836  do jj=mmax,1,-1
837    if (ff(jj,1) > tol13) then
838      rchrg=rad(jj)
839      exit
840    end if
841  end do
842 
843 !Check that input rchrg is consistent with last grid point
844  if(rchrg>rad(mmax)) then
845    write(message, '(5a)' )&
846 &   'Pseudopotential input file core charge mesh',ch10,&
847 &   'is inconsistent with rchrg in header.',ch10,&
848 &   'Action: check your pseudopotential input file.'
849    MSG_ERROR(message)
850  end if
851 
852 !Factors for unit range scaling
853  do jj = 1, 5
854    rscale(jj)=rchrg**(jj-1)
855  end do
856 
857 !Generate uniform mesh xx in the box cut by rchrg
858 !and interpolate the core charge and derivatives
859 !Cubic polynomial interpolation is used which is consistent
860 !with the original interpolation of these functions from
861 !a log grid to the input linear grid.
862 
863  dri=1.d0/amesh
864  do i1xccc=1,n1xccc
865    xx=(i1xccc-1)* rchrg/dble(n1xccc-1)
866 
867 !  index to find bracketing input mesh points
868    irad = int(dri * xx) + 1
869    irad = max(irad,2)
870    irad = min(irad,mmax-2)
871 !  interpolation coefficients
872    xp = dri * (xx - rad(irad))
873    xpp1 = xp + one
874    xpm1 = xp - one
875    xpm2 = xp - two
876    c1 = -xp * xpm1 * xpm2 * sixth
877    c2 = xpp1 * xpm1 * xpm2 * half
878    c3 = - xp * xpp1 * xpm2 * half
879    c4 = xp * xpp1 * xpm1 * sixth
880 !  Now do the interpolation on all derivatives for this grid point
881 !  Include 1/4pi normalization and unit range scaling
882    do jj=1,5
883      tff =  c1 * ff(irad - 1, jj) &
884 &     + c2 * ff(irad    , jj) &
885 &     + c3 * ff(irad + 1, jj) &
886 &     + c4 * ff(irad + 2, jj)
887      xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff
888    end do
889  end do
890 
891 !5th derivative is apparently not in use, so set to zero
892  xccc1d(:,6)=zero
893 
894  ABI_DEALLOCATE(ff)
895 
896 end subroutine psp9cc

ABINIT/psp9in [ Functions ]

[ Top ] [ Functions ]

NAME

 psp9in

FUNCTION

 Initialize pspcod=9 (pseudopotentials from the PSML XML format):
 continue to read the corresponding file, then compute the
 local and non-local potentials.

INPUTS

  filpsp=filename of the PSML pseudopotential
  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,
             {{\ \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)
  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,ps_corecharge_get
      ps_destroy,ps_nonlocalprojectors_filter,ps_projector_get
      ps_provenance_get,ps_pseudoatomspec_get,ps_valenceconfiguration_get
      ps_valenceshell_get,psml_reader,psp8lo,psp8nl,psp9cc,spline,wrtout

SOURCE

120 subroutine psp9in(filpsp,ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
121 &                  mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,&
122 &                  useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad)
123 
124 
125 !This section has been created automatically by the script Abilint (TD).
126 !Do not modify the following lines by hand.
127 #undef ABI_FUNC
128 #define ABI_FUNC 'psp9in'
129 !End of the abilint section
130 
131  implicit none
132 
133 !Arguments ------------------------------------
134 !scalars
135  integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mpsang,mpssoang,mqgrid,mqgrid_vl
136  integer,intent(in) :: pspso,n1xccc,useylm
137  integer,intent(out) :: mmax
138  real(dp),intent(in) :: zion,znucl
139  real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad
140  type(nctab_t),intent(inout) :: nctab
141  character(len=fnlen),intent(in) :: filpsp
142 !arrays
143  integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang)
144  real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl)
145  real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2)
146  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
147 
148 !Local variables-------------------------------
149 !scalars
150 #if defined HAVE_PSML
151  integer :: iln,pspindex,ipsang,irad,kk,ll
152  integer :: mm,nn,nso,ii,ir,il
153  integer :: nshells
154  integer :: iproj,irelt,nders
155  integer :: np_dn, np_lj, np_nr, np_so, np_sr, np_up, val_l, val_n
156  real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn,zval
157  real(dp) :: rmax,rmatch,z,chgvps
158  real(dp) :: val_occ
159  logical :: has_nlcc,has_spin
160  logical :: has_tvale,oncvpsp
161  character(len=500) :: message
162  character(len=30)  :: creator
163  character(len=7), parameter  :: oncvpsp_name = "ONCVPSP"
164  type(pawrad_type) :: mesh
165 #endif
166 !arrays
167 #if defined HAVE_PSML
168  integer, allocatable :: idx_so(:),idx_sr(:)
169  real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_spl(:)
170  type(ps_t) :: psxml
171 #endif
172 
173 ! ***************************************************************************
174 
175 #if defined HAVE_PSML
176 
177  call ps_destroy(psxml)
178  call psml_reader(filpsp,psxml,debug=.true.)
179 
180 !Identify the atomic code that generated the pseudopotential
181  call ps_Provenance_Get(psxml, 1, creator=creator)
182 !Check whether the pseudopotential has been created with ONCVPSP,
183 !Don Hamann's code
184  oncvpsp = (trim(creator(1:7)) .eq. trim(oncvpsp_name))
185 !DEBUG
186 !write(std_out,*)' psp9in : creator : ', creator
187 !write(std_out,*)' psp9in : oncvpsp : ', oncvpsp
188 !ENDDEBUG
189 
190 ! SIESTA's ATOM uses spherical harmonics, while ONCVPSP uses Legendre
191 ! polynomials, which means we have to check the consistency of input variables
192 ! wrt the pseudos
193 !
194 ! Note: commented because NC pseudos do not have non-diagonal terms
195 !
196 ! if ( oncvpsp ) then
197 !   if ( useylm /= 0 ) then
198 !     write(message,'(3a)') "ONCVPSP pseudos use Legendre polynomials but we use spherical harmonics", &
199 !&      ch10, "ACTION: set useylm to 0 in your input file"
200 !     MSG_ERROR(message)
201 !   endif
202 ! else
203 !   if ( useylm == 0 ) then
204 !     write(message,'(3a)') "ATOM pseudos use spherical harmonics but we use Legendre polynomials", &
205 !&      ch10, "ACTION: set useylm to 1 in your input file"
206 !     MSG_ERROR(message)
207 !   endif
208 ! endif
209 
210 ! The atomic number is a real number instead of a simple integer
211 ! z (in Abinit), atomic-number in the header of the PSML file.
212 ! z      = ps_AtomicNumber(psxml)
213 !
214 ! The difference between the number of protons in the nucleus and the
215 ! sum of the populations of the core shells is the effective atomic number
216 ! of the pseudo-atom, Zval (in Abinit), z-pseudo in the header of the
217 ! PSML file.
218 ! zval   = ps_Zpseudo(psxml)
219 
220  call ps_PseudoAtomSpec_Get(psxml, &
221 & atomic_number=z, z_pseudo=zval, &
222 & spin_dft=has_spin, core_corrections=has_nlcc)
223 
224 !---
225 
226 !Feb 2015: shifted to Hamann grid for convenience - libpsml interpolates anyway
227 !
228 ! The following lines are taken from the oncvpsp.f90 subroutine of the oncvpsp
229 ! code implemented by D. Hamann
230 ! The atomic number of the element is read from the header of the XML file
231 ! Logarithmic grid defined by Hamann in oncvpsp code
232 ! z    = psxml%header%z
233 ! amesh = 1.012d0
234 ! al    = dlog(amesh)
235 ! rr1   = .0005d0/z
236 ! mmax  = dlog(45.0d0 /rr1)/al
237 !
238 ! ABI_ALLOCATE( rad,(mmax) )
239 !
240 ! do ir = 1, mmax
241 !   rad(ir) = rr1 * dexp(al*(ir-1))
242 ! end do
243 
244 !Determine the maximum number of points in the grid ---
245  rmax  = 6.0_dp
246  amesh = 0.01_dp
247  mmax  = int(rmax/amesh)
248 ! if(mod(mmax,2) .eq. 0) mmax = mmax + 1
249 
250 !Print core charge info, for compatibility with psp8
251  rchrg  = zero
252  fchrg  = zero
253  if (has_nlcc) then
254    rchrg = amesh * (mmax - 2)
255 !  PSML does not store fchrg for now but we know we have core corrections,
256 !  then let's set it arbitrarily to 1.0
257    fchrg = one
258  else
259    write(message, '(a)' ) '- psp9in: No XC core correction.'
260    call wrtout(std_out,message,'COLL')
261  end if
262  write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,zero,'rchrg,fchrg,qchrg'
263  call wrtout(ab_out,message,'COLL')
264  call wrtout(std_out,message,'COLL')
265 
266 !Do we have a valence charge?
267  call ps_ValenceConfiguration_Get(psxml, nshells=nshells)
268  has_tvale = (nshells > 0)
269 
270 ! Compute the valence charge of the reference configuration used to
271 ! generate the pseudopotential
272  chgvps = 0.0_dp
273 ! Loop on all the shells included in the valence
274  do il = 1, nshells
275 !  Sum the corresponding occupation of each shell
276 !  FIXME: What if there is spin?
277    call ps_ValenceShell_Get(psxml, il, n=val_n, l=val_l, occupation=val_occ)
278    chgvps = chgvps + val_occ
279    write(std_out,*)' psp9in : n, l, occupation = ',   &
280 &   val_n, val_l, val_occ
281  end do
282 
283 !DEBUG
284 !write(std_out,*)' psp9in : atomic number'
285 !write(std_out,*)' psp9in :   z = ', z
286 !write(std_out,*)' psp9in : valence charge of the reference configuration'
287 !write(std_out,*)' psp9in :   chgvps = ', chgvps
288 !write(std_out,*)' psp9in : nominal valence charge'
289 !write(std_out,*)' psp9in :   zval = ', zval
290 !write(std_out,*)' psp9in :   mqgrid_vl = ', mqgrid_vl
291 !write(std_out,*)' psp9in : parameters to define the points of the grid'
292 !write(std_out,*)' psp9in :   amesh = ', amesh
293 !write(std_out,*)' psp9in :   rmax = ', rmax
294 !write(std_out,*)' psp9in :   mmax = ', mmax
295 !ENDDEBUG
296 
297 ! TODO: should be simple to average these and get difference for SREL+SOC,
298 ! but also the Ekb etc...
299  call ps_NonlocalProjectors_Filter(psxml, set=SET_DOWN, number=np_dn)
300  call ps_NonlocalProjectors_Filter(psxml, set=SET_LJ, number=np_lj)
301  call ps_NonlocalProjectors_Filter(psxml, set=SET_NONREL, number=np_nr)
302  call ps_NonlocalProjectors_Filter(psxml, set=SET_SO, number=np_so)
303  call ps_NonlocalProjectors_Filter(psxml, set=SET_SREL, number=np_sr)
304  call ps_NonlocalProjectors_Filter(psxml, set=SET_UP, number=np_up)
305  if (np_lj > 0) then
306    message = 'For the moment LJ format projectors are not supported; SREL + SO is the internal abinit format'
307    MSG_BUG(message)
308  end if
309 
310  if (np_up > 0 .or. np_dn > 0) then
311    write (message,'(3a)') 'For the moment separate spin up and down format projectors are not supported;',ch10,&
312 &   ' spin average is the internal abinit format'
313    MSG_BUG(message)
314  end if
315 
316 !--------------------------------------------------------------------
317 
318 !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn
319  if(pspso==2) then
320    nso=2
321  else
322    nso=1
323  end if
324 
325 !Find the number of projectors per angular momentum shell
326  nproj(:)=0
327  if (np_nr > 0) then
328    call ps_NonlocalProjectors_Filter(psxml, set=SET_NONREL, indexes=idx_sr)
329    do iproj = 1, np_nr
330      call ps_Projector_Get(psxml, idx_sr(iproj), l=il)
331      nproj(il+1) = nproj(il+1) + 1
332    end do
333  else
334    if (np_sr > 0) then
335      call ps_NonlocalProjectors_Filter(psxml, set=SET_SREL, indexes=idx_sr)
336      do iproj = 1, np_sr
337        call ps_Projector_Get(psxml, idx_sr(iproj), l=il)
338        nproj(il+1) = nproj(il+1) + 1
339      end do
340    else ! this should not happen
341      MSG_BUG('Your psml potential should have either scalar- or non- relativistic projectors')
342    end if
343  end if
344 
345  write(message, '(a,5i6)' ) '     nproj',nproj(1:lmax+1)
346  call wrtout(ab_out,message,'COLL')
347  call wrtout(std_out,  message,'COLL')
348 
349  irelt = 0
350  if (nso == 2) then
351    call ps_NonlocalProjectors_Filter(psxml, set=SET_SO, indexes=idx_so)
352    do iproj = 1, np_so
353      call ps_Projector_Get(psxml, idx_so(iproj), l=il)
354      nproj(il+lmax+2) = nproj(il+lmax+2) + 1
355      irelt = 1
356    end do
357  end if
358 
359  pspindex=0;iln=0;indlmn(:,:)=0
360  do nn=1,nso
361    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
362      ll=ipsang-(nn-1)*lmax-1
363      if (nproj(ipsang)>0) then
364        do kk=1,nproj(ipsang)
365          iln=iln+1
366          do mm=1,2*ll*useylm+1
367            pspindex=pspindex+1
368            indlmn(1,pspindex)=ll                      ! l angular momentum channel
369            indlmn(2,pspindex)=mm-ll*useylm-1          ! hash of position in m
370            indlmn(3,pspindex)=kk                      ! index of projector
371            indlmn(4,pspindex)=ll*ll+(1-useylm)*ll+mm  ! hash of position in l(l+1) array
372            indlmn(5,pspindex)=iln                     ! absolute index of l, n disregarding m values
373            indlmn(6,pspindex)=nn                      ! spin orbit index!!! NOT the n shell index
374          end do
375        end do
376      end if
377    end do
378  end do
379 
380 ! Determine whether the atomic calculation to generate the pseudopotential
381 ! is relativistic or not
382 
383 !DEBUG
384 !write(std_out,*)' psp9in : pseudopotential generation relativity ', ps_Relativity(psxml)
385 !write(std_out,*)' psp9in : SOC pseudopotential? (1=yes, 0 =no) '
386 !write(std_out,*)' psp9in : irelt = ', irelt
387 !write(ab_out,*)' psp9in : irelt = ', irelt
388 !ENDDEBUG
389 
390 !Can now allocate grids, potentials and projectors
391  ABI_ALLOCATE(rad,(mmax))
392  ABI_ALLOCATE(vloc,(mmax))
393  ABI_ALLOCATE(vpspll,(mmax,lnmax))
394 
395 !Feb 2015: shifted to Hamann grid for convenience - libpsml interpolates anyway
396  do ir=1,mmax
397    rad(ir) = amesh * (ir - 1)
398  end do
399 !! DEBUG
400 ! do ir = 2, mmax
401 !   write(std_out,'(i5,f20.12)')ir, rad(ir)
402 ! end do
403 !! ENDDEBUG
404 !---
405  write(message, '(a,i5,es16.6,es16.6)')'  psp9in : mmax, amesh, rad(mmax) = ', mmax, amesh, rad(mmax)
406  call wrtout(ab_out,message,'COLL')
407  call wrtout(std_out,message,'COLL')
408 
409 !Check that rad grid is linear starting at zero
410  amesh=rad(2)-rad(1)
411  damesh=zero
412  do irad=2,mmax-1
413    damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
414  end do
415  if(damesh>tol8 .or. rad(1)/=zero) then
416    write(message, '(5a)' )&
417 &   'Pseudopotential input file requires linear radial mesh',ch10,&
418 &   'starting at zero.',ch10,&
419 &   'Action: check your pseudopotential input file.'
420    MSG_ERROR(message)
421  end if
422 
423 !Take care of the non-linear core corrections
424 !----------------------------------------------------------------------------
425 ! xcccrc           : XC core correction cutoff radius (bohr)
426 !                    It is defined as the radius where the pseudo-core
427 !                    charge density becomes zero
428 !                    (here we have set up a tolerance of 1.d-12).
429 
430  rmatch = zero
431  nders  = 0
432  if (has_nlcc) then
433 
434 !    In Abinit, at least for the Troullier-Martins pseudopotential,
435 !    the pseudocore charge density and its derivatives (xccc1d)
436 !    are introduced in a linear grid.
437 !    This grid is normalized, so the radial coordinates run between
438 !    from 0 and 1 (from 0 to xcccrc, where xcccrc is the radius
439 !    where the pseudo-core becomes zero).
440 
441    call ps_CoreCharge_get(psxml, rc=rmatch, nderivs=nders)
442    write (message,'(1X,A,A,5X,A,1X,F8.3,A,5X,A,I8,A)') &
443 &   "Reading pseudocore charge",ch10, &
444 &   "- matching radius:",rmatch,ch10, &
445 &   "- number of continuous derivatives",nders,ch10
446    call wrtout(std_out,message,'COLL')
447 
448 !Get core charge function and derivatives, if needed
449    if(fchrg>1.0d-15)then
450      call psp9cc(psxml,mmax,n1xccc,rad,rchrg,xccc1d)
451 !  The core charge function for pspcod=9
452 !  becomes zero beyond rchrg. Thus xcccrc must be set
453 !  equal to rchrg.
454      xcccrc=rchrg
455    else
456      xccc1d(:,:) = zero
457      xcccrc = zero
458      fchrg = zero
459      qchrg = zero
460    end if
461 
462    maxrad = rad(mmax)
463 
464  end if ! has_nlcc
465 
466 !!   DEBUG
467 !    write(std_out,*)' xcccrc = ', xcccrc, rchrg
468 !    write(std_out,*)
469 !    write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc
470 !    do ii = 1, n1xccc
471 !    write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),&
472 ! &         xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6)
473 !    enddo
474 !    write(std_out,*)
475 !    stop
476 !!   ENDDEBUG
477 
478 
479 !--------------------------------------------------------------------
480 !Carry out calculations for local (lloc) pseudopotential.
481 !Obtain Fourier transform (1-d sine transform)
482 !to get q^2 V(q).
483 
484 !Read and process vlocal:
485 !The local potential is given by a <radfunc> element under the <local-potential>
486 !element.
487 !After reading, this is a copy of the treatment to the
488 !local part carry out in psp8
489 !i.e. (as in Hamann pseudopotential)
490 !
491 !Read the local component of the pseudopotential
492  vloc = zero
493  do ir = 1, mmax
494    vloc(ir) = ps_LocalPotential_Value(psxml, rad(ir))
495  end do
496 
497  call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,&
498 & vlspl(:,1),rad,vloc,yp1,ypn,zion)
499 
500 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
501  ABI_ALLOCATE(work_spl,(mqgrid))
502  call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
503  vlspl(:,2)=work_spl(:)
504  ABI_DEALLOCATE(work_spl)
505 
506 !!  DEBUG
507 ! write(std_out,*)'# Vlocal = '
508 ! write(std_out,*)' amesh  = ', amesh
509 ! write(std_out,*)' epsatm = ', epsatm
510 ! write(std_out,*)' mmax   = ', mmax
511 ! write(std_out,*)' mqgrid = ', mqgrid
512 ! do ir = 1, mqgrid
513 !   write(std_out,*)'   qgrid = ', ir, qgrid(ir)
514 ! enddo
515 ! do ir = 1, mqgrid
516 !   write(std_out,'(a,i5,2f20.12)')'   iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2)
517 ! enddo
518 ! write(std_out,*)
519 ! do ir = 1, mmax
520 !   write(std_out,*)'   rad   = ', rad(ir), vloc(ir)
521 ! enddo
522 ! write(std_out,*)
523 ! write(std_out,*)' yp1    = ', yp1
524 ! write(std_out,*)' ypn    = ', ypn
525 ! write(std_out,*)' zion   = ', zion
526 ! stop
527 !!  ENDDEBUG
528 
529 
530 !--------------------------------------------------------------------
531 !Take care of non-local part
532 
533 !Zero out all Kleinman-Bylander energies to initialize
534  do ii = 1, lmnmax ! loop over all possible projectors
535    if (indlmn(6,ii) == 1) then
536      call ps_Projector_Get(psxml, idx_sr(indlmn(5,ii)), ekb=ekb(indlmn(5,ii)))
537    else if (indlmn(6,ii) == 2) then
538      call ps_Projector_Get(psxml, idx_so(indlmn(5,ii)), ekb=ekb(indlmn(5,ii)))
539    end if
540  end do
541 
542 !Read the KB projectors from the PSML file
543 !Note than in the PSML file the radial part of the projector is stored,
544 !while Abinit expects the radial part times the radii.
545 !We have to multiply by r after reading it.
546 !Note than in Hamann's format (psp8), Abinit directly reads r * radial_part_KB
547  vpspll = zero
548  do ii = 1, lmnmax
549    if (indlmn(6,ii) == 1) then
550      do ir = 1, mmax
551        vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_sr(indlmn(5,ii)), rad(ir))
552        vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii))
553      end do
554    else if (indlmn(6,ii) == 2) then
555      do ir = 1, mmax
556        vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_so(indlmn(5,ii)), rad(ir))
557        vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii))
558      end do
559    end if
560  end do
561 
562 !Allow for option of no nonlocal corrections (lloc=lmax=0)
563  if (lloc==0.and.lmax==0) then
564    write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
565    call wrtout(ab_out,message,'COLL')
566    call wrtout(std_out,message,'COLL')
567  else
568 
569 !  ----------------------------------------------------------------------
570 !  Compute Vanderbilt-KB form factors and fit splines
571 
572    call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,&
573 &   mqgrid,qgrid,rad,vpspll)
574 
575  end if
576 
577 !!  DEBUG
578 ! write(std_out,*)'# KB Projectors = '
579 ! write(std_out,*)' amesh  = ', amesh
580 ! do ir = 1, mqgrid
581 !   do il = 1, lnmax
582 !     write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il)
583 !   enddo
584 ! enddo
585 ! do il = 1, lmnmax
586 !   write(std_out,*)' indlmn = ', il, indlmn(:,il)
587 ! enddo
588 ! write(std_out,*)' lmax   = ', lmax
589 ! write(std_out,*)' lmnmax = ', lmnmax
590 ! write(std_out,*)' lnmax  = ', lnmax
591 ! write(std_out,*)' mmax   = ', mmax
592 ! write(std_out,*)' mqgrid = ', mqgrid
593 ! do ir = 1, mqgrid
594 !   write(std_out,*)'   qgrid = ', ir, qgrid(ir)
595 ! enddo
596 ! do il = 1, lnmax
597 !   write(std_out,*)
598 !   write(std_out,*)'# il = ', il
599 !   do ir = 1, mmax
600 !     write(std_out,*)'   rad   = ', rad(ir), vpspll(ir,il)
601 !   enddo
602 ! enddo
603 ! stop
604 !!  ENDDEBUG
605 
606 ! Read pseudo valence charge in real space on the linear mesh
607 ! and transform it to reciprocal space on a regular grid. Use vloc as workspace.
608  vloc(:) = zero
609  if (has_tvale) then
610    do irad=1,mmax
611      vloc(irad) = ps_ValenceCharge_Value(psxml,rad(irad))
612      vloc(irad) = vloc(irad) / four_pi
613    end do
614 
615 !! DEBUG
616 !  do irad = 1, mmax
617 !    write(std_out,*)' Valence Charge  = ', rad(irad), vloc(irad)
618 !  enddo
619 !  stop
620 !! ENDDEBUG
621 
622 
623    ! Check that rad grid is linear starting at zero
624    amesh=rad(2)-rad(1)
625    damesh=zero
626    do irad=2,mmax-1
627      damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1)))
628    end do
629    if(damesh>tol8 .or. rad(1)/=zero) then
630      write(message, '(5a)' )&
631 &     'Pseudopotential input file requires linear radial mesh',ch10,&
632 &     'starting at zero.',ch10,&
633 &     'Action: check your pseudopotential input file.'
634      MSG_ERROR(message)
635    end if
636 
637    !  Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space.
638    call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh)
639    call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl)
640    call pawrad_free(mesh)
641  end if
642 
643  ABI_DEALLOCATE(vpspll)
644  ABI_DEALLOCATE(vloc)
645  ABI_DEALLOCATE(rad)
646  if (allocated(idx_sr)) then
647    ABI_FREE_NOCOUNT(idx_sr)
648  end if
649  if (allocated(idx_so)) then
650    ABI_FREE_NOCOUNT(idx_so)
651  end if
652 
653  call ps_destroy(psxml)
654 
655 !--------------------------------------------------------------------
656 
657 #else
658  ABI_UNUSED(mpsang)
659  ABI_UNUSED(pspso)
660  ABI_UNUSED(qgrid_vl)
661  ABI_UNUSED(nctab%mqgrid_vl)
662 !Initialize some arguments, for portability at compile time
663  indlmn=0 ; mmax=0 ; nproj=0
664  ekb=zero ; epsatm=zero ; ffspl=zero ; qchrg=zero ; vlspl=zero ; xcccrc=zero ; xccc1d=zero
665 
666  if(.false.)write(std_out,*)filpsp ! Just to keep filpsp when HAVE_PSML is false
667  if(.false.)write(std_out,*)lloc   ! Just to keep lloc when HAVE_PSML is false
668  if(.false.)write(std_out,*)lmax   ! Just to keep lmax when HAVE_PSML is false
669  if(.false.)write(std_out,*)mpsang ! Just to keep mpsang when HAVE_PSML is false
670  if(.false.)write(std_out,*)pspso  ! Just to keep pspso when HAVE_PSML is false
671  if(.false.)write(std_out,*)qgrid  ! Just to keep qgrid when HAVE_PSML is false
672  if(.false.)write(std_out,*)qgrid_vl ! Just to keep qgrid_vl when HAVE_PSML is false
673  if(.false.)write(std_out,*)useylm ! Just to keep useylm when HAVE_PSML is false
674  if(.false.)write(std_out,*)zion   ! Just to keep zion when HAVE_PSML is false
675  if(.false.)write(std_out,*)znucl  ! Just to keep znucl when HAVE_PSML is false
676 #endif
677 
678 end subroutine psp9in