Logo Search packages:      
Sourcecode: abinit version File versions  Download package

psp7in.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/psp7in
!! NAME
!! psp7in
!!
!! FUNCTION
!! Initialize pspcod=7 ("PAW pseudopotentials"):
!! continue to read the corresponding file and compute the form factors
!!
!! COPYRIGHT
!! Copyright (C) 1998-2007 ABINIT group (GJ, FJ, MT, FB)
!! 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
!!  ixc=exchange-correlation choice from main routine data file
!!  lmax=value of lmax mentioned at the second line of the psp file
!!  lmnmax=max number of (l,m,n) components over all type of psps
!!  lnmax=max number of (l,n)   components over all type of psps
!!  mmax=max number of pts in real space grid (already read in the psp file header)
!!  mqgrid_ff=dimension of q grid for array ffspl (spline fit of nl form factors)
!!  mqgrid_vl=dimension of q grid for array vlspl (spline fit of Vloc)
!!  pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1=dev. on moments)
!!  pspso= spin orbit signal
!!  qgrid_ff(mqgrid_ff)=values of q on grid (from 0 to qmax) for nl form factors
!!  qgrid_vl(mqgrid_vl)=values of q on grid (from 0 to qmax) for Vloc
!!  zion=nominal valence of atom as specified in psp file
!!
!! OUTPUT
!!  dnqdq0= 1/q d(tNcore(q))/dq for q=0 (tNcore(q) = FT of pseudo core density)
!!  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
!!  ffspl(mqgrid_ff,2,lnmax)=form factor f_l(q) and second derivative
!!   from spline fit for each angular momentum and each projector;
!!   if any, spin-orbit components begin at l=mpsang+1
!!  indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn
!!  ncspl(mqgrid_vl,2)=tncore(q) and second derivatives from spline fit
!!  pawrad <type(pawrad_type)>=paw radial mesh and related data
!!  pawtab <type(pawtab_type)>=paw tabulated starting data
!!  vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
!!  xcccrc=XC core correction cutoff radius (bohr) from psp file
!!
!! NOTES
!!  Spin-orbit not yet implemented (to be done)
!!
!! PARENTS
!!      pspatm
!!
!! CHILDREN
!!      compmesh,copymesh,deducer0,jbessel,leave_new,pawxc,pawxcm,poisson
!!      psp7cg,psp7lo,psp7nl,shapebes,simp_gen,spline,splint,wrtout
!!
!! SOURCE

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

subroutine psp7in(dnqdq0,epsatm,ffspl,indlmn,ixc,lmax,lmnmax,lnmax,mmax,&
&                  mqgrid_ff,mqgrid_vl,ncspl,pawrad,pawtab,pawxcdev,&
&                  pspso,qgrid_ff,qgrid_vl,vlspl,xcccrc,zion)

 use defs_basis
 use defs_datatypes

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_01manage_mpi
 use interfaces_11util
 use interfaces_13paw
 use interfaces_13psp, except_this_one => psp7in
#endif
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: ixc,lmax,lmnmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,pawxcdev
 integer,intent(in) :: pspso
 real(dp),intent(in) :: zion
 real(dp),intent(out) :: dnqdq0,epsatm,xcccrc
 type(pawrad_type),intent(out) :: pawrad
 type(pawtab_type),intent(out) :: pawtab
!arrays
 integer,intent(out) :: indlmn(6,lmnmax)
 real(dp),intent(in) :: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl)
 real(dp),intent(out) :: ffspl(mqgrid_ff,2,lnmax),ncspl(mqgrid_vl,2)
 real(dp),intent(out) :: vlspl(mqgrid_vl,2)

!Local variables ------------------------------
!scalars
 integer,parameter :: reduced_mshsz=2501
 integer :: creatorid,ib,icoremesh,ii,il,ilm,ilmn,iln,imainmesh,imsh,iprojmesh
 integer :: ir,iread1,iread2,ishpfmesh,isnotzero,itest,ivlocmesh,j0lmn,jlm,jlmn,jln,klmn,jj
 integer :: msz,nmesh,pspversion
 real(dp),parameter :: reduced_rstep=0.00025_dp,rm_vloc=20.0_dp
 real(dp) :: arg,argj1,argj2,intg,intvh,jbes1,jbes2,qcore,rc,rread1,rread2
 real(dp) :: shapefunc1,shapefunc2,shapefunc3,yp1,ypn
 logical :: reduced_ncor,reduced_vloc
 character(len=500) :: message,submessage
 character(len=80) :: pspline
 character :: blank=' ',numb=' '
 type(pawang_type) :: pawang_tmp
 type(pawrad_type) :: core_mesh,rcore_mesh,rvloc_mesh,shpf_mesh,tproj_mesh
 type(pawrad_type) :: vloc_mesh
!arrays
 integer :: tmp_lmselect(1)
 integer,allocatable :: nprj(:),orbitals(:)
 real(dp),allocatable :: dvlocr(:),ncore(:),nwk(:),r2k(:),rtncor(:),rvlocr(:),shpf(:,:)
 real(dp),allocatable :: tncore(:),tproj(:,:),vbare(:),vh(:),vlocr(:),work1(:)
 real(dp),allocatable :: work2(:),work3(:),work4(:)
 type(pawrad_type),allocatable :: radmesh(:)
!no_abirules
!Statement functions -----------------------------------
!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifndef HAVE_FORTRAN_INTERFACES
 integer :: ifromr
#endif
!End of the abilint section

 shapefunc1(arg)= exp(-(arg/pawtab%shape_sigma)**pawtab%shape_lambda)
 shapefunc2(arg)= (sin(pi*arg/pawtab%rshp)/(pi*arg/pawtab%rshp))**2
 shapefunc3(argj1,argj2)= pawtab%shape_alpha(1,1)*argj1+pawtab%shape_alpha(2,1)*argj2

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

!  File format of formatted PAW psp input (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
!  (4) psp_version, creatorID
!  (5) basis_size, lmn_size
!  (6) orbitals (for l=1 to basis_size)
!  (7) number_of_meshes
!  For imsh=1 to number_of_meshes
!      (8)  mesh_index, mesh_type ,mesh_size, rad_step[, log_step]
!  (9) r_cut(SPH)
!  (10) shape_type, r_shape[, shapefunction arguments]
!  For iln=1 to basis_size
!      (11) comment(character)
!      (12) radial mesh index for phi
!      (13) phi(r) (for ir=1 to phi_meshsz)
!  For iln=1 to basis_size
!      (14) comment(character)
!      (15) radial mesh index for tphi
!      (16) tphi(r) (for ir=1 to phi_mesh_size)
!  For iln=1 to basis_size
!      (17) comment(character)
!      (18) radial mesh index for tproj
!      (19) tproj(r) (for ir=1 to proj_mesh_size)
!  (20) comment(character)
!  (21) radial mesh index for core_density
!  (22) core_density (for ir=1 to phi_mesh_size)
!  (23) comment(character)
!  (24) radial mesh index for tcore_density
!  (25) tcore_density (for ir=1 to phi_mesh_size)
!  (26) comment(character)
!  (27) Dij0 (for ij=1 to lmn_size*(lmn_size+1)/2)
!  (28) comment(character)
!  (29) Rhoij0 (for ij=1 to lmn_size*(lmn_size+1)/2)
!  (30) comment(character)
!  (31) radial mesh index for Vloc, format of Vloc (0=Vbare, 1=VH(tnzc))
!  (32) Vloc(r) (for ir=1 to vloc_mesh_size)
!  ===== Following lines only if shape_type=-1 =====
!  For il=1 to 2*max(orbitals)+1
!      (33) comment(character)
!      (34) radial mesh index for shapefunc
!      (35) shapefunc(r)*gnorm(l)*r**l (for ir=1 to phi_meshsz)
!
!  Comments:
!  * psp_version= ID of PAW_psp version
!     4 characters string of the form 'pawn' (with n varying)
!  * creatorID= ID of psp generator
!     creatorid=1xyz : psp generated from Holzwarth AtomPAW generator version x.yz
!     creatorid=2xyz : psp generated from Vanderbilt ultra-soft generator version x.yz
!     creatorid=-1: psp for tests (for developpers only)
!  * mesh_type= type of radial mesh
!     mesh_type=1 (regular grid): rad(i)=(i-1)*AA
!     mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
!     mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
!     mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
!  * radial shapefunction type
!     shape_type=-1 ; g(r)=numeric (read from psp file)
!     shape_type= 1 ; g(r)=exp[-(r/sigma)**lambda]
!     shape_type= 2 ; g(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2 if r<=rshp, zero if r>rshp
!     shape_type= 3 ; gl(r)=Alpha(1,l)*jl(q(1,l)*r)+Alpha(2,l)*jl(q(2,l)*r) for each l (not implemented)


!==========================================================

!DEBUG
!write(6,*)' psp7in : enter'
!ENDDEBUG

!==========================================================
!Read psp version in line 4 of the header

 pspversion=1
 read (tmp_unit,'(a80)') pspline;pspline=adjustl(pspline)
 if (pspline(1:3)=="paw".or.pspline(1:3)=="PAW") &
& read(unit=pspline(4:80),fmt=*) pspversion
 if (pspversion/=1.and.pspversion/=2.and.pspversion/=3) then
  write(message, '(a,a,a,a,i2,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  This version of PAW psp file (',pspversion,') is not compatible with',ch10,&
&   '  current version of Abinit.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!==========================================================
!Read lines 4 to 11 of the header

!Have to maintain compatibility with Abinit v4.2.x
 if (pspversion==1) then
  read (unit=pspline,fmt=*) pawtab%basis_size,pawtab%lmn_size
  allocate(orbitals(pawtab%basis_size))
  read(tmp_unit,*) (orbitals(ib), ib=1,pawtab%basis_size)
  pawtab%l_size=2*maxval(orbitals)+1
  nmesh=3;allocate(radmesh(nmesh))
  read(tmp_unit,'(a80)') pspline
  radmesh(1)%lstep=zero
  read(unit=pspline,fmt=*,err=10,end=10) radmesh(1)%mesh_type,&
&                  radmesh(1)%rstep,radmesh(1)%lstep
10 read(tmp_unit,*) pawtab%rpaw
  read(tmp_unit,*) radmesh(1)%mesh_size,radmesh(2)%mesh_size,&
&                  radmesh(3)%mesh_size
  read(tmp_unit,'(a80)') pspline
  pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99
  read(unit=pspline,fmt=*,err=11,end=11) pawtab%shape_type,&
&                          pawtab%shape_lambda,pawtab%shape_sigma
11 read(tmp_unit,*) creatorid
  if (pawtab%shape_type==3) pawtab%shape_type=-1
  radmesh(2)%mesh_type=radmesh(1)%mesh_type
  radmesh(3)%mesh_type=radmesh(1)%mesh_type
  radmesh(2)%rstep=radmesh(1)%rstep
  radmesh(3)%rstep=radmesh(1)%rstep
  radmesh(2)%lstep=radmesh(1)%lstep
  radmesh(3)%lstep=radmesh(1)%lstep
 else

!Here psp file for Abinit 4.3+
  read (unit=pspline(5:80),fmt=*) creatorid
  read (tmp_unit,*) pawtab%basis_size,pawtab%lmn_size
  allocate(orbitals(pawtab%basis_size))
  read(tmp_unit,*) (orbitals(ib), ib=1,pawtab%basis_size)
  pawtab%l_size=2*maxval(orbitals)+1
  read(tmp_unit,*) nmesh
  allocate(radmesh(nmesh))
  do imsh=1,nmesh
   rread2=zero
   read(tmp_unit,'(a80)') pspline
   read(unit=pspline,fmt=*,err=20,end=20) ii,iread1,iread2,rread1,rread2
20 continue
   if (ii<=nmesh) then
    radmesh(ii)%mesh_type=iread1
    radmesh(ii)%mesh_size=iread2
    radmesh(ii)%rstep=rread1
    radmesh(ii)%lstep=rread2
   else
    write(message, '(a,a,a,a,a,a)' ) ch10,&
&     ' psp7in: ERROR -',ch10,&
&     '  Index of mesh out of range !',ch10,&
&     '  Action : check your pseudopotential file.'
    call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
    call leave_new('COLL')
   end if
  end do
  read(tmp_unit,*) pawtab%rpaw
  read(tmp_unit,'(a80)') pspline
  read(unit=pspline,fmt=*) pawtab%shape_type
  pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99
 end if

!Here reading shapefunction parameters for Abinit 4.3...4.5
 if (pspversion==2) then
  if (pawtab%shape_type==1) read(unit=pspline,fmt=*) ii,pawtab%shape_lambda,pawtab%shape_sigma
  if (pawtab%shape_type==3) pawtab%shape_type=-1
  pawtab%rshp=zero

!Here reading shapefunction parameters for Abinit 4.6+
 else if (pspversion>=3) then
  pawtab%rshp=zero
  if (pawtab%shape_type==-1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp
  if (pawtab%shape_type== 1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp, &
&                            pawtab%shape_lambda,pawtab%shape_sigma
  if (pawtab%shape_type== 2) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp
  if (pawtab%shape_type== 3) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp
 end if
21 continue

!If shapefunction type is Bessel, deduce here its parameters from rc
 if (pawtab%shape_type== 3) then
  rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw
  do il=1,pawtab%l_size
   call shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc)
  end do
 end if
!==========================================================
!Mirror pseudopotential parameters to the output and log files

 write(message,'(a,i1)')' Pseudopotential format is: paw',pspversion
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 write(message,'(2(a,i3),a,64i4)') &
&     ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',&
&     pawtab%lmn_size,'), orbitals=',orbitals(1:pawtab%basis_size)
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 write(message,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 write(message,'(a,i1,a)')' ',nmesh,' radial meshes are used:'
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 do imsh=1,nmesh
  if (radmesh(imsh)%mesh_type==1) &
&  write(message,'(a,i1,a,i4,a,g12.5)') &
&   '  - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,&
&   ' , step=',radmesh(imsh)%rstep
  if (radmesh(imsh)%mesh_type==2) &
&  write(message,'(a,i1,a,i4,2(a,g12.5))') &
&   '  - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,&
&   ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep
  if (radmesh(imsh)%mesh_type==3) &
&  write(message,'(a,i1,a,i4,2(a,g12.5))') &
&   '  - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,&
&   ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep
  if (radmesh(imsh)%mesh_type==4) &
&  write(message,'(a,i1,a,i4,a,g12.5)') &
&   '  - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,&
&   ' , AA=',radmesh(imsh)%rstep
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 end do
 if (pawtab%shape_type==-1) then
  write(message,'(a)')&
      ' Shapefunction is NUMERIC type: directly read from atomic data file'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 end if
 if (pawtab%shape_type==1) then
  write(message,'(2a,a,f6.3,a,i3)')&
&     ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,&
&     '                            with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 end if
 if (pawtab%shape_type==2) then
  write(message,'(a)')&
      ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 end if
 if (pawtab%shape_type==3) then
  write(message,'(a)')&
&     ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 end if
 if (pawtab%rshp<1.d-8) then
  write(message,'(a)') ' Radius for shape functions = sphere core radius'
 else
  write(message,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp
 end if
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')

!==========================================================
!Perfom tests

!Spin-orbit not yet implemented
 if (pspso>1) then
  write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Spin-orbit not yet implemented for PAW psps !',ch10,&
&   '  Action : check your pseudopotential or input file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!Are lmax and orbitals compatibles ?
 if (lmax/=maxval(orbitals)) then
  write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  lmax /= MAX(orbitals) !',ch10,&
&   '  Action: check your pseudopotential file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!Are lmnmax and pawtab%lmn_size compatibles ?
 if (lmnmax<pawtab%lmn_size) then
  write(message, '(a,a,a,a)' ) ch10,&
&   ' psp7in: BUG -',ch10,&
&   '  lmnmax < lmn size (read from pseudo) !'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!Only mesh_type=1,2, 3 or 4 allowed
 do imsh=1,nmesh
  if (radmesh(imsh)%mesh_type>4) then
   write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Only mesh types 1,2, 3 or 4 allowed !',ch10,&
&   '  Action : check your pseudopotential or input file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
  end if
 end do

!==========================================================
!Initialize radial meshes

 do imsh=1,nmesh
  allocate(radmesh(imsh)%rad (radmesh(imsh)%mesh_size))
  allocate(radmesh(imsh)%radfact(radmesh(imsh)%mesh_size))
  allocate(radmesh(imsh)%simfact(radmesh(imsh)%mesh_size))
  call compmesh(radmesh(imsh),-1._dp)
 end do

!==========================================================
!Read tabulated atomic data

!---------------------------------
!Read wave-functions (phi)
 do ib=1,pawtab%basis_size
  read (tmp_unit,*)
  if (pspversion==1) iread1=1
  if (pspversion>1) read (tmp_unit,*) iread1
  if (ib==1) then
   pawrad%mesh_type=radmesh(iread1)%mesh_type
   pawrad%mesh_size=radmesh(iread1)%mesh_size
   pawrad%rstep=radmesh(iread1)%rstep
   pawrad%lstep=radmesh(iread1)%lstep
   call compmesh(pawrad,pawtab%rpaw)
   imainmesh=iread1
  else if (iread1/=imainmesh) then
   write(message, '(a,a,a,a,a,a)' ) ch10,&
&    ' psp7in: ERROR -',ch10,&
&    '  All Phi and tPhi must be given on the same radial mesh !',ch10,&
&    '  Action: check your pseudopotential file.'
   call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
   call leave_new('COLL')
  end if
  read (tmp_unit,*) (pawtab%phi(ir,ib),ir=1,pawrad%mesh_size)
 end do

!---------------------------------
!Read pseudo wave-functions (tphi)
 do ib=1,pawtab%basis_size
  read (tmp_unit,*)
  if (pspversion==1) iread1=1
  if (pspversion>1) read (tmp_unit,*) iread1
  if (iread1/=imainmesh) then
   write(message, '(a,a,a,a,a,a)' ) ch10,&
&    ' psp7in: ERROR -',ch10,&
&    '  All Phi and tPhi must be given on the same radial mesh !',ch10,&
&    '  Action: check your pseudopotential file.'
   call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
   call leave_new('COLL')
  end if
  read (tmp_unit,*) (pawtab%tphi(ir,ib),ir=1,pawrad%mesh_size)
 end do
 write(message,'(a,i1)') &
& ' Radial grid used for partial waves is grid ',imainmesh
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')

!---------------------------------
!Read projectors (tproj)
 do ib=1,pawtab%basis_size
  read (tmp_unit,*)
  if (pspversion==1) iread1=2
  if (pspversion>1) read (tmp_unit,*) iread1
  if (ib==1) then
   call copymesh(radmesh(iread1),tproj_mesh)
   iprojmesh=iread1
   allocate(tproj(tproj_mesh%mesh_size,pawtab%basis_size))
  else if (iread1/=iprojmesh) then
   write(message, '(a,a,a,a,a,a)' ) ch10,&
&    ' psp7in: ERROR -',ch10,&
&    '  All tprojectors must be given on the same radial mesh !',ch10,&
&    '  Action: check your pseudopotential file.'
   call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
   call leave_new('COLL')
  end if
  read (tmp_unit,*) (tproj(ir,ib),ir=1,tproj_mesh%mesh_size)
 end do
 write(message,'(a,i1)') &
& ' Radial grid used for projectors is grid ',iprojmesh
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')

!---------------------------------
!Read core density (coredens)
 read (tmp_unit,*)
 if (pspversion==1) iread1=1
 if (pspversion>1) read (tmp_unit,*) iread1
 icoremesh=iread1;call copymesh(radmesh(iread1),core_mesh)
 if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.&
&    (radmesh(icoremesh)%rstep    /=pawrad%rstep)    .or.&
&    (radmesh(icoremesh)%lstep    /=pawrad%lstep)) then
  write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Ncore must be given on a radial mesh with the same',ch10,&
&   '  type and step(s) than the main radial mesh (mesh for Phi) !',ch10,&
&   '  Action: check your pseudopotential file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if
 allocate(ncore(core_mesh%mesh_size))
 read (tmp_unit,*) (ncore(ir),ir=1,core_mesh%mesh_size)
 pawtab%coredens(1:pawrad%mesh_size)=ncore(1:pawrad%mesh_size)

!---------------------------------
!Read tcore density (tcoredens)
 read (tmp_unit,*)
 if (pspversion==1) iread1=1
 if (pspversion>1) read (tmp_unit,*) iread1
 if (iread1/=icoremesh) then
  write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Pseudized core density (tNcore) must be given',ch10,&
&   '  on the same radial mesh as core density (Ncore) !',ch10,&
&   '  Action: check your pseudopotential file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if
 allocate(tncore(core_mesh%mesh_size))
 read (tmp_unit,*) (tncore(ir),ir=1,core_mesh%mesh_size)
 if (maxval(abs(tncore(:)))<tol6) then
  pawtab%usetcore=0
  pawtab%tcoredens(1:pawrad%mesh_size)=zero
 else
  pawtab%usetcore=1
  pawtab%tcoredens(1:pawrad%mesh_size)=tncore(1:pawrad%mesh_size)
 end if
 write(message,'(a,i1)') &
& ' Radial grid used for (t)core density is grid ',icoremesh
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
!---------------------------------
!Read frozen part of Dij terms (dij0)
 read (tmp_unit,*)
 read (tmp_unit,*) (pawtab%dij0(ib),ib=1,pawtab%lmn2_size)

!---------------------------------
!Read initial guess of rhoij (rhoij0)
 read (tmp_unit,*)
 read (tmp_unit,*) (pawtab%rhoij0(ib),ib=1,pawtab%lmn2_size)

!---------------------------------
!Read local pseudopotential=Vh(tn_zc) or Vbare
 read (tmp_unit,*)
 if (pspversion==1) ivlocmesh=3
 pawtab%vlocopt=1
 if (pspversion==2) then
  read (tmp_unit,*) ivlocmesh
 else if (pspversion>2) then
! read (tmp_unit,fmt=*,err=30,end=30) ivlocmesh,pawtab%vlocopt
  message=blank
  read (tmp_unit,fmt=*) message
  read (message,fmt=*) ivlocmesh
  write(numb,'(i1)')ivlocmesh
  ii=index(message,numb)
  submessage=trim(message(ii+1:))
  if(len_trim(submessage)/=0)then
   do ii=1,len_trim(submessage) 
    numb=submessage(ii:ii)
    if(numb==blank)cycle
    jj=index('123456789',numb)
    if(jj<1 .or. jj>9)exit
    pawtab%vlocopt=jj
   enddo
  endif
 end if
 call copymesh(radmesh(ivlocmesh),vloc_mesh)
 allocate(vlocr(vloc_mesh%mesh_size))
 read (tmp_unit,*) (vlocr(ir),ir=1,vloc_mesh%mesh_size)
 write(message,'(a,i1)') &
& ' Radial grid used for Vloc is grid ',ivlocmesh
 call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
!---------------------------------
!Eventually read "numeric" shapefunctions (if shape_type=-1)
 if (pawtab%shape_type==-1) then
  do il=1,pawtab%l_size
   read (tmp_unit,*)
   if (pspversion==1) iread1=1
   if (pspversion>1) read (tmp_unit,*) iread1
   if (il==1) then
    call copymesh(radmesh(iread1),shpf_mesh)
    ishpfmesh=iread1
    allocate(shpf(shpf_mesh%mesh_size,pawtab%l_size))
   else if (iread1/=ishpfmesh) then
    write(message, '(a,a,a,a,a,a)' ) ch10,&
&     ' psp7in: ERROR -',ch10,&
&     '  All shape functions must be given on the same radial mesh !',ch10,&
&     '  Action: check your pseudopotential file.'
    call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
    call leave_new('COLL')
   end if
   read (tmp_unit,*) (shpf(ir,il),ir=1,shpf_mesh%mesh_size)
  end do
  write(message,'(a,i1)') &
&  ' Radial grid used for shape functions is grid ',iread1
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')

! Has to spline shape functions if mesh is not the "main" mesh
  if (ishpfmesh/=imainmesh) then
   msz=shpf_mesh%mesh_size
   allocate(work1(msz),work2(msz),work3(msz),work4(pawrad%mesh_size))
   work3(:)=shpf_mesh%rad(:)
   work4(:)=pawrad%rad(:)
   do il=1,pawtab%l_size
    if (shpf_mesh%radfact(1)>zero) then
     yp1=1._dp/12._dp/shpf_mesh%stepint/shpf_mesh%radfact(1) &
&     *(-25._dp*shpf(1,il)+48._dp*shpf(2,il)-36._dp*shpf(3,il)+16._dp*shpf(4,il)-3._dp*shpf(5,il))
    else
     yp1=(shpf(2,il)-shpf(1,il))/(shpf_mesh%rad(2)-shpf_mesh%rad(1))
    end if
    ypn=1._dp/12._dp/shpf_mesh%stepint &
&    *(  3._dp*shpf(msz-4,il)-16._dp*shpf(msz-3,il) &
&      +36._dp*shpf(msz-2,il)-48._dp*shpf(msz-1,il) &
&      +25._dp*shpf(msz  ,il)) / shpf_mesh%radfact(msz)
    call spline(work3,shpf(:,il),msz,yp1,ypn,work1,work2)
    call splint(msz,work3,shpf(:,il),work1,pawrad%mesh_size,work4,pawtab%shapefunc(:,il))
   end do
   deallocate(work1,work2,work3,work4)
  else
   pawtab%shapefunc(:,:)=shpf(:,:)
  end if
  deallocate(shpf)
 end if
!==========================================================
! Perfom tests on meshes

!Are radial meshes for Phi and Vloc compatibles ?
 if (vloc_mesh%rmax<pawrad%rmax) then
  write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Rmax for Vloc < Rmax for Phi !',ch10,&
&   '  Action : check your pseudopotential (increase Vloc meshSize).'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!Are mmax and mesh_size for partial waves compatibles ?
 if (mmax/=pawrad%mesh_size) then
  write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  mmax /= phi_mesh_size in psp file !',ch10,&
&   '  Action: check your pseudopotential file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if
!Are radial meshes for (t)Ncore and Phi compatibles ?
 if (core_mesh%mesh_size<pawrad%mesh_size) then
  write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Mesh size for core density must be equal or larger',ch10,&
&   '  than mesh size for spheres (partial waves) !',ch10,&
&   '  Action : check your pseudopotential (increase Ncore meshSize).'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!Is PAW radius included inside radial mesh ?
 if (pawtab%rpaw>pawrad%rmax+tol8) then
  write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Radius of PAW sphere is outside the radial mesh !',ch10,&
&   '  Action: check your pseudopotential file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!Max. radius of mesh for Vloc has to be "small" in order to avoid numeric noise ?
 if (vloc_mesh%rmax>rm_vloc) then
  vloc_mesh%mesh_size=ifromr(vloc_mesh,rm_vloc)
  vloc_mesh%rmax=vloc_mesh%rad(vloc_mesh%mesh_size)
  write(message, '(a,a,a,a,f6.2,a,a,a,a,a,i4,a)' ) ch10,&
&   ' psp7in: WARNING -',ch10,&
&   '  Max. radius for Vloc was too large (>',rm_vloc,' a.u.) !',ch10,&
&   '  Numeric noise was possible.',ch10,&
&   '  Mesh size for Vloc has been set to ',vloc_mesh%mesh_size,'.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
 end if
!This test has been disable... MT 2006-25-10
!For Simpson rule, it is better to have odd mesh sizes
! itest=0
! do imsh=1,nmesh
!  if (mod(radmesh(imsh)%mesh_size,2)==0.and.radmesh(imsh)%mesh_type==1) itest=1
! end do
! if (itest==1) then
!  write(message, '(8a)' ) ch10,&
!&   ' psp7in: WARNING -',ch10,&
!&   '  Regular radial meshes should have odd number of points ',ch10,&
!&   '  for better accuracy of integration sheme (Simpson rule).',ch10,&
!&   '  Althought it''s not compulsory, you should change mesh sizes in psp file.'
!  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
! end if

!Test the compatibilty between Rpaw and mesh for (t)Phi
 if (pspversion>=3) then
  itest=ifromr(radmesh(imainmesh),pawtab%rpaw)
  if (itest+2>radmesh(imainmesh)%mesh_size) then
   write(message, '(12a)' ) ch10,&
&   ' psp7in: WARNING -',ch10,&
&   '  Atomic data could produce inaccurate results:',ch10,&
&   '    Wavefunctions and pseudo-wavefunctions should',ch10,&
&   '    be given on a radial mesh larger than the PAW',ch10,&
&   '    spheres (at least 2 additional points) !',ch10,&
&   '  Action: check your pseudopotential file.'
   call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  end if
  if (abs(pawtab%rpaw-radmesh(imainmesh)%rad(itest))<tol8) itest=itest-1
  ib=0;isnotzero=0
  do while ((isnotzero==0).and.(ib<pawtab%basis_size))
   ib=ib+1;ir=itest
   do while ((isnotzero==0).and.(ir<radmesh(imainmesh)%mesh_size))
    ir=ir+1;if (abs(pawtab%phi(ir,ib)-pawtab%tphi(ir,ib))>tol8) isnotzero=1
   end do
  end do
  if (isnotzero>0) then
   write(message, '(10a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Atomic data are inconsistent:',ch10,&
&   '  For r>=r_paw, pseudo wavefunctions are not',ch10,&
&   '  equal to wave functions (Phi(r)/=tPhi(r)) !',ch10,&
&   '  Action: check your pseudopotential file.'
   call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
   call leave_new('COLL')
  end if
 else
! For compatibility reasons set PAW radius at the end of mesh (older versions)
  if (pawtab%rpaw/=pawrad%rmax) then
   pawtab%rpaw=pawrad%rmax
   call compmesh(pawrad,-1._dp)
  end if
 end if
!DEBUG
 write(6,*)' pawtab%vlocopt=',pawtab%vlocopt
 write(6,*)' vloc_mesh%rmax=',vloc_mesh%rmax
 write(6,*)' pawtab%rpaw=',pawtab%rpaw
!ENDDEBUG
!If Vloc is a "Vbare" potential, it has to be localized inside PAW spheres
 if (pawtab%vlocopt==0.and.(vloc_mesh%rmax>pawtab%rpaw+tol10)) then
  write(message, '(8a)' ) ch10,&
&  ' psp7in: ERROR -',ch10,&
&  '  Atomic data are inconsistent:',ch10,&
&  '  Local potential is a "Vbare" potential',ch10,&
&  '  and is not localized inside PAW sphere !'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if

!==========================================================
!Initialize various parameters in structured data

!Some dims
 pawtab%l_size=2*maxval(orbitals)+1
 pawtab%mesh_size=pawrad%mesh_size
 pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2

!indlmn calculation (indices for (l,m,n) basis)
 allocate(nprj(0:maxval(orbitals)))
 ilmn=0;iln=0;nprj=0
 do ib=1,pawtab%basis_size
  il=orbitals(ib)
  nprj(il)=nprj(il)+1
  iln=iln+1
  do ilm=1,2*il+1
   indlmn(1,ilmn+ilm)=il
   indlmn(2,ilmn+ilm)=ilm-(il+1)
   indlmn(3,ilmn+ilm)=nprj(il)
   indlmn(4,ilmn+ilm)=il*il+ilm
   indlmn(5,ilmn+ilm)=iln
   indlmn(6,ilmn+ilm)=1
  end do
  ilmn=ilmn+2*il+1
 end do
 deallocate(nprj,orbitals)

!Are ilmn (found here) and pawtab%lmn_size compatibles ?
 if (ilmn/=pawtab%lmn_size) then
  write(message, '(a,a,a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Calculated lmn size differs from',ch10,&
&   '  lmn_size read from pseudo !',ch10,&
&   ' Action: check your pseudopotential file.'
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
  call leave_new('COLL')
 end if
!==========================================================
!Compute ffspl(q) (and derivatives)

 ffspl=zero
 call psp7nl(ffspl,indlmn,pawtab%lmn_size,lnmax,mqgrid_ff,qgrid_ff,tproj_mesh,tproj)
 deallocate(tproj)
!==========================================================
!Compute eventually compensation charge radius (i.e. radius for shape functions)

 if (pawtab%shape_type>0.and.pawtab%rshp<1.d-8) then
  pawtab%rshp=pawtab%rpaw
 else if (pawtab%shape_type==-1) then
  ir=ifromr(radmesh(imainmesh),pawtab%rpaw)+1;isnotzero=0
  do while ((isnotzero==0).and.(ir>1))
   ir=ir-1;il=0
   do while ((isnotzero==0).and.(il<pawtab%l_size))
    il=il+1;if (pawtab%shapefunc(ir,il)>tol16) isnotzero=1
   end do
  end do
  ir=min(ir+1,radmesh(imainmesh)%mesh_size)
  pawtab%rshp=radmesh(imainmesh)%rad(ir)
  do il=1,pawtab%l_size
   if (pawtab%shapefunc(ir,il)>tol6) then
    write(message, '(a,a,a,a,a,a)' ) ch10,&
&    ' psp7in: ERROR -',ch10,&
&    '  Shape function is not zero at PAW radius !',ch10,&
&    '  Action: check your pseudopotential file.'
    call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
    call leave_new('COLL')
   end if
  end do
 end if
!==========================================================
!If read Vloc potential is in "Vbare" format,
!   translate it into VH(tnzc) format

 if (pawtab%vlocopt==0) then

  write(message,'(a)') ' Local potential is in "Vbare" format... '
  call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')

  allocate(vh(core_mesh%mesh_size),vbare(core_mesh%mesh_size))
  allocate(r2k(radmesh(imainmesh)%mesh_size))

! Compute r2.k(r)
  if (pawtab%shape_type==-1) then
   do ir=1,radmesh(imainmesh)%mesh_size
    r2k(ir)=pawtab%shapefunc(ir,1)*core_mesh%rad(ir)**2
   end do
  else if (pawtab%shape_type==1) then
   do ir=1,radmesh(imainmesh)%mesh_size
    r2k(ir)=shapefunc1(core_mesh%rad(ir))*core_mesh%rad(ir)**2
   end do
  else if (pawtab%shape_type==2) then
   r2k(1)=zero
   do ir=2,radmesh(imainmesh)%mesh_size
    r2k(ir)=shapefunc2(core_mesh%rad(ir))*core_mesh%rad(ir)**2
   end do
  else if (pawtab%shape_type==3) then
   do ir=1,radmesh(imainmesh)%mesh_size
    CALL jbessel(jbes1,yp1,ypn,0,0,pawtab%shape_q(1,1)*core_mesh%rad(ir))
    CALL jbessel(jbes2,yp1,ypn,0,0,pawtab%shape_q(2,1)*core_mesh%rad(ir))
    r2k(ir)=shapefunc3(jbes1,jbes2)*core_mesh%rad(ir)**2
   end do
  end if
  call simp_gen(intg,r2k,radmesh(imainmesh))
  intg=1._dp/intg;r2k(:)=r2k(:)*intg
! Compute VH[4pi.r2.n(r)=4pi.r2.tncore(r)+(Qcore-Z).r2.k(r)]
  allocate(nwk(core_mesh%mesh_size))
  nwk(:)=tncore(:)*four_pi*core_mesh%rad(:)**2
  call simp_gen(qcore,nwk,core_mesh)
  nwk(1:radmesh(imainmesh)%mesh_size)=nwk(1:radmesh(imainmesh)%mesh_size) &
&                                    -r2k(1:radmesh(imainmesh)%mesh_size)*(qcore+zion)
  call poisson(nwk,0,intg,core_mesh,vh)
  vh(2:core_mesh%mesh_size)=half*vh(2:core_mesh%mesh_size)/core_mesh%rad(2:core_mesh%mesh_size)
  call deducer0(vh,core_mesh%mesh_size,core_mesh)
  deallocate(nwk)

!Eventually spline Vbare
  if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.&
&     (core_mesh%rstep    /=vloc_mesh%rstep)    .or.&
&     (core_mesh%lstep    /=vloc_mesh%lstep)) then
   msz=core_mesh%mesh_size
   if (vloc_mesh%rmax<core_mesh%rmax) msz=ifromr(core_mesh,vloc_mesh%rmax)
   if (vloc_mesh%radfact(1)>zero) then
    yp1=1._dp/12._dp/vloc_mesh%stepint/vloc_mesh%radfact(1) &
&       *(-25._dp*vlocr(1)+48._dp*vlocr(2)-36._dp*vlocr(3)+16._dp*vlocr(4)-3._dp*vlocr(5))
   else
    yp1=(vlocr(2)-vlocr(1))/(vloc_mesh%rad(2)-vloc_mesh%rad(1))
   end if
   ypn=1._dp/12._dp/vloc_mesh%stepint &
&      *(  3._dp*vlocr(vloc_mesh%mesh_size-4)-16._dp*vlocr(vloc_mesh%mesh_size-3) &
&        +36._dp*vlocr(vloc_mesh%mesh_size-2)-48._dp*vlocr(vloc_mesh%mesh_size-1) &
&        +25._dp*vlocr(vloc_mesh%mesh_size)) / vloc_mesh%radfact(vloc_mesh%mesh_size)
   allocate(work1(vloc_mesh%mesh_size),work2(vloc_mesh%mesh_size),work3(vloc_mesh%mesh_size))
   work3(:)=vloc_mesh%rad(:)
   call spline(work3,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1,work2)
   call splint(vloc_mesh%mesh_size,work3,vlocr,work1,msz,core_mesh%rad(1:msz),vbare)
   deallocate(work1,work2,work3)
  else
   msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size)
   vbare(1:msz)=vlocr(1:msz)
  end if
! Build VH(tnzc) from Vbare
  deallocate(vlocr,vloc_mesh%rad,vloc_mesh%radfact,vloc_mesh%simfact)
  call copymesh(core_mesh,vloc_mesh)
  allocate(vlocr(core_mesh%mesh_size))
  vlocr(:)=vh(:)

! Compute <tPhi_i|VH(tnzc)|tPhi_j> part of Dij0
! Compute int[VH(tnzc)*Qijhat(r)dr] part of Dij0
  msz=radmesh(imainmesh)%mesh_size;allocate(work1(msz))
  work1(1:msz)=vlocr(1:msz)*r2k(1:msz)
  call simp_gen(intvh,work1,radmesh(imainmesh))
  do jlmn=1,pawtab%lmn_size
   j0lmn=jlmn*(jlmn-1)/2;jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
   do ilmn=1,jlmn
    klmn=j0lmn+ilmn;ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
    if (jlm==ilm) then
     work1(1:msz)=pawtab%tphi(1:msz,iln)*pawtab%tphi(1:msz,jln)*(vlocr(1:msz)-intvh) &
&                -pawtab%phi (1:msz,iln)*pawtab%phi (1:msz,jln)*intvh
     call simp_gen(intg,work1,radmesh(imainmesh))
     pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
    end if
   end do
  end do
  deallocate(work1)
  deallocate(r2k)

 end if
!==========================================================
!Try to optimize CPU time:
!If Vloc mesh size is big, spline Vloc into a smaller log. mesh

 reduced_vloc=(vloc_mesh%mesh_size>int(reduced_mshsz))
 if (reduced_vloc) then
  msz=vloc_mesh%mesh_size
  rvloc_mesh%mesh_type=3
  rvloc_mesh%mesh_size=reduced_mshsz
  rvloc_mesh%rstep    =reduced_rstep
  rvloc_mesh%lstep    =log(0.9999999_dp*vloc_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2)
  allocate(rvloc_mesh%rad (reduced_mshsz))
  allocate(rvloc_mesh%radfact(reduced_mshsz))
  allocate(rvloc_mesh%simfact(reduced_mshsz))
  allocate(rvlocr(reduced_mshsz))
  call compmesh(rvloc_mesh,-1._dp)
  if (vloc_mesh%radfact(1)>zero) then
   yp1=1._dp/12._dp/vloc_mesh%stepint/vloc_mesh%radfact(1) &
&   *(-25._dp*vlocr(1)+48._dp*vlocr(2)-36._dp*vlocr(3)+16._dp*vlocr(4)-3._dp*vlocr(5))
  else
   yp1=(vlocr(2)-vlocr(1))/(vloc_mesh%rad(2)-vloc_mesh%rad(1))
  end if
  ypn=1._dp/12._dp/vloc_mesh%stepint &
&      *(  3._dp*vlocr(msz-4)-16._dp*vlocr(msz-3) &
&        +36._dp*vlocr(msz-2)-48._dp*vlocr(msz-1) &
&        +25._dp*vlocr(msz)) / vloc_mesh%radfact(msz)
  allocate(work1(msz),work2(msz),work3(msz))
  work3(:)=vloc_mesh%rad(:)
  call spline(work3,vlocr,msz,yp1,ypn,work1,work2)
  call splint(msz,work3,vlocr,work1,reduced_mshsz,rvloc_mesh%rad,rvlocr)
  deallocate(work1,work2,work3)
 end if
!==========================================================
!Try to optimize CPU time:
!If ncore mesh size is big, spline tncore into a smaller log. mesh

 reduced_ncor=(core_mesh%mesh_size>int(reduced_mshsz)).and.(pawtab%usetcore/=0)
 if (reduced_ncor) then
  msz=core_mesh%mesh_size
  rcore_mesh%mesh_type=3
  rcore_mesh%mesh_size=reduced_mshsz
  rcore_mesh%rstep    =reduced_rstep
  rcore_mesh%lstep    =log(0.9999999_dp*core_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2)
  allocate(rcore_mesh%rad (reduced_mshsz))
  allocate(rcore_mesh%radfact(reduced_mshsz))
  allocate(rcore_mesh%simfact(reduced_mshsz))
  allocate(rtncor(reduced_mshsz))
  call compmesh(rcore_mesh,-1._dp)
  if (core_mesh%radfact(1)>zero) then
   yp1=1._dp/12._dp/core_mesh%stepint/core_mesh%radfact(1) &
&   *(-25._dp*tncore(1)+48._dp*tncore(2)-36._dp*tncore(3)+16._dp*tncore(4)-3._dp*tncore(5))
  else
   yp1=(tncore(2)-tncore(1))/(core_mesh%rad(2)-core_mesh%rad(1))
  end if
  ypn=1._dp/12._dp/core_mesh%stepint &
&      *(  3._dp*tncore(msz-4)-16._dp*tncore(msz-3) &
&        +36._dp*tncore(msz-2)-48._dp*tncore(msz-1) &
&        +25._dp*tncore(msz)) / core_mesh%radfact(msz)
  allocate(work1(msz),work2(msz),work3(msz))
  work3(:)=core_mesh%rad(:)
  call spline(work3,tncore,msz,yp1,ypn,work1,work2)
  call splint(msz,work3,tncore,work1,reduced_mshsz,rcore_mesh%rad,rtncor)
  deallocate(work1,work2,work3)
 end if
!==========================================================
!Compute Vlspl(q) (and second derivative) from Vloc(r)

!Compute Vlspl(q)=-Zv/pi + q^2.Vloc(q) from vloc(r)
 if (reduced_vloc) then
  call psp7lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),rvloc_mesh,rvlocr,yp1,ypn,zion)
 else
  call psp7lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),vloc_mesh,vlocr,yp1,ypn,zion)
 end if

!Compute second derivative of Vlspl(q)
 allocate(work1(mqgrid_vl))
 call spline(qgrid_vl,vlspl(:,1),mqgrid_vl,yp1,ypn,vlspl(:,2),work1)
 deallocate(work1)

!==========================================================
!Compute Ncspl(q) (and second derivative) from tNcore(r)

 if (pawtab%usetcore/=0) then
! Compute Ncspl(q)=Nc(q) from Ncore(r)
  if (reduced_ncor) then
   call psp7cg(dnqdq0,mqgrid_vl,qgrid_vl,ncspl(:,1),rcore_mesh,rtncor,yp1,ypn)
  else
   call psp7cg(dnqdq0,mqgrid_vl,qgrid_vl,ncspl(:,1),core_mesh,tncore,yp1,ypn)
  end if
! Compute second derivative of Ncspl(q)
  allocate(work1(mqgrid_vl))
  call spline(qgrid_vl,ncspl(:,1),mqgrid_vl,yp1,ypn,ncspl(:,2),work1)
  deallocate (work1)
 else
  ncspl=zero
  dnqdq0=zero
 end if
!OBSOLETE: call psp7cc(core_mesh,pawtab%usetcore,xcccrc,tncore,xccc1d)
 xcccrc=core_mesh%rmax

!==========================================================
!Compute first  and third derivative of Vloc(r) (only if needed)
 if (pawtab%usedvloc==1) then
! --- First derivative
  if (ivlocmesh==imainmesh) then
   call nderiv_gen(pawtab%dvlocdr(:,1),vlocr,1,pawrad)
  else if (radmesh(ivlocmesh)%mesh_type==radmesh(imainmesh)%mesh_type.and.&
 &         abs(radmesh(ivlocmesh)%rstep-radmesh(imainmesh)%rstep)<tol10.and.&
 &         abs(radmesh(ivlocmesh)%lstep-radmesh(imainmesh)%lstep)<tol10) then
   allocate(work1(radmesh(ivlocmesh)%mesh_size))
   call nderiv_gen(work1,vlocr,1,radmesh(ivlocmesh))
   pawtab%dvlocdr(1:pawrad%mesh_size,1)=work1(1:pawrad%mesh_size)
   deallocate(work1)
  else
   msz=vloc_mesh%mesh_size;allocate(dvlocr(msz))
   call nderiv_gen(dvlocr,vlocr,1,vloc_mesh)
   if (vloc_mesh%radfact(1)>zero) then
    yp1=1._dp/12._dp/vloc_mesh%stepint/vloc_mesh%radfact(1) &
&    *(-25._dp*dvlocr(1)+48._dp*dvlocr(2)-36._dp*dvlocr(3)+16._dp*dvlocr(4)-3._dp*dvlocr(5))
   else
    yp1=(dvlocr(2)-dvlocr(1))/(vloc_mesh%rad(2)-vloc_mesh%rad(1))
   end if
   ypn=1._dp/12._dp/vloc_mesh%stepint &
&      *(3._dp*dvlocr(msz-4)-16._dp*dvlocr(msz-3)+36._dp*dvlocr(msz-2)-48._dp*dvlocr(msz-1) &
&       +25._dp*dvlocr(msz)) / vloc_mesh%radfact(msz)
   allocate(work1(msz),work2(msz),work3(msz));work3(:)=vloc_mesh%rad(:)
   call spline(work3,dvlocr,msz,yp1,ypn,work1,work2)
   call splint(msz,work3,dvlocr,work1,pawrad%mesh_size,pawrad%rad,pawtab%dvlocdr(:,1))
   deallocate(work1,work2,work3)
   deallocate(dvlocr)
  end if
  if (abs(pawtab%dvlocdr(1,1))>0.001_dp) then
   write(message, '(a,a,a,a,a,a)' ) ch10,&
&   ' psp7in: ERROR -',ch10,&
&   '  Vloc derivative must be zero at r=0 !',ch10,&
&   '  Action: check your pseudopotential file.'
   call wrtout(ab_out,message,'COLL');call wrtout(06,  message,'COLL')
   call leave_new('COLL')
  end if
! --- Third derivative
  msz=pawrad%mesh_size
  if (pawrad%radfact(1)>zero) then
   yp1=1._dp/12._dp/pawrad%stepint/pawrad%radfact(1) &
&   *(-25._dp*pawtab%dvlocdr(1,1)+48._dp*pawtab%dvlocdr(2,1)-36._dp*pawtab%dvlocdr(3,1)&
&     +16._dp*pawtab%dvlocdr(4,1)-3._dp*pawtab%dvlocdr(5,1))
  else
   yp1=(pawtab%dvlocdr(2,1)-pawtab%dvlocdr(1,1))/(pawrad%rad(2)-pawrad%rad(1))
  end if
  ypn=1._dp/12._dp/pawrad%stepint &
&     *( 3._dp*pawtab%dvlocdr(msz-4,1)-16._dp*pawtab%dvlocdr(msz-3,1)&
&      +36._dp*pawtab%dvlocdr(msz-2,1)-48._dp*pawtab%dvlocdr(msz-1,1) &
&      +25._dp*pawtab%dvlocdr(msz  ,1)) / pawrad%radfact(msz)
  allocate(work1(msz),work2(msz));work1(:)=pawrad%rad(:)
  call spline(work1,pawtab%dvlocdr(:,1),msz,yp1,ypn,pawtab%dvlocdr(:,2),work2)
  deallocate(work1,work2)
 end if

!==================================================
!Compute Ex-correlation energy for the core density

 allocate(work1(core_mesh%mesh_size),work2(core_mesh%mesh_size))
 allocate(work3(core_mesh%mesh_size))
 work1(:)=zero
 if (pawxcdev/=0) then
  call pawxcm(ncore,pawtab%exccore,yp1,ixc,1,tmp_lmselect,work3,1,4,&
&             core_mesh,work1,1,0,work2)
 else
  pawang_tmp%l_size_max=1;pawang_tmp%angl_size=1;tmp_lmselect(1)=1
  allocate(pawang_tmp%angwgth(1));pawang_tmp%angwgth(1)=1._dp
  allocate(pawang_tmp%ylmr(1,1));pawang_tmp%ylmr(1,1)=1._dp/sqrt(four_pi)
  call pawxc (ncore,pawtab%exccore,yp1,ixc,1,tmp_lmselect,work3,1,4,&
&             pawang_tmp,core_mesh,work1,1,0,work2)
  deallocate(pawang_tmp%angwgth,pawang_tmp%ylmr)
 end if
 deallocate(work1,work2,work3)

!==========================================================
!Free temporary allocated space

 do imsh=1,nmesh
  deallocate(radmesh(imsh)%rad)
  deallocate(radmesh(imsh)%radfact)
  deallocate(radmesh(imsh)%simfact)
 end do
 deallocate(radmesh)
 deallocate(vlocr,ncore,tncore)
 deallocate(tproj_mesh%rad,tproj_mesh%radfact,tproj_mesh%simfact)
 deallocate(core_mesh%rad,core_mesh%radfact,core_mesh%simfact)
 deallocate(vloc_mesh%rad,vloc_mesh%radfact,vloc_mesh%simfact)
 if(pawtab%shape_type==-1)deallocate(shpf_mesh%rad, shpf_mesh%radfact,shpf_mesh%simfact)
 if (reduced_vloc) deallocate(rvloc_mesh%rad,rvloc_mesh%radfact,rvloc_mesh%simfact,rvlocr)
 if (reduced_ncor) deallocate(rcore_mesh%rad,rcore_mesh%radfact,rcore_mesh%simfact,rtncor)

!DEBUG
!write(6,*)' psp7in : exit'
!ENDDEBUG

end subroutine psp7in
!!***

Generated by  Doxygen 1.6.0   Back to index