TABLE OF CONTENTS


ABINIT/m_pawpsp [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pawpsp

FUNCTION

  Module to read PAW atomic data

COPYRIGHT

  Copyright (C) 2012-2018 ABINIT group (MT, FJ,TR, GJ, FB, FrD, AF, GMR, DRH)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  FOR DEVELOPERS: in order to preserve the portability of libPAW library,
  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt

SOURCE

21 #include "libpaw.h"
22 
23 module m_pawpsp
24 
25  USE_DEFS
26  USE_MSG_HANDLING
27  USE_MPI_WRAPPERS
28  USE_MEMORY_PROFILING
29 
30  use m_libpaw_libxc
31 #if defined LIBPAW_HAVE_FOX
32  use fox_sax
33 #endif
34 
35  use m_libpaw_tools, only : libpaw_basename, libpaw_get_free_unit
36 
37  use m_pawang, only: pawang_type
38  use m_pawtab, only: pawtab_type, wvlpaw_type, wvlpaw_allocate, wvlpaw_rholoc_free, &
39 &                    pawtab_free, wvlpaw_free, wvlpaw_rholoc_nullify, pawtab_bcast, &
40 &                    pawtab_set_flags, wvlpaw_allocate, wvlpaw_free, wvlpaw_rholoc_nullify, &
41 &                    wvlpaw_rholoc_free
42  use m_pawxmlps, only: rdpawpsxml_core, paw_setup_t, paw_setuploc, paw_setup_free
43  use m_pawrad, only: pawrad_type, pawrad_init, pawrad_free, pawrad_copy, &
44 &      pawrad_bcast, pawrad_ifromr, simp_gen, nderiv_gen, bound_deriv, pawrad_deducer0, poisson
45  use m_paw_numeric, only: paw_splint, paw_spline, paw_smooth, paw_jbessel_4spline
46  use m_paw_atom, only: atompaw_shapebes, atompaw_vhnzc, atompaw_shpfun, &
47 &                     atompaw_dij0, atompaw_kij
48  use m_pawxc, only: pawxc, pawxcm
49  use m_paw_gaussfit, only: gaussfit_projector
50 
51  implicit none
52 
53  private
54 
55  public:: pawpsp_calc_d5         !calculate up to the 5th derivative
56  public:: pawpsp_main            !main routine to read psp
57  public:: pawpsp_nl              !make paw projector form factors f_l(q)
58  public:: pawpsp_read            !read psp from file
59  public:: pawpsp_read_header     !read header of psp file
60  public:: pawpsp_read_corewf     !read core wavefunction
61  public:: pawpsp_read_header_2   !reads pspversion, basis_size and lmn_size
62  public:: pawpsp_rw_atompaw      !read and writes ATOMPAW psp with gaussian |p>
63  public:: pawpsp_wvl             !wavelet and icoulomb>0 related operations
64  public:: pawpsp_wvl_calc        !wavelet related operations
65  public:: pawpsp_7in             !reads non-XML atomic data
66  public:: pawpsp_17in            !reads XML atomic data
67  public:: pawpsp_calc            !calculates atomic quantities from psp info
68  public:: pawpsp_read_header_xml !read header of psp file for XML
69  public:: pawpsp_read_pawheader  !read header variables from XML objects
70  public:: pawpsp_bcast           ! broadcast PAW psp data
71  public:: pawpsp_cg              !compute sin FFT transform of a density
72  public:: pawpsp_lo              !compute sin FFT transform of local potential
73 
74 ! Private procedures
75  private:: pawpsp_wvl_sin2gauss  !convert sin/cos to gaussians

m_pawpsp/pawpsp_17in [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_17in

FUNCTION

 Initialize pspcod=17 ("PAW  XML pseudopotentials"):
 continue to read the corresponding file and compute the form factors

INPUTS

  ipsp= id in the array of the currently read pseudo.
  ixc=exchange-correlation choice from main routine data file
  lmax=value of lmax mentioned at the second line of the psp file
  lnmax=max. number of (l,n) components over all type of psps
            angular momentum of nonlocal pseudopotential
  mmax=max number of pts in real space grid (already read in the psp file header)
  mqgrid_ff=dimension of q (or G) grid for nl form factors (array ffspl)
  mqgrid_vl=dimension of q (or G) grid for Vloc (array vlspl)
  pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  pspheads= header of the current pseudopotential
  qgrid_ff(psps%mqgrid_ff)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
  qgrid_vl(psps%mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for Vloc
  xclevel= XC functional level
  zion=nominal valence of atom as specified in psp file
  znucl=atomic number of atom as specified in input file to main routine

OUTPUT

  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
  ffspl(psps%mqgrid_ff,2,psps%lnmax)=form factor f_l(q) and second derivative
   from spline fit for each angular momentum and each projector;
  pawrad <type(pawrad_type)>=paw radial mesh and related data
  pawtab <type(pawtab_type)>=paw tabulated starting data
  vlspl(psps%mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
  wvl_crmult,wvl_frmult= variables definining the fine and coarse grids in a wavelets calculation
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
  xcccrc=XC core correction cutoff radius (bohr) from psp file

NOTES

  Spin-orbit not yet implemented (to be done)
  Comments:
  * 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 ; gl(r)=numeric (read from psp file)
  shape_type= 1 ; gl(r)=k(r).r^l; k(r)=exp[-(r/sigma)**lambda]
  shape_type= 2 ; gl(r)=k(r).r^l; k(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2 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

PARENTS

      m_pawpsp,pspatm

CHILDREN

SOURCE

2951 subroutine pawpsp_17in(epsatm,ffspl,icoulomb,ipsp,ixc,lmax,&
2952 & lnmax,mmax,mqgrid_ff,mqgrid_vl,pawpsp_header,pawrad,pawtab,&
2953 & pawxcdev, qgrid_ff,qgrid_vl,usewvl,usexcnhat_in,vlspl,xcccrc,&
2954 & xclevel,xc_denpos,zion,znucl)
2955 
2956 
2957 !This section has been created automatically by the script Abilint (TD).
2958 !Do not modify the following lines by hand.
2959 #undef ABI_FUNC
2960 #define ABI_FUNC 'pawpsp_17in'
2961 !End of the abilint section
2962 
2963  implicit none
2964 
2965 !Arguments ------------------------------------
2966 !scalars
2967  integer,intent(in) :: ipsp,ixc,lmax,lnmax,mqgrid_ff,mqgrid_vl,pawxcdev,usexcnhat_in
2968  integer,intent(inout) ::mmax
2969  integer,intent(in) :: xclevel,icoulomb,usewvl
2970  real(dp),intent(in) :: xc_denpos,zion,znucl
2971  real(dp),intent(out) :: epsatm,xcccrc
2972  type(pawpsp_header_type),intent(in) :: pawpsp_header
2973  type(pawrad_type),intent(inout) :: pawrad
2974  type(pawtab_type),intent(inout) :: pawtab
2975 !arrays
2976  real(dp),intent(in) :: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl)
2977  real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax)
2978  real(dp),intent(out) :: vlspl(mqgrid_vl,2)
2979 
2980 !Local variables ------------------------------
2981 !scalars
2982  integer :: has_v_minushalf,ib,icoremesh,il,ilm,ilmn,ilmn0,iln,imainmesh,imsh,iprojmesh
2983  integer :: ir,iread1,ishpfmesh,ivalemesh,ivlocmesh,j0lmn,jlm,pngau
2984  integer :: jlmn,jln,klmn,msz,nmesh,nval,pspversion,shft,sz10,usexcnhat,vlocopt
2985  real(dp), parameter :: rmax_vloc=10.0_dp
2986  real(dp) :: fourpi,occ,rc,yp1,ypn
2987  logical :: save_core_msz
2988  character(len=500) :: msg
2989  type(pawrad_type) :: core_mesh,shpf_mesh,tproj_mesh,vale_mesh,vloc_mesh
2990 !arrays
2991  integer,allocatable :: mesh_shift(:),nprj(:)
2992  real(dp),allocatable :: kij(:),ncore(:)
2993  real(dp),allocatable :: shpf(:,:),tncore(:),tnvale(:),tproj(:,:),vhnzc(:),vlocr(:)
2994  real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:)
2995  type(pawrad_type),allocatable :: radmesh(:)
2996 
2997 !************************************************************************
2998 
2999 !==========================================================
3000 !Destroy everything in pawtab but optional flags
3001  call pawtab_free(pawtab)
3002 !Destroy everything in pawrad
3003  call pawrad_free(pawrad)
3004 
3005 !==========================================================
3006 !Initialize useful data
3007 
3008  pawtab%usexcnhat=usexcnhat_in
3009  fourpi=4*acos(-1.d0)
3010  pspversion=pawpsp_header%pawver
3011  save_core_msz=(usewvl==1 .or. icoulomb .ne. 0)
3012 
3013 !==========================================================
3014 !Initialize partial waves quantum numbers
3015 
3016  pawtab%basis_size=pawpsp_header%basis_size
3017  LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size))
3018  do ib=1,pawtab%basis_size
3019    pawtab%orbitals(ib)=paw_setuploc%valence_states%state(ib)%ll
3020  end do
3021 
3022 !==========================================================
3023 !Initialize various dims and indexes
3024 
3025  pawtab%lmn_size=pawpsp_header%lmn_size
3026  pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2
3027  pawtab%l_size=2*maxval(pawtab%orbitals)+1
3028  pawtab%ij_size=pawtab%basis_size*(pawtab%basis_size+1)/2
3029 
3030 !indlmn calculation (indices for (l,m,n) basis)
3031  if (allocated(pawtab%indlmn)) then
3032    LIBPAW_DEALLOCATE(pawtab%indlmn)
3033  end if
3034  LIBPAW_ALLOCATE(pawtab%indlmn,(6,pawtab%lmn_size))
3035  pawtab%indlmn(:,:)=0
3036  LIBPAW_BOUND1_ALLOCATE(nprj,BOUNDS(0,maxval(pawtab%orbitals)))
3037  ilmn=0;iln=0;nprj=0
3038  do ib=1,pawtab%basis_size
3039    il=pawtab%orbitals(ib)
3040    nprj(il)=nprj(il)+1
3041    iln=iln+1
3042    do ilm=1,2*il+1
3043      pawtab%indlmn(1,ilmn+ilm)=il
3044      pawtab%indlmn(2,ilmn+ilm)=ilm-(il+1)
3045      pawtab%indlmn(3,ilmn+ilm)=nprj(il)
3046      pawtab%indlmn(4,ilmn+ilm)=il*il+ilm
3047      pawtab%indlmn(5,ilmn+ilm)=iln
3048      pawtab%indlmn(6,ilmn+ilm)=1
3049    end do
3050    ilmn=ilmn+2*il+1
3051  end do
3052  LIBPAW_DEALLOCATE(nprj)
3053 !Are ilmn (found here) and pawtab%lmn_size compatibles ?
3054  if (ilmn/=pawtab%lmn_size) then
3055    write(msg, '(a,a,a,a,a)' )&
3056 &   'Calculated lmn size differs from',ch10,&
3057 &   'lmn_size read from pseudo !',ch10,&
3058 &   'Action: check your pseudopotential file.'
3059    MSG_ERROR(msg)
3060  end if
3061 
3062 !==========================================================
3063 !Read and initialize radial meshes
3064 
3065  nmesh=paw_setuploc%ngrid
3066  LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh))
3067  LIBPAW_ALLOCATE(mesh_shift,(nmesh))
3068  do imsh=1,nmesh
3069    radmesh(imsh)%mesh_type=-1
3070    radmesh(imsh)%rstep=zero
3071    radmesh(imsh)%lstep=zero
3072    mesh_shift(imsh)=0
3073    select case(trim(paw_setuploc%radial_grid(imsh)%eq))
3074      case("r=a*exp(d*i)")
3075        mesh_shift(imsh)=1
3076        radmesh(imsh)%mesh_type=3
3077        radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend &
3078 &                             -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh)
3079        radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa
3080        radmesh(imsh)%lstep=paw_setuploc%radial_grid(imsh)%dd
3081      case("r=a*i/(1-b*i)")
3082        write(msg, '(3a)' )&
3083 &       'The grid r=a*i/(1-b*i) is not implemented in ABINIT !',ch10,&
3084 &       'Action: check your psp file.'
3085        MSG_ERROR(msg)
3086      case("r=a*i/(n-i)")
3087        mesh_shift(imsh)=0
3088        radmesh(imsh)%mesh_type=5
3089        radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend &
3090 &                             -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh)
3091        radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa
3092        radmesh(imsh)%lstep=dble(paw_setuploc%radial_grid(imsh)%nn)
3093      case("r=a*(exp(d*i)-1)")
3094        mesh_shift(imsh)=0
3095        radmesh(imsh)%mesh_type=2
3096        radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend &
3097 &                             -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh)
3098        if(paw_setuploc%radial_grid(imsh)%istart==1)radmesh(imsh)%mesh_size=radmesh(imsh)%mesh_size+1
3099        radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa
3100        radmesh(imsh)%lstep=paw_setuploc%radial_grid(imsh)%dd
3101      case("r=d*i")
3102        mesh_shift(imsh)=0
3103        radmesh(imsh)%mesh_type=1
3104        radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend &
3105 &                             -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh)
3106        if(paw_setuploc%radial_grid(imsh)%istart==1)radmesh(imsh)%mesh_size=radmesh(imsh)%mesh_size+1
3107        radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%dd
3108      case("r=(i/n+a)^5/a-a^4")
3109        write(msg, '(3a)' )&
3110 &       'The grid r=(i/n+a)^5/a-a^4 is not implemented in ABINIT !',ch10,&
3111 &       'Action: check your psp file.'
3112        MSG_ERROR(msg)
3113    end select
3114  end do
3115 
3116 !Initialize radial meshes
3117  do imsh=1,nmesh
3118    call pawrad_init(radmesh(imsh))
3119  end do
3120 
3121  pawtab%rpaw=pawpsp_header%rpaw
3122 
3123 !==========================================================
3124 !Here reading shapefunction parameters
3125 
3126  pawtab%shape_type=pawpsp_header%shape_type
3127  pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99
3128  pawtab%rshp=pawpsp_header%rshp
3129  pawtab%shape_lambda=paw_setuploc%shape_function%lamb
3130  if(trim(paw_setuploc%shape_function%gtype)=="gauss")pawtab%shape_lambda=2
3131  pawtab%shape_sigma=paw_setuploc%shape_function%rc
3132 !If shapefunction type is gaussian, check exponent
3133  if (pawtab%shape_type==1) then
3134    if (pawtab%shape_lambda<2) then
3135      write(msg, '(3a)' )&
3136 &     'For a gaussian shape function, exponent lambda must be >1 !',ch10,&
3137 &     'Action: check your psp file.'
3138      MSG_ERROR(msg)
3139    end if
3140  end if
3141 
3142 !If shapefunction type is Bessel, deduce here its parameters from rc
3143  if (pawtab%shape_type==3) then
3144    LIBPAW_ALLOCATE(pawtab%shape_alpha,(2,pawtab%l_size))
3145    LIBPAW_ALLOCATE(pawtab%shape_q,(2,pawtab%l_size))
3146    rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw
3147    do il=1,pawtab%l_size
3148      call atompaw_shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc)
3149    end do
3150  end if
3151 
3152 !==========================================================
3153 !Mirror pseudopotential parameters to the output and log files
3154 
3155  write(msg,'(a,i2)')' Pseudopotential format is: paw',pspversion
3156  call wrtout(ab_out,msg,'COLL')
3157  call wrtout(std_out,  msg,'COLL')
3158  write(msg,'(2(a,i3),a,64i4)') &
3159 & ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',&
3160 & pawtab%lmn_size,'), orbitals=',pawtab%orbitals(1:pawtab%basis_size)
3161  call wrtout(ab_out,msg,'COLL')
3162  call wrtout(std_out,  msg,'COLL')
3163  write(msg,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw
3164  call wrtout(ab_out,msg,'COLL')
3165  call wrtout(std_out,  msg,'COLL')
3166  write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:'
3167  call wrtout(ab_out,msg,'COLL')
3168  call wrtout(std_out,  msg,'COLL')
3169 
3170  do imsh=1,nmesh
3171    if (radmesh(imsh)%mesh_type==1) &
3172 &   write(msg,'(a,i1,a,i4,a,g12.5)') &
3173 &   '  - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,&
3174 &   ' , step=',radmesh(imsh)%rstep
3175    if (radmesh(imsh)%mesh_type==2) &
3176 &   write(msg,'(a,i1,a,i4,2(a,g12.5))') &
3177 &   '  - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,&
3178 &   ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep
3179    if (radmesh(imsh)%mesh_type==3) &
3180 &   write(msg,'(a,i1,a,i4,2(a,g12.5))') &
3181 &   '  - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,&
3182 &   ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep
3183    if (radmesh(imsh)%mesh_type==4) &
3184 &   write(msg,'(a,i1,a,i4,a,g12.5)') &
3185 &   '  - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,&
3186 &   ' , AA=',radmesh(imsh)%rstep
3187    if (radmesh(imsh)%mesh_type==5) &
3188 &   write(msg,'(a,i1,a,i4,2(a,g12.5))') &
3189 &   '  - mesh ',imsh,': r(i)=-AA*i/(NN-i)), n=size=',radmesh(imsh)%mesh_size,&
3190 &   ' , AA=',radmesh(imsh)%rstep,' NN=',radmesh(imsh)%lstep
3191    call wrtout(ab_out,msg,'COLL')
3192    call wrtout(std_out,  msg,'COLL')
3193  end do
3194  if (pawtab%shape_type==-1) then
3195    write(msg,'(a)')&
3196    ' Shapefunction is NUMERIC type: directly read from atomic data file'
3197    call wrtout(ab_out,msg,'COLL')
3198    call wrtout(std_out,  msg,'COLL')
3199  end if
3200  if (pawtab%shape_type==1) then
3201    write(msg,'(2a,a,f6.3,a,i3)')&
3202 &   ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,&
3203 &   '                            with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda
3204    call wrtout(ab_out,msg,'COLL')
3205    call wrtout(std_out,  msg,'COLL')
3206  end if
3207  if (pawtab%shape_type==2) then
3208    write(msg,'(a)')&
3209    ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2'
3210    call wrtout(ab_out,msg,'COLL')
3211    call wrtout(std_out,  msg,'COLL')
3212  end if
3213  if (pawtab%shape_type==3) then
3214    write(msg,'(a)')&
3215 &   ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)'
3216    call wrtout(ab_out,msg,'COLL')
3217    call wrtout(std_out,  msg,'COLL')
3218  end if
3219  if (pawtab%rshp<1.d-8) then
3220    write(msg,'(a)') ' Radius for shape functions = sphere core radius'
3221  else
3222    write(msg,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp
3223  end if
3224  call wrtout(ab_out,msg,'COLL')
3225  call wrtout(std_out,  msg,'COLL')
3226 
3227 !==========================================================
3228 !Perfom tests
3229 
3230 !Are lmax and orbitals compatibles ?
3231  if (lmax/=maxval(pawtab%orbitals)) then
3232    write(msg, '(a,a,a)' )&
3233 &   'lmax /= MAX(orbitals) !',ch10,&
3234 &   'Action: check your pseudopotential file.'
3235    MSG_ERROR(msg)
3236  end if
3237 
3238 !Only mesh_type=1,2, 3 or 5 allowed
3239  do imsh=1,nmesh
3240    if (radmesh(imsh)%mesh_type>5) then
3241      write(msg, '(a,a,a)' )&
3242 &     'Only mesh types 1,2,3 or 5 allowed !',ch10,&
3243 &     'Action : check your pseudopotential or input file.'
3244      MSG_ERROR(msg)
3245    end if
3246  end do
3247 
3248 !==========================================================
3249 !Read tabulated atomic data
3250 
3251 !---------------------------------
3252 !Read wave-functions (phi)
3253 
3254  do ib=1,pawtab%basis_size
3255 
3256    if (ib==1) then
3257      do imsh=1,nmesh
3258        if(trim(paw_setuploc%ae_partial_wave(1)%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3259          mmax=radmesh(imsh)%mesh_size
3260          call pawrad_init(pawrad,mesh_size=mmax,mesh_type=radmesh(imsh)%mesh_type, &
3261 &         rstep=radmesh(imsh)%rstep,lstep=radmesh(imsh)%lstep,r_for_intg=pawtab%rpaw)
3262          pawtab%partialwave_mesh_size=pawrad%mesh_size
3263          pawtab%mesh_size=pawrad_ifromr(pawrad,pawtab%rpaw)+5
3264          pawtab%mesh_size=min(pawtab%mesh_size,pawrad%mesh_size)
3265          imainmesh=imsh
3266          exit
3267        end if
3268      end do
3269      LIBPAW_ALLOCATE(pawtab%phi,(pawtab%partialwave_mesh_size,pawtab%basis_size))
3270    else if (trim(paw_setuploc%ae_partial_wave(ib)%grid)/=trim(paw_setuploc%radial_grid(imainmesh)%id)) then
3271      write(msg, '(a,a,a)' )&
3272 &     'All Phi and tPhi must be given on the same radial mesh !',ch10,&
3273 &     'Action: check your pseudopotential file.'
3274      MSG_ERROR(msg)
3275    end if
3276    shft=mesh_shift(imainmesh)
3277    pawtab%phi(1+shft:pawtab%partialwave_mesh_size,ib)= &
3278 &         paw_setuploc%ae_partial_wave(ib)%data(1:pawtab%partialwave_mesh_size-shft) &
3279 &        *pawrad%rad(1+shft:pawtab%partialwave_mesh_size)
3280    if (shft==1) pawtab%phi(1,ib)=zero
3281  end do
3282  write(msg,'(a,i4)') &
3283 & ' mmax= ',mmax
3284  call wrtout(ab_out,msg,'COLL')
3285  call wrtout(std_out,  msg,'COLL')
3286 
3287 !---------------------------------
3288 !Read pseudo wave-functions (tphi)
3289 
3290  LIBPAW_ALLOCATE(pawtab%tphi,(pawtab%partialwave_mesh_size,pawtab%basis_size))
3291  do ib=1,pawtab%basis_size
3292 
3293    if(trim(paw_setuploc%pseudo_partial_wave(ib)%grid)/=trim(paw_setuploc%radial_grid(imainmesh)%id)) then
3294      write(msg, '(a,a,a)' )&
3295 &     'All Phi and tPhi must be given on the same radial mesh !',ch10,&
3296 &     'Action: check your pseudopotential file.'
3297      MSG_ERROR(msg)
3298    end if
3299    shft=mesh_shift(imainmesh)
3300    pawtab%tphi(1+shft:pawtab%partialwave_mesh_size,ib)=&
3301 &         paw_setuploc%pseudo_partial_wave(ib)%data(1:pawtab%partialwave_mesh_size-shft) &
3302 &        *pawrad%rad(1+shft:pawtab%partialwave_mesh_size)
3303    if (shft==1) pawtab%tphi(1,ib)=zero
3304  end do
3305  write(msg,'(a,i1)') &
3306 & ' Radial grid used for partial waves is grid ',imainmesh
3307  call wrtout(ab_out,msg,'COLL')
3308  call wrtout(std_out,  msg,'COLL')
3309 
3310 !---------------------------------
3311 !Read projectors (tproj)
3312 
3313  if (allocated(paw_setuploc%projector_fit)) then
3314     call wvlpaw_allocate(pawtab%wvl)
3315     LIBPAW_ALLOCATE(pawtab%wvl%pngau,(pawtab%basis_size))
3316     do ib=1,pawtab%basis_size
3317        pawtab%wvl%pngau(ib) = paw_setuploc%projector_fit(ib)%ngauss
3318     end do
3319     pawtab%wvl%ptotgau = sum(pawtab%wvl%pngau) * 2
3320     LIBPAW_ALLOCATE(pawtab%wvl%parg,(2,pawtab%wvl%ptotgau))
3321     LIBPAW_ALLOCATE(pawtab%wvl%pfac,(2,pawtab%wvl%ptotgau))
3322     pngau = 1
3323     do ib=1,pawtab%basis_size
3324        ! Complex gaussian
3325        pawtab%wvl%parg(:,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = &
3326             & paw_setuploc%projector_fit(ib)%expos(:,1:pawtab%wvl%pngau(ib))
3327        pawtab%wvl%pfac(:,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = &
3328             & paw_setuploc%projector_fit(ib)%factors(:,1:pawtab%wvl%pngau(ib))
3329        pngau = pngau + pawtab%wvl%pngau(ib)
3330        ! Conjugate gaussian
3331        pawtab%wvl%parg(1,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = &
3332             & paw_setuploc%projector_fit(ib)%expos(1,1:pawtab%wvl%pngau(ib))
3333        pawtab%wvl%parg(2,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = &
3334             & -paw_setuploc%projector_fit(ib)%expos(2,1:pawtab%wvl%pngau(ib))
3335        pawtab%wvl%pfac(1,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = &
3336             & paw_setuploc%projector_fit(ib)%factors(1,1:pawtab%wvl%pngau(ib))
3337        pawtab%wvl%pfac(2,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = &
3338             & -paw_setuploc%projector_fit(ib)%factors(2,1:pawtab%wvl%pngau(ib))
3339        pngau = pngau + pawtab%wvl%pngau(ib)
3340        pawtab%wvl%pngau(ib) = pawtab%wvl%pngau(ib) * 2
3341     end do
3342     pawtab%has_wvl=2
3343  else
3344     !Nullify wavelet objects for safety:
3345     pawtab%has_wvl=0
3346     call wvlpaw_free(pawtab%wvl)
3347  end if
3348 
3349  do ib=1,pawtab%basis_size
3350    if (ib==1) then
3351      do imsh=1,nmesh
3352        if(trim(paw_setuploc%projector_function(1)%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3353          iprojmesh=imsh
3354          exit
3355        end if
3356      end do
3357      call pawrad_copy(radmesh(iprojmesh),tproj_mesh)
3358      LIBPAW_ALLOCATE(tproj,(tproj_mesh%mesh_size,pawtab%basis_size))
3359    else if (trim(paw_setuploc%projector_function(ib)%grid)/=trim(paw_setuploc%radial_grid(iprojmesh)%id)) then
3360      write(msg, '(a,a,a)' )&
3361 &     'All tprojectors must be given on the same radial mesh !',ch10,&
3362 &     'Action: check your pseudopotential file.'
3363      MSG_ERROR(msg)
3364    end if
3365    shft=mesh_shift(iprojmesh)
3366    tproj(1+shft:tproj_mesh%mesh_size,ib)=paw_setuploc%projector_function(ib)%data(1:tproj_mesh%mesh_size-shft)&
3367 &   *tproj_mesh%rad(1+shft:tproj_mesh%mesh_size)
3368    if (shft==1) tproj(1,ib)=zero
3369  end do
3370  write(msg,'(a,i1)') &
3371 & ' Radial grid used for projectors is grid ',iprojmesh
3372  call wrtout(ab_out,msg,'COLL')
3373  call wrtout(std_out,  msg,'COLL')
3374 
3375 
3376 !---------------------------------
3377 !Read core density (coredens)
3378  do imsh=1,nmesh
3379    if(trim(paw_setuploc%ae_core_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3380      icoremesh=imsh
3381      exit
3382    end if
3383  end do
3384  call pawrad_copy(radmesh(icoremesh),core_mesh)
3385  if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.&
3386 & (radmesh(icoremesh)%rstep    /=pawrad%rstep)    .or.&
3387 & (radmesh(icoremesh)%lstep    /=pawrad%lstep)) then
3388    write(msg, '(a,a,a,a,a)' )&
3389 &   'Ncore must be given on a radial mesh with the same',ch10,&
3390 &   'type and step(s) than the main radial mesh (mesh for Phi) !',ch10,&
3391 &   'Action: check your pseudopotential file.'
3392    MSG_ERROR(msg)
3393  end if
3394  LIBPAW_ALLOCATE(ncore,(core_mesh%mesh_size))
3395  shft=mesh_shift(icoremesh)
3396  ncore(1+shft:core_mesh%mesh_size)=paw_setuploc%ae_core_density%data(1:core_mesh%mesh_size-shft)/sqrt(fourpi)
3397  if (shft==1) call pawrad_deducer0(ncore,core_mesh%mesh_size,core_mesh)
3398 
3399 !Construct and save VH[z_NC] if requested
3400  if (pawtab%has_vhnzc==1) then
3401    LIBPAW_ALLOCATE(pawtab%VHnZC,(pawtab%mesh_size))
3402    LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size))
3403    call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl)
3404    pawtab%VHnZC(1:pawtab%mesh_size)=vhnzc(1:pawtab%mesh_size)
3405    pawtab%has_vhnzc=2
3406    LIBPAW_DEALLOCATE(vhnzc)
3407  end if
3408 
3409  pawtab%core_mesh_size=pawtab%mesh_size
3410  if(save_core_msz) pawtab%core_mesh_size=core_mesh%mesh_size
3411  LIBPAW_ALLOCATE(pawtab%coredens,(pawtab%core_mesh_size))
3412  pawtab%rcore=core_mesh%rad(pawtab%core_mesh_size)
3413  pawtab%coredens(1:pawtab%core_mesh_size)=ncore(1:pawtab%core_mesh_size)
3414 
3415 !---------------------------------
3416 !Read pseudo core density (tcoredens)
3417  do imsh=1,nmesh
3418    if(trim(paw_setuploc%pseudo_core_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3419      iread1=imsh
3420      exit
3421    end if
3422  end do
3423  if (iread1/=icoremesh) then
3424    write(msg, '(a,a,a,a,a,a,a,a)' )&
3425 &   'Pseudized core density (tNcore) must be given',ch10,&
3426 &   'on the same radial mesh as core density (Ncore) !',ch10,&
3427 &   'Action: check your pseudopotential file.'
3428    MSG_ERROR(msg)
3429  end if
3430  LIBPAW_ALLOCATE(tncore,(core_mesh%mesh_size))
3431  shft=mesh_shift(icoremesh)
3432  tncore(1+shft:core_mesh%mesh_size)=paw_setuploc%pseudo_core_density%data(1:core_mesh%mesh_size-shft)/sqrt(fourpi)
3433  if (shft==1) call pawrad_deducer0(tncore,core_mesh%mesh_size,core_mesh)
3434  if(save_core_msz)  then
3435    LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,6))
3436  else
3437    LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,1))
3438  end if
3439  if (maxval(abs(tncore(:)))<tol6) then
3440    pawtab%usetcore=0
3441    pawtab%tcoredens(1:pawtab%core_mesh_size,:)=zero
3442  else
3443    pawtab%usetcore=1
3444    pawtab%tcoredens(1:pawtab%core_mesh_size,1)=tncore(1:pawtab%core_mesh_size)
3445  end if
3446  write(msg,'(a,i1)') &
3447 & ' Radial grid used for (t)core density is grid ',icoremesh
3448  call wrtout(ab_out,msg,'COLL')
3449  call wrtout(std_out,  msg,'COLL')
3450 
3451 
3452 !---------------------------------
3453 !Read local pseudopotential=Vh(tn_zc) or Vbare
3454 
3455  if ((paw_setuploc%blochl_local_ionic_potential%tread).and.&
3456 & (pawtab%usexcnhat==-1.or.pawtab%usexcnhat==0.or.(pawtab%usexcnhat==1.and.&
3457 & ((.not.paw_setuploc%zero_potential%tread).or.(.not.paw_setuploc%kresse_joubert_local_ionic_potential%tread))))) then
3458    usexcnhat=0;vlocopt=2
3459    do imsh=1,nmesh
3460      if(trim(paw_setuploc%blochl_local_ionic_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3461        iread1=imsh
3462        exit
3463      end if
3464    end do
3465    ivlocmesh=iread1
3466    call pawrad_copy(radmesh(ivlocmesh),vloc_mesh)
3467    LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size))
3468    shft=mesh_shift(ivlocmesh)
3469    vlocr(1+shft:vloc_mesh%mesh_size)=paw_setuploc%blochl_local_ionic_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi)
3470    if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh)
3471  else if((paw_setuploc%kresse_joubert_local_ionic_potential%tread).and.&
3472 &   (pawtab%usexcnhat==-1.or.pawtab%usexcnhat==1.or.(pawtab%usexcnhat==0.and.&
3473 &   (.not.paw_setuploc%zero_potential%tread)))) then
3474    usexcnhat=1;vlocopt=1
3475    do imsh=1,nmesh
3476      if(trim(paw_setuploc%kresse_joubert_local_ionic_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3477        iread1=imsh
3478        exit
3479      end if
3480    end do
3481    ivlocmesh=iread1
3482    call pawrad_copy(radmesh(ivlocmesh),vloc_mesh)
3483    LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size))
3484    shft=mesh_shift(ivlocmesh)
3485    vlocr(1+shft:vloc_mesh%mesh_size)= &
3486 &   paw_setuploc%kresse_joubert_local_ionic_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi)
3487    if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh)
3488  else if(paw_setuploc%zero_potential%tread) then
3489    usexcnhat=0;vlocopt=0
3490    do imsh=1,nmesh
3491      if(trim(paw_setuploc%zero_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3492        iread1=imsh
3493        exit
3494      end if
3495    end do
3496    ivlocmesh=iread1
3497    vloc_mesh%mesh_type=radmesh(ivlocmesh)%mesh_type
3498    vloc_mesh%rstep=radmesh(ivlocmesh)%rstep
3499    vloc_mesh%lstep=radmesh(ivlocmesh)%lstep
3500    vloc_mesh%mesh_size=pawrad_ifromr(radmesh(ivlocmesh),rmax_vloc)
3501    call pawrad_init(vloc_mesh)
3502 !   call pawrad_copy(radmesh(ivlocmesh),vloc_mesh)
3503    LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size))
3504    vlocr=zero
3505    shft=mesh_shift(ivlocmesh)
3506    vlocr(1+shft:vloc_mesh%mesh_size)=paw_setuploc%zero_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi)
3507    if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh)
3508  else
3509    write(msg, '(a,a,a,a,a)' )&
3510 &   'At least one local potential must be given',ch10,&
3511 &   'Action: check your pseudopotential file.'
3512    MSG_ERROR(msg)
3513  end if
3514 
3515  write(msg,'(a,i1)') &
3516 & ' Radial grid used for Vloc is grid ',ivlocmesh
3517  call wrtout(ab_out,msg,'COLL')
3518  call wrtout(std_out,  msg,'COLL')
3519 
3520 !-------------------------------------------------
3521 !Read LDA-1/2 potential
3522  if (paw_setuploc%LDA_minus_half_potential%tread) then
3523    do imsh=1,nmesh
3524      if(trim(paw_setuploc%LDA_minus_half_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3525        iread1=imsh
3526        exit
3527      end if
3528    end do
3529    if(iread1/=ivlocmesh) then
3530      write(msg, '(a)' )&
3531 &     'The LDA-1/2 potential must be given on the same grid as the local potential.'
3532      MSG_ERROR(msg)
3533    end if
3534    has_v_minushalf=1
3535    LIBPAW_ALLOCATE(pawtab%vminushalf,(vloc_mesh%mesh_size))
3536    shft=mesh_shift(ivlocmesh)
3537    pawtab%vminus_mesh_size=vloc_mesh%mesh_size
3538    pawtab%vminushalf(1+shft:vloc_mesh%mesh_size)= &
3539 &   paw_setuploc%LDA_minus_half_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi)
3540    if (shft==1) call pawrad_deducer0(pawtab%vminushalf,vloc_mesh%mesh_size,vloc_mesh)
3541    write(msg,'(a,i1)') &
3542 &   ' Radial grid used for LDA-1/2 potential is grid ',ivlocmesh
3543    call wrtout(ab_out,msg,'COLL')
3544    call wrtout(std_out,  msg,'COLL')
3545  else
3546    LIBPAW_ALLOCATE(pawtab%vminushalf,(0))
3547    has_v_minushalf=0
3548  end if
3549  if(has_v_minushalf==0.and.pawtab%has_vminushalf==1) then
3550    write(msg, '(a)' )&
3551 &     'The LDA-1/2 potential must be given in the XML PAW datafile.'
3552    MSG_ERROR(msg)
3553  end if
3554 
3555 !---------------------------------
3556 !Eventually read "numeric" shapefunctions (if shape_type=-1)
3557  if (pawtab%shape_type==-1) then
3558    LIBPAW_ALLOCATE(pawtab%shapefunc,(pawtab%mesh_size,pawtab%l_size))
3559    do imsh=1,nmesh
3560      if(trim(paw_setuploc%shape_function%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3561        iread1=imsh
3562        exit
3563      end if
3564    end do
3565    call pawrad_copy(radmesh(iread1),shpf_mesh)
3566    ishpfmesh=iread1
3567    LIBPAW_ALLOCATE(shpf,(shpf_mesh%mesh_size,pawtab%l_size))
3568    shft=mesh_shift(ishpfmesh)
3569    shpf(1,1)=one
3570    do ir=2,shpf_mesh%mesh_size
3571      shpf(ir,1)=paw_setuploc%shape_function%data(ir-shft,1)
3572    end do
3573    sz10=size(paw_setuploc%shape_function%data,2)
3574    if(sz10>=2) then
3575      do il=2,pawtab%l_size
3576        shpf(1,il)=zero
3577        do ir=2,shpf_mesh%mesh_size
3578          shpf(ir,il)=paw_setuploc%shape_function%data(ir-shft,il)
3579        end do
3580      end do
3581    else
3582      do il=2,pawtab%l_size
3583        shpf(1,il)=zero
3584        do ir=2,shpf_mesh%mesh_size
3585          shpf(ir,il)=paw_setuploc%shape_function%data(ir-shft,1)*shpf_mesh%rad(ir)**(il-1)
3586        end do
3587      end do
3588    end if
3589    write(msg,'(a,i1)') &
3590 &   ' Radial grid used for shape functions is grid ',iread1
3591    call wrtout(ab_out,msg,'COLL')
3592    call wrtout(std_out,  msg,'COLL')
3593 
3594 !  Has to spline shape functions if mesh is not the "main" mesh
3595    if (ishpfmesh/=imainmesh) then
3596      msz=shpf_mesh%mesh_size
3597      LIBPAW_ALLOCATE(work1,(msz))
3598      LIBPAW_ALLOCATE(work2,(msz))
3599      LIBPAW_ALLOCATE(work3,(msz))
3600      LIBPAW_ALLOCATE(work4,(pawrad%mesh_size))
3601      work3(1:msz)=shpf_mesh%rad(1:msz)
3602      work4(1:pawrad%mesh_size)=pawrad%rad(1:pawrad%mesh_size)
3603      do il=1,pawtab%l_size
3604        call bound_deriv(shpf(1:msz,il),shpf_mesh,msz,yp1,ypn)
3605        call paw_spline(work3,shpf(:,il),msz,yp1,ypn,work1)
3606        call paw_splint(msz,work3,shpf(:,il),work1,pawrad%mesh_size,work4,pawtab%shapefunc(:,il))
3607      end do
3608      LIBPAW_DEALLOCATE(work1)
3609      LIBPAW_DEALLOCATE(work2)
3610      LIBPAW_DEALLOCATE(work3)
3611      LIBPAW_DEALLOCATE(work4)
3612    else
3613      pawtab%shapefunc(:,:)=shpf(:,:)
3614    end if
3615    LIBPAW_DEALLOCATE(shpf)
3616  end if
3617 
3618 !---------------------------------
3619 !Read pseudo valence density
3620  if (paw_setuploc%pseudo_valence_density%tread) then
3621    do imsh=1,nmesh
3622      if(trim(paw_setuploc%pseudo_valence_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then
3623        iread1=imsh
3624        exit
3625      end if
3626    end do
3627    ivalemesh=iread1
3628    call pawrad_copy(radmesh(iread1),vale_mesh)
3629    LIBPAW_ALLOCATE(tnvale,(vale_mesh%mesh_size))
3630    shft=mesh_shift(ivalemesh)
3631    tnvale(1+shft:vale_mesh%mesh_size)=paw_setuploc%pseudo_valence_density%data(1:vale_mesh%mesh_size-shft)/sqrt(fourpi)
3632    if (shft==1) call pawrad_deducer0(tnvale,vale_mesh%mesh_size,vale_mesh)
3633    pawtab%has_tvale=1
3634    write(msg,'(a,i1)') &
3635 &   ' Radial grid used for pseudo valence density is grid ',ivalemesh
3636    call wrtout(ab_out,msg,'COLL')
3637    call wrtout(std_out,  msg,'COLL')
3638  else
3639    pawtab%has_tvale=0
3640    LIBPAW_ALLOCATE(tnvale,(0))
3641  end if
3642 
3643 !---------------------------------
3644 !Read initial guess of rhoij (rhoij0)
3645 
3646  LIBPAW_ALLOCATE(pawtab%rhoij0,(pawtab%lmn2_size))
3647  pawtab%rhoij0=zero
3648  ilmn0=0
3649  do ib=1,pawtab%basis_size
3650    il=2*pawtab%orbitals(ib)+1
3651    occ=paw_setuploc%valence_states%state(ib)%ff
3652    if (occ<zero)occ=zero
3653    do ilmn=ilmn0+1,ilmn0+il
3654      pawtab%rhoij0(ilmn*(ilmn+1)/2)=occ/dble(il)
3655    end do
3656    ilmn0=ilmn0+il
3657  end do
3658 
3659 !---------------------------------
3660 !Read Kij terms (kij0) and deduce eventually Dij0
3661 
3662  LIBPAW_ALLOCATE(kij,(pawtab%lmn2_size))
3663  kij=zero
3664  nval=paw_setuploc%valence_states%nval
3665  do jlmn=1,pawtab%lmn_size
3666    j0lmn=jlmn*(jlmn-1)/2
3667    jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn)
3668    do ilmn=1,jlmn
3669      klmn=j0lmn+ilmn
3670      ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn)
3671      if (ilm==jlm) kij(klmn)=paw_setuploc%kinetic_energy_differences%data(jln+(iln-1)*nval)
3672    end do
3673  end do
3674  if (vlocopt>0) then
3675    LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size))
3676    if (size(pawtab%vminushalf)>0.and.pawtab%has_vminushalf==1) then
3677      vlocr(1:vloc_mesh%mesh_size)=vlocr(1:vloc_mesh%mesh_size)+pawtab%vminushalf(1:vloc_mesh%mesh_size)
3678    end if
3679    call atompaw_dij0(pawtab%indlmn,kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,&
3680 &                    vloc_mesh,vlocr,znucl)
3681  end if
3682 
3683 !Keep eventualy Kij in memory
3684  if (pawtab%has_kij==1.or.vlocopt==0) then
3685    LIBPAW_ALLOCATE(pawtab%kij,(pawtab%lmn2_size))
3686    pawtab%kij(:)=kij(:)
3687    if (vlocopt> 0) pawtab%has_kij=2
3688 !  This -1 means that pawtab%kij will be freed later
3689    if (vlocopt==0) pawtab%has_kij=-1
3690  end if
3691 
3692  LIBPAW_DEALLOCATE(kij)
3693 
3694 !---------------------------------
3695 !Read exact-exchange Fock terms for core-valence interactions (ex_cvij)
3696  if (paw_setuploc%exact_exchange_matrix%tread.eqv..true.) then
3697    pawtab%has_fock=2
3698    LIBPAW_ALLOCATE(pawtab%ex_cvij,(pawtab%lmn2_size))
3699    pawtab%ex_cvij=zero
3700    nval=paw_setuploc%valence_states%nval
3701    do jlmn=1,pawtab%lmn_size
3702      j0lmn=jlmn*(jlmn-1)/2
3703      jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn)
3704      do ilmn=1,jlmn
3705        klmn=j0lmn+ilmn
3706        ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn)
3707        if (ilm==jlm) pawtab%ex_cvij(klmn)=paw_setuploc%exact_exchange_matrix%data(jln+(iln-1)*nval)
3708      end do
3709    end do
3710    pawtab%ex_cc=paw_setuploc%ex_cc
3711  end if
3712 
3713 !==========================================================
3714 !Compute additional atomic data only depending on present DATASET
3715 
3716  call pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,ixc,lnmax,&
3717 &     mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,&
3718 &     qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,&
3719 &     vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl)
3720 
3721  if(usewvl==1 .or. icoulomb > 0) then
3722 !  Calculate up to the 5th derivative of tcoredens
3723    call pawpsp_calc_d5(core_mesh,pawtab%core_mesh_size,pawtab%tcoredens)
3724 !  Other wvl related operations
3725    call pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr)
3726  else if (pawtab%has_wvl>0) then
3727    call wvlpaw_rholoc_nullify(pawtab%wvl%rholoc)
3728  end if
3729 
3730 !==========================================================
3731 !Free temporary allocated space
3732 
3733  call pawrad_free(radmesh)
3734  LIBPAW_DATATYPE_DEALLOCATE(radmesh)
3735  LIBPAW_DEALLOCATE(mesh_shift)
3736 
3737  call pawrad_free(tproj_mesh)
3738  call pawrad_free(core_mesh)
3739  call pawrad_free(vloc_mesh)
3740 
3741  LIBPAW_DEALLOCATE(vlocr)
3742  LIBPAW_DEALLOCATE(ncore)
3743  LIBPAW_DEALLOCATE(tncore)
3744  LIBPAW_DEALLOCATE(tproj)
3745 
3746  if(pawtab%shape_type==-1) then
3747    call pawrad_free(shpf_mesh)
3748  end if
3749  if (paw_setuploc%pseudo_valence_density%tread) then
3750    call pawrad_free(vale_mesh)
3751  end if
3752 
3753  LIBPAW_DEALLOCATE(tnvale)
3754 
3755 end subroutine pawpsp_17in

m_pawpsp/pawpsp_7in [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_7in

FUNCTION

 Initialize pspcod=7 ("PAW pseudopotentials"):
 continue to read the corresponding file and compute the form factors

INPUTS

  icoulomb==0 : usual reciprocal space computation
           =1 : free boundary conditions are used
  ipsp=id in the array of the currently read pseudo.
  ixc=exchange-correlation choice from main routine data file
  lloc=angular momentum choice of local pseudopotential
  lmax=value of lmax mentioned at the second line of the psp file
  pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  xclevel= XC functional level
  zion=nominal valence of atom as specified in psp file

OUTPUT

  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
  ffspl(psps%mqgrid_ff,2,psps%lnmax)=form factor f_l(q) and second derivative
   from spline fit for each angular momentum and each projector;
  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
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
  xcccrc=XC core correction cutoff radius (bohr) from psp file

NOTES

  Spin-orbit not yet implemented (to be done)

PARENTS

      m_pawpsp,pspatm

CHILDREN

SOURCE

3799 subroutine pawpsp_7in(epsatm,ffspl,icoulomb,ixc,&
3800 & lmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,&
3801 & pawrad,pawtab,pawxcdev,qgrid_ff,qgrid_vl,&
3802 & usewvl,usexcnhat_in,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl)
3803 
3804 
3805 !This section has been created automatically by the script Abilint (TD).
3806 !Do not modify the following lines by hand.
3807 #undef ABI_FUNC
3808 #define ABI_FUNC 'pawpsp_7in'
3809 !End of the abilint section
3810 
3811  implicit none
3812 
3813 !Arguments ------------------------------------
3814 !scalars
3815  integer, intent(in):: icoulomb,ixc
3816  integer, intent(in):: lmax,lnmax,mmax
3817  integer, intent(in):: mqgrid_ff,mqgrid_vl,pawxcdev
3818  integer, intent(in):: usewvl,usexcnhat_in,xclevel
3819  real(dp), intent(in):: xc_denpos,zion,znucl
3820  real(dp), intent(out):: epsatm,xcccrc
3821  type(pawrad_type), intent(inout):: pawrad
3822  type(pawtab_type), intent(inout) :: pawtab
3823 !arrays
3824  real(dp),intent(in):: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl)
3825  real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax)
3826  real(dp),intent(out) :: vlspl(mqgrid_vl,2)
3827 
3828 !Local variables ------------------------------
3829 !scalars
3830  integer :: imainmesh
3831  integer :: nmesh
3832  integer :: pspversion,usexcnhat,vlocopt
3833  logical :: save_core_msz
3834  type(pawrad_type) :: core_mesh,tproj_mesh,vale_mesh,vloc_mesh
3835  real(dp),pointer :: ncore(:),tncore(:),tnvale(:),tproj(:,:),vlocr(:)
3836  type(pawrad_type),pointer :: radmesh(:)
3837 
3838 !************************************************************************
3839 
3840 !Destroy everything in pawtab but optional flags
3841  call pawtab_free(pawtab)
3842 !Destroy everything in pawrad
3843  call pawrad_free(pawrad)
3844 
3845  save_core_msz=(usewvl==1 .or. icoulomb .ne. 0)
3846  nullify(ncore);nullify(tncore);nullify(tnvale)
3847  nullify(tproj);nullify(vlocr)
3848  nullify(radmesh)
3849 
3850  call pawpsp_read(core_mesh,tmp_unit,imainmesh,lmax,&
3851 &  ncore,nmesh,pawrad,pawtab,pspversion,radmesh,save_core_msz,&
3852 &  tncore,tnvale,tproj,tproj_mesh,usexcnhat_in,usexcnhat,vale_mesh,vlocopt,vlocr,&
3853 &  vloc_mesh,znucl)
3854 
3855  call pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,ixc,lnmax,&
3856 &     mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,&
3857 &     qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,&
3858 &     vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl)
3859 
3860  if(usewvl==1 .or. icoulomb > 0) then
3861 !  Calculate up to the 5th derivative of tcoredens
3862    call pawpsp_calc_d5(core_mesh,pawtab%core_mesh_size,pawtab%tcoredens)
3863 !  Other wvl related operations
3864    call pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr)
3865  else if (pawtab%has_wvl>0) then
3866    call wvlpaw_rholoc_nullify(pawtab%wvl%rholoc)
3867  end if
3868 
3869 !==========================================================
3870 !Free temporary allocated space
3871  call pawrad_free(radmesh)
3872  call pawrad_free(tproj_mesh)
3873  call pawrad_free(core_mesh)
3874  call pawrad_free(vloc_mesh)
3875  LIBPAW_DATATYPE_DEALLOCATE(radmesh)
3876  if (associated(vlocr)) then
3877    LIBPAW_POINTER_DEALLOCATE(vlocr)
3878  end if
3879  if (associated(ncore)) then
3880    LIBPAW_POINTER_DEALLOCATE(ncore)
3881  end if
3882  if (associated(tncore)) then
3883    LIBPAW_POINTER_DEALLOCATE(tncore)
3884  end if
3885  if (associated(tnvale)) then
3886    LIBPAW_POINTER_DEALLOCATE(tnvale)
3887  end if
3888  if (associated(tproj)) then
3889    LIBPAW_POINTER_DEALLOCATE(tproj)
3890  end if
3891  if (pspversion>=4)  then
3892    call pawrad_free(vale_mesh)
3893  end if
3894 
3895 end subroutine pawpsp_7in

m_pawpsp/pawpsp_bcast [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_bcast

FUNCTION

 Communicate paw data to all processors

INPUTS

 comm_mpi= communicator used to broadcast data
 lnmax= Max. number of (l,n) components over all type of psps
 mqgrid_ff= dimension of ffspl
 mqgrid_vl= dimension of vlspl

OUTPUT

  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
  ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and derivative
  pawrad=<type pawrad_type>
  pawtab=<type pawtab_type>
  vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
  xcccrc=XC core correction cutoff radius (bohr) from psp file

PARENTS

      m_pawpsp,pspatm

CHILDREN

SOURCE

4654 subroutine pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc)
4655 
4656 
4657 !This section has been created automatically by the script Abilint (TD).
4658 !Do not modify the following lines by hand.
4659 #undef ABI_FUNC
4660 #define ABI_FUNC 'pawpsp_bcast'
4661 !End of the abilint section
4662 
4663  implicit none
4664 
4665 !Arguments ------------------------------------
4666  integer,intent(in) :: comm_mpi
4667  real(dp),intent(inout) :: epsatm,xcccrc
4668  real(dp),intent(inout) :: ffspl(:,:,:),vlspl(:,:)
4669  type(pawrad_type),intent(inout) :: pawrad
4670  type(pawtab_type),intent(inout) :: pawtab
4671 
4672 !Local variables-------------------------------
4673  integer :: ierr,ii,me,nn_dpr
4674  integer :: siz_ffspl,siz1_ffspl,siz2_ffspl,siz3_ffspl,siz_vlspl,siz1_vlspl,siz2_vlspl
4675  integer,allocatable :: list_int(:)
4676  real(dp),allocatable :: list_dpr(:)
4677 
4678 !*************************************************************************
4679 
4680  me=xmpi_comm_rank(comm_mpi)
4681 
4682 !Broadcast pawrad
4683  call pawrad_bcast(pawrad,comm_mpi)
4684 
4685 !Broadcast pawtab (only data read from file)
4686  call pawtab_bcast(pawtab,comm_mpi,only_from_file=.true.)
4687 
4688 !Broadcast the sizes of the arrays
4689  LIBPAW_ALLOCATE(list_int,(5))
4690  if (me==0) then
4691    siz1_vlspl=size(vlspl,1); list_int(1)=siz1_vlspl
4692    siz2_vlspl=size(vlspl,2); list_int(2)=siz2_vlspl
4693    siz1_ffspl=size(ffspl,1); list_int(3)=siz1_ffspl
4694    siz2_ffspl=size(ffspl,2); list_int(4)=siz2_ffspl
4695    siz3_ffspl=size(ffspl,3); list_int(5)=siz3_ffspl
4696  end if
4697  call xmpi_bcast(list_int,0,comm_mpi,ierr)
4698  if (me/=0) then
4699    siz1_vlspl=list_int(1)
4700    siz2_vlspl=list_int(2)
4701    siz1_ffspl=list_int(3)
4702    siz2_ffspl=list_int(4)
4703    siz3_ffspl=list_int(5)
4704  end if
4705  siz_vlspl=siz1_vlspl*siz2_vlspl
4706  siz_ffspl=siz1_ffspl*siz2_ffspl*siz3_ffspl
4707  LIBPAW_DEALLOCATE(list_int)
4708 
4709 !Broadcast the reals
4710  nn_dpr=2+siz_vlspl+siz_ffspl
4711  LIBPAW_ALLOCATE(list_dpr,(nn_dpr))
4712  if (me==0) then
4713    ii=1
4714    list_dpr(ii)=epsatm ;ii=ii+1
4715    list_dpr(ii)=xcccrc ;ii=ii+1
4716    list_dpr(ii:ii+siz_vlspl-1)=reshape(vlspl,(/siz_vlspl/)) ;ii=ii+siz_vlspl
4717    list_dpr(ii:ii+siz_ffspl-1)=reshape(ffspl,(/siz_ffspl/)) ;ii=ii+siz_ffspl
4718  end if
4719  call xmpi_bcast(list_dpr,0,comm_mpi,ierr)
4720  if (me/=0) then
4721    ii=1
4722    epsatm=list_dpr(ii) ;ii=ii+1
4723    xcccrc=list_dpr(ii) ;ii=ii+1
4724    vlspl=reshape(list_dpr(ii:ii+siz_vlspl-1),(/siz1_vlspl,siz2_vlspl/))
4725    ii=ii+siz_vlspl
4726    ffspl=reshape(list_dpr(ii:ii+siz_ffspl-1),(/siz1_ffspl,siz2_ffspl,siz3_ffspl/))
4727    ii=ii+siz_ffspl
4728  end if
4729  LIBPAW_DEALLOCATE(list_dpr)
4730 
4731 end subroutine pawpsp_bcast

m_pawpsp/pawpsp_calc [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_calc

FUNCTION

 Performs tests and compute data related to pspcod=7 or 17 ("PAW pseudopotentials")

INPUTS

  core_mesh<type(pawrad_type)>= radial mesh for the core density
  imainmesh= serial number of the main mesh
  ixc=exchange-correlation choice from main routine data file
  lnmax=max. number of (l,n) components over all type of psps
            angular momentum of nonlocal pseudopotential
  mqgrid_ff=dimension of q (or G) grid for nl form factors (array ffspl)
  mqgrid_vl=dimension of q (or G) grid for Vloc (array vlspl)
  ncore(core_mesh%mesh_size)= core density
  nmesh= number of radial meshes
  pawrad <type(pawrad_type)>=paw radial mesh and related data
  pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments)
  pspversion= version of the atompaw code used to generate paw data.
  qgrid_ff(mqgrid_ff)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
  qgrid_vl(mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for Vloc
  radmesh(nmesh)<type(pawrad_type)>=paw radial meshes and related data
  tncore(core_mesh%mesh_size)= pseudo core density
  tproj(tproj_mesh%mesh_size)= non-local projectors in real space
  tproj_mesh<type(pawrad_type)>= radial mesh for the projectors
  usexcnhat=0 if compensation charge density is not included in XC terms
            1 if compensation charge density is included in XC terms
  vale_mesh<type(pawrad_type)>= radial mesh for the valence density
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
  vlocopt= option for the local potential.(0=Vbare, 1=VH(tnzc) with hat in XC, 2=VH(tnzc) w/o hat in XC)
  vlocr(vloc_mesh%mesh_size)= local potential according to vlocopt.
  xclevel= XC functional level
  zion=nominal valence of atom as specified in psp file
  znucl=atomic number of atom as specified in input file to main routine

OUTPUT

  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;
  vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
  xcccrc=XC core correction cutoff radius (bohr) from psp file

SIDE EFFECTS

  pawtab <type(pawtab_type)>=paw tabulated starting data
  tnvale(vale_mesh%mesh_size)= pseudo valence density (+ nhat in output)
  vloc_mesh<type(pawrad_type)>= radial mesh for the local potential

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

1757 subroutine pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,ixc,lnmax,&
1758 &          mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,&
1759 &          qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,&
1760 &          vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl)
1761 
1762 
1763 !This section has been created automatically by the script Abilint (TD).
1764 !Do not modify the following lines by hand.
1765 #undef ABI_FUNC
1766 #define ABI_FUNC 'pawpsp_calc'
1767 !End of the abilint section
1768 
1769  implicit none
1770 
1771 !Arguments ------------------------------------
1772 !scalars
1773  integer,intent(in) :: imainmesh,ixc,lnmax,mqgrid_ff,mqgrid_vl
1774  integer,intent(in) :: nmesh,pawxcdev,pspversion,usexcnhat,vlocopt
1775  integer,intent(in) ::mmax
1776  integer,intent(in) :: xclevel
1777  real(dp),intent(in) :: xc_denpos,zion,znucl
1778  real(dp),intent(out) :: epsatm,xcccrc
1779  type(pawrad_type),intent(in) :: core_mesh,tproj_mesh,vale_mesh
1780  type(pawrad_type),intent(inout) ::pawrad,vloc_mesh
1781  type(pawtab_type),intent(inout) :: pawtab
1782 !arrays
1783  real(dp),intent(in) :: ncore(core_mesh%mesh_size),tncore(core_mesh%mesh_size)
1784  real(dp),intent(in) :: qgrid_vl(mqgrid_vl),qgrid_ff(mqgrid_ff)
1785  real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax)
1786  real(dp),intent(out) :: vlspl(mqgrid_vl,2)
1787  real(dp),intent(inout) :: tnvale(vale_mesh%mesh_size*pawtab%has_tvale)
1788  real(dp),intent(inout) :: tproj(tproj_mesh%mesh_size,pawtab%basis_size)
1789  real(dp),intent(inout) :: vlocr(vloc_mesh%mesh_size)
1790  type(pawrad_type),intent(in) :: radmesh(nmesh)
1791 
1792 !Local variables ------------------------------
1793 !scalars
1794  integer,parameter :: reduced_mshsz=2501
1795  integer :: ib,il,ilm,ilmn,iln,ir,isnotzero,itest
1796  integer :: j0lmn,jlm,jlmn,jln,klmn,msz,msz1,msz_tmp,mst_tmp,nspden
1797  logical :: has_dij0,non_magnetic_xc,reduced_ncor,reduced_nval,reduced_vloc,testval
1798  real(dp),parameter :: reduced_rstep=0.00025_dp,rm_vloc=20.0_dp
1799  real(dp) :: d2nvdq0,intg,intvh,lstep_tmp,qcore,qq,rstep_tmp,yp1,ypn
1800  character(len=500) :: msg
1801  type(pawang_type) :: pawang_tmp
1802  type(pawrad_type) :: rcore_mesh,rvale_mesh,rvloc_mesh,tproj_mesh_new
1803 !arrays
1804  real(dp) :: tmp_qgrid(1),tmp_q2vq(1)
1805  real(dp),allocatable :: ncorwk(:),nhat(:),nhatwk(:),nwk(:),r2k(:)
1806  real(dp),allocatable :: rtncor(:),rtnval(:),rvlocr(:)
1807  real(dp),allocatable :: vbare(:),vh(:),vhnzc(:)
1808  real(dp),allocatable :: vxc1(:),vxc2(:),work1(:),work1b(:),work2(:),work3(:),work4(:)
1809  logical :: tmp_lmselect(1)
1810 
1811 ! *************************************************************************
1812 
1813 !==========================================================
1814 !Perfom tests on meshes
1815 
1816 ! initialise logical
1817  non_magnetic_xc=.false.
1818 
1819 !Are radial meshes for Phi and Vloc compatibles ?
1820 ! if (vloc_mesh%rmax<pawrad%rmax) then
1821 !   write(msg, '(a,a,a)' )&
1822 !&   'Rmax for Vloc < Rmax !',ch10,&
1823 !&   'Action : check your pseudopotential (increase Vloc meshSize).'
1824 !   MSG_ERROR(msg)
1825 ! end if
1826 
1827 !Are mmax and mesh_size for partial waves compatibles ?
1828  if (mmax/=pawtab%partialwave_mesh_size) then
1829    write(msg, '(a,a,a)' )&
1830 &   'mmax /= phi_mesh_size in psp file !',ch10,&
1831 &   'Action: check your pseudopotential file.'
1832    MSG_ERROR(msg)
1833  end if
1834 
1835 !Are radial meshes for (t)Ncore and Phi compatibles ?
1836  if (core_mesh%mesh_size<pawtab%mesh_size) then
1837    write(msg, '(a,a,a,a,a)' )&
1838 &   'Mesh size for core density must be equal or larger',ch10,&
1839 &   'than mesh size for PAW augmentation regions !',ch10,&
1840 &   'Action : check your pseudopotential (increase Ncore meshSize).'
1841    MSG_ERROR(msg)
1842  end if
1843 
1844 !Are radial meshes for (t)Nvale and Phi compatibles ?
1845  if ((pawtab%has_tvale==1).and.(vale_mesh%rmax<pawrad%rmax)) then
1846    write(msg, '(a,a,a)' )&
1847 &   'Rmax for tNvale < Rmax for Phi !',ch10,&
1848 &   'Action : check your pseudopotential (increase tNvale meshSize).'
1849    MSG_ERROR(msg)
1850  end if
1851 
1852 !Is PAW radius included inside radial mesh ?
1853  if (pawtab%rpaw>pawrad%rmax+tol8) then
1854    write(msg, '(a,a,a)' )&
1855 &   'Radius of PAW sphere is outside the radial mesh !',ch10,&
1856 &   'Action: check your pseudopotential file.'
1857    MSG_ERROR(msg)
1858  end if
1859 
1860 !Max. radius of mesh for Vloc has to be "small" in order to avoid numeric noise ?
1861  if (vloc_mesh%rmax>rm_vloc) then
1862    msz_tmp=pawrad_ifromr(vloc_mesh,rm_vloc);mst_tmp=vloc_mesh%mesh_type
1863    rstep_tmp=vloc_mesh%rstep;lstep_tmp=vloc_mesh%lstep
1864    call pawrad_free(vloc_mesh)
1865    call pawrad_init(vloc_mesh,mesh_size=msz_tmp,mesh_type=mst_tmp,&
1866 &   rstep=rstep_tmp,lstep=lstep_tmp,r_for_intg=rm_vloc)
1867    write(msg, '(a,i4,a)' ) ' Mesh size for Vloc has been set to ', &
1868 &    vloc_mesh%mesh_size,' to avoid numerical noise.'
1869    call wrtout(std_out,msg,'COLL')
1870    call wrtout(ab_out,msg,'COLL')
1871  end if
1872 
1873 !This test has been disable... MT 2006-25-10
1874 !For Simpson rule, it is better to have odd mesh sizes
1875 !itest=0
1876 !do imsh=1,nmesh
1877 !if (mod(radmesh(imsh)%mesh_size,2)==0.and.radmesh(imsh)%mesh_type==1) itest=1
1878 !end do
1879 !if (itest==1) then
1880 !  write(msg, '(5a)' ) &
1881 !&   'Regular radial meshes should have odd number of points ',ch10,&
1882 !&   'for better accuracy of integration sheme (Simpson rule).',ch10,&
1883 !&   'Althought it''s not compulsory, you should change mesh sizes in psp file.'
1884 !  MSG_WARNING(msg)
1885 !end if
1886 
1887 !Test the compatibilty between Rpaw and mesh for (t)Phi
1888  if (pspversion>=3) then
1889    itest=pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw)
1890 !  This test has been disable... MT 2015-02-12
1891 !   if (itest+2>radmesh(imainmesh)%mesh_size) then
1892 !     write(msg, '(9a)' ) &
1893 !&     'Atomic data could produce inaccurate results:',ch10,&
1894 !&     'Wavefunctions and pseudo-wavefunctions should',ch10,&
1895 !&     'be given on a radial mesh larger than the PAW',ch10,&
1896 !&     'spheres (at least 2 additional points) !',ch10,&
1897 !&     'Action: check your pseudopotential file.'
1898 !     MSG_WARNING(msg)
1899 !   end if
1900    if (abs(pawtab%rpaw-radmesh(imainmesh)%rad(itest))<tol8) itest=itest-1
1901    ib=0;isnotzero=0
1902    do while ((isnotzero==0).and.(ib<pawtab%basis_size))
1903      ib=ib+1;ir=itest
1904      do while ((isnotzero==0).and.(ir<pawtab%mesh_size))
1905        ir=ir+1;if (abs(pawtab%phi(ir,ib)-pawtab%tphi(ir,ib))>tol8) isnotzero=1
1906      end do
1907    end do
1908    if (isnotzero>0) then
1909      write(msg, '(7a)' )&
1910 &     'Atomic data are inconsistent:',ch10,&
1911 &     'For r>=r_paw, pseudo wavefunctions are not',ch10,&
1912 &     'equal to wave functions (Phi(r)/=tPhi(r)) !',ch10,&
1913 &     'Action: check your pseudopotential file.'
1914      MSG_ERROR(msg)
1915    end if
1916  else
1917 !  For compatibility reasons set PAW radius at the end of mesh (older versions)
1918    if (pawtab%rpaw/=pawrad%rmax) then
1919      msz_tmp=pawrad%mesh_size;mst_tmp=pawrad%mesh_type
1920      rstep_tmp=pawrad%rstep;lstep_tmp=pawrad%lstep
1921      call pawrad_free(pawrad)
1922      call pawrad_init(pawrad,mesh_size=msz_tmp,mesh_type=mst_tmp,rstep=rstep_tmp,lstep=lstep_tmp)
1923      pawtab%rpaw=pawrad%rmax
1924    end if
1925  end if
1926 !If Vloc is a "Vbare" potential, it has to be localized inside PAW spheres
1927  if (vlocopt==0.and.(vloc_mesh%rmax>pawtab%rpaw+tol10)) then
1928    if(vlocr(pawrad_ifromr(vloc_mesh,pawtab%rpaw))>tol10) then
1929      write(msg, '(7a)' )&
1930 &     'Atomic data are inconsistent:',ch10,&
1931 &     'Local potential is a "Vbare" potential',ch10,&
1932 &     'and is not localized inside PAW sphere !',ch10,&
1933 &     'Vbare is set to zero if r>rpaw.'
1934      MSG_WARNING(msg)
1935      do ir=pawrad_ifromr(vloc_mesh,pawtab%rpaw),vloc_mesh%mesh_size
1936        vlocr(ir)=zero
1937      end do
1938    end if
1939  end if
1940 
1941 !==========================================================
1942 !Initializations
1943 
1944  has_dij0=(allocated(pawtab%dij0))
1945 
1946 !Allocate/initialize some dummy variables
1947  tmp_lmselect(1)=.true.
1948  if (pawxcdev==0) then
1949    pawang_tmp%l_size_max=1;pawang_tmp%angl_size=1;pawang_tmp%ylm_size=1
1950    pawang_tmp%use_ls_ylm=0;pawang_tmp%gnt_option=0;pawang_tmp%ngnt=0;pawang_tmp%nsym=0
1951    LIBPAW_ALLOCATE(pawang_tmp%angwgth,(1))
1952    pawang_tmp%angwgth(1)=one
1953    LIBPAW_ALLOCATE(pawang_tmp%anginit,(3,1))
1954    pawang_tmp%anginit(1,1)=one
1955    pawang_tmp%anginit(2:3,1)=zero
1956    LIBPAW_ALLOCATE(pawang_tmp%ylmr,(1,1))
1957    pawang_tmp%ylmr(1,1)=1._dp/sqrt(four_pi)
1958    LIBPAW_ALLOCATE(pawang_tmp%ylmrgr,(3,1,1))
1959    pawang_tmp%ylmrgr(1:3,1,1)=zero
1960  end if
1961 
1962 !==========================================================
1963 !Compute ffspl(q) (and derivatives)
1964 
1965  ffspl=zero
1966  if (mqgrid_ff>0) then
1967    call pawpsp_nl(ffspl,pawtab%indlmn,pawtab%lmn_size,lnmax,mqgrid_ff,qgrid_ff,&
1968 &                 tproj_mesh,tproj)
1969  end if
1970 
1971 !==========================================================
1972 !Compute eventually compensation charge radius (i.e. radius for shape functions)
1973 
1974  if (pawtab%shape_type>0.and.pawtab%rshp<1.d-8) then
1975    pawtab%rshp=pawtab%rpaw
1976  else if (pawtab%shape_type==-1) then
1977    ir=pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw)+1;isnotzero=0
1978    do while ((isnotzero==0).and.(ir>1))
1979      ir=ir-1;il=0
1980      do while ((isnotzero==0).and.(il<pawtab%l_size))
1981        il=il+1;if (pawtab%shapefunc(ir,il)>tol16) isnotzero=1
1982      end do
1983    end do
1984    ir=min(ir+1,pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw))
1985    pawtab%rshp=radmesh(imainmesh)%rad(ir)
1986    do il=1,pawtab%l_size
1987      if (pawtab%shapefunc(ir,il)>tol6) then
1988        write(msg, '(a,a,a)' )&
1989 &       'Shape function is not zero at PAW radius !',ch10,&
1990 &       'Action: check your pseudopotential file.'
1991        MSG_ERROR(msg)
1992      end if
1993    end do
1994  end if
1995 
1996 !==========================================================
1997 !Compute compensation charge density (nhat)
1998 !Add it to pseudo valence density
1999 
2000  if (pawtab%has_tvale==1) then
2001    msz=vale_mesh%mesh_size
2002    LIBPAW_ALLOCATE(nhat,(msz))
2003 !  A-Has to compute norm of nhat (Int[n-tild_n])
2004    testval=(abs(tnvale(msz))<tol9)
2005 !  A1-If tnvale is not given with enough points,
2006 !  try to compute it from rhoij0 and tphi
2007    if (.not.testval) then
2008      msz1=pawtab%mesh_size
2009 !    Compute n and tild_n from phi and tphi
2010      LIBPAW_ALLOCATE(work1,(msz1))
2011      LIBPAW_ALLOCATE(work2,(msz1))
2012      work1=zero
2013      work2=zero
2014      do jlmn=1,pawtab%lmn_size
2015        j0lmn=jlmn*(jlmn-1)/2;jln=pawtab%indlmn(5,jlmn)
2016        do ilmn=1,jlmn
2017          klmn=j0lmn+ilmn;iln=pawtab%indlmn(5,ilmn)
2018          yp1=two;if (ilmn==jlmn) yp1=one
2019          work1(1:msz1)=work1(1:msz1)+yp1*pawtab%rhoij0(klmn) &
2020 &         *pawtab% phi(1:msz1,iln)*pawtab% phi(1:msz1,jln)
2021          work2(1:msz1)=work2(1:msz1)+yp1*pawtab%rhoij0(klmn) &
2022 &         *pawtab%tphi(1:msz1,iln)*pawtab%tphi(1:msz1,jln)
2023        end do
2024      end do
2025 !    Spline tnvale onto pawrad if needed
2026      LIBPAW_ALLOCATE(nwk,(msz1))
2027      if ((vale_mesh%mesh_type/=pawrad%mesh_type).or.(vale_mesh%rstep/=pawrad%rstep).or.&
2028 &     (vale_mesh%lstep/=pawrad%lstep)) then
2029        LIBPAW_ALLOCATE(work3,(vale_mesh%mesh_size))
2030        call bound_deriv(tnvale(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn)
2031        call paw_spline(vale_mesh%rad,tnvale,vale_mesh%mesh_size,yp1,ypn,work3)
2032        call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,tnvale,work3,msz1,pawrad%rad(1:msz1),nwk(1:msz1))
2033        LIBPAW_DEALLOCATE(work3)
2034      else
2035        nwk(1:msz1)=tnvale(1:msz1)
2036      end if
2037 !    Compare tild_n and tnvale (inside aug. region)
2038      if (maxval(abs((nwk(1:msz1)*four_pi*pawrad%rad(1:msz1)**2)-work2(1:msz1)))<tol6) then
2039 !      If equality then compute Int[n-tild_n]
2040        work1=work1-work2
2041        call simp_gen(qq,work1,pawrad)
2042        qq=qq/four_pi
2043      else
2044 !      If not equality, will use tnvale
2045        testval=.true.
2046        write(msg, '(3a)' ) &
2047 &       'Valence density is not given with enough points',ch10,&
2048 &       'in psp file. Some charge estimations will be coarse.'
2049        MSG_WARNING(msg)
2050      end if
2051      LIBPAW_DEALLOCATE(nwk)
2052      LIBPAW_DEALLOCATE(work1)
2053      LIBPAW_DEALLOCATE(work2)
2054    end if
2055 !  A2-If tnvale is given with enough points, use it
2056    if (testval) then
2057      nhat(1:msz)=tnvale(1:msz)*vale_mesh%rad(1:msz)**2
2058      call simp_gen(qq,nhat,vale_mesh)
2059      qq=zion/four_pi-qq
2060    end if
2061 !  B-Compute nhat and add it to pseudo valence density
2062    call atompaw_shpfun(0,vale_mesh,intg,pawtab,nhat)
2063    nhat(1:msz)=qq*nhat(1:msz)
2064    tnvale(1:msz)=tnvale(1:msz)+nhat(1:msz)
2065  end if
2066 
2067 !==========================================================
2068 !If Vloc potential is in "Vbare" format, translate it into VH(tnzc) format
2069 
2070  if (vlocopt==0) then
2071    write(msg,'(a)') ' Local potential is in "Vbare" format... '
2072    call wrtout(ab_out,msg,'COLL')
2073    call wrtout(std_out,  msg,'COLL')
2074    msz=core_mesh%mesh_size
2075    LIBPAW_ALLOCATE(r2k,(msz))
2076    call atompaw_shpfun(0,core_mesh,intg,pawtab,r2k)
2077    r2k(1:msz)=r2k(1:msz)*core_mesh%rad(1:msz)**2
2078 !  Compute VH[4pi.r2.n(r)=4pi.r2.tncore(r)+(Qcore-Z).r2.k(r)]
2079    LIBPAW_ALLOCATE(nwk,(core_mesh%mesh_size))
2080    LIBPAW_ALLOCATE(vh,(core_mesh%mesh_size))
2081    if (core_mesh%mesh_type==5) then
2082      nwk(1:msz)=tncore(1:msz)*four_pi*core_mesh%rad(1:msz)**2
2083      call simp_gen(qcore,nwk,core_mesh)
2084      qcore=znucl-zion-qcore
2085    else
2086      nwk(1:msz)=(ncore(1:msz)-tncore(1:msz))*four_pi*core_mesh%rad(1:msz)**2
2087      ib=1
2088      do ir=msz,2,-1
2089        if(abs(nwk(ir))<tol14)ib=ir
2090      end do
2091      call simp_gen(qcore,nwk,core_mesh,r_for_intg=core_mesh%rad(ib))
2092      nwk(1:msz)=tncore(1:msz)*four_pi*core_mesh%rad(1:msz)**2
2093    end if
2094    nwk(1:msz)=nwk(1:msz)+r2k(1:msz)*(qcore-znucl)
2095    call poisson(nwk,0,core_mesh,vh)
2096    vh(2:msz)=vh(2:msz)/core_mesh%rad(2:msz)
2097    call pawrad_deducer0(vh,msz,core_mesh)
2098 
2099    LIBPAW_DEALLOCATE(nwk)
2100 !  Eventually spline Vbare
2101    LIBPAW_ALLOCATE(vbare,(core_mesh%mesh_size))
2102    if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.&
2103 &   (core_mesh%rstep    /=vloc_mesh%rstep)    .or.&
2104 &   (core_mesh%lstep    /=vloc_mesh%lstep)) then
2105      msz=core_mesh%mesh_size;if (vloc_mesh%rmax<core_mesh%rmax) msz=pawrad_ifromr(core_mesh,vloc_mesh%rmax)
2106      call bound_deriv(vlocr(1:vloc_mesh%mesh_size),vloc_mesh,vloc_mesh%mesh_size,yp1,ypn)
2107      LIBPAW_ALLOCATE(work1,(vloc_mesh%mesh_size))
2108      LIBPAW_ALLOCATE(work2,(vloc_mesh%mesh_size))
2109      call paw_spline(vloc_mesh%rad,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1)
2110      call paw_splint(vloc_mesh%mesh_size,vloc_mesh%rad,vlocr,work1,msz,core_mesh%rad(1:msz),vbare)
2111      LIBPAW_DEALLOCATE(work1)
2112      LIBPAW_DEALLOCATE(work2)
2113    else
2114      msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size)
2115      vbare(1:msz)=vlocr(1:msz)
2116    end if
2117 !  Build VH(tnzc) from Vbare
2118    vlocr(1:msz)=vbare(1:msz)+vh(1:msz)
2119    if(vloc_mesh%mesh_size>msz)then
2120      vlocr(msz+1:vloc_mesh%mesh_size)=vh(msz)*vloc_mesh%rad(msz)/vloc_mesh%rad(msz+1:vloc_mesh%mesh_size)
2121    end if
2122    LIBPAW_DEALLOCATE(vbare)
2123    LIBPAW_DEALLOCATE(vh)
2124 
2125 !  Compute <tPhi_i|VH(tnzc)|tPhi_j> and int[VH(tnzc)*Qijhat(r)dr] parts of Dij0
2126 !  Note: it is possible as core_mesh and radmesh(imainmesh) have the same steps
2127    if (has_dij0) then
2128      msz=radmesh(imainmesh)%mesh_size
2129      LIBPAW_ALLOCATE(work1,(msz))
2130      work1(1:msz)=vlocr(1:msz)*r2k(1:msz)
2131      call simp_gen(intvh,work1,radmesh(imainmesh))
2132      do jlmn=1,pawtab%lmn_size
2133        j0lmn=jlmn*(jlmn-1)/2;jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn)
2134        do ilmn=1,jlmn
2135          klmn=j0lmn+ilmn;ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn)
2136          if (jlm==ilm) then
2137            work1(1:msz)=pawtab%tphi(1:msz,iln)*pawtab%tphi(1:msz,jln)*(vlocr(1:msz)-intvh) &
2138 &           -pawtab%phi (1:msz,iln)*pawtab%phi (1:msz,jln)*intvh
2139            call simp_gen(intg,work1,radmesh(imainmesh))
2140            pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
2141          end if
2142        end do
2143      end do
2144      LIBPAW_DEALLOCATE(work1)
2145    end if
2146    LIBPAW_DEALLOCATE(r2k)
2147  end if
2148 
2149 !==========================================================
2150 !If usexcnhat in psp file is different from usexcnhat chosen
2151 !by user, convert VH(tnzc) and Dij0
2152 
2153  if (pawtab%usexcnhat==-1) then
2154    pawtab%usexcnhat=usexcnhat
2155  else if (usexcnhat/=pawtab%usexcnhat) then
2156    if (pawtab%has_tvale==0) then
2157      write(msg, '(5a)' ) &
2158 &     'It is only possible to modify the use of compensation charge density',ch10,&
2159 &     'for a file format greater or equal than paw4 !',ch10,&
2160 &     'Action: use usexcnhat=-1 in input file or change psp file format.'
2161      MSG_WARNING(msg)
2162    else
2163      msz=vloc_mesh%mesh_size
2164 !    Retrieve tvale and nhat onto vloc mesh
2165      LIBPAW_ALLOCATE(nwk,(msz))
2166      LIBPAW_ALLOCATE(ncorwk,(msz))
2167      LIBPAW_ALLOCATE(nhatwk,(msz))
2168      nwk=zero;ncorwk=zero;nhatwk=zero
2169      if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.&
2170 &        (core_mesh%rstep    /=vloc_mesh%rstep)    .or.&
2171 &        (core_mesh%lstep    /=vloc_mesh%lstep)) then
2172        LIBPAW_ALLOCATE(work1,(core_mesh%mesh_size))
2173        msz1=msz;if (core_mesh%rmax<vloc_mesh%rmax) msz1=pawrad_ifromr(vloc_mesh,core_mesh%rmax)
2174        call bound_deriv(tncore(1:core_mesh%mesh_size),core_mesh,core_mesh%mesh_size,yp1,ypn)
2175        call paw_spline(core_mesh%rad,tncore,core_mesh%mesh_size,yp1,ypn,work1)
2176        call paw_splint(core_mesh%mesh_size,core_mesh%rad,tncore,work1,msz1,vloc_mesh%rad(1:msz1),ncorwk(1:msz1))
2177        LIBPAW_DEALLOCATE(work1)
2178      else
2179        msz1=min(core_mesh%mesh_size,msz)
2180        ncorwk(1:msz1)=tncore(1:msz1)
2181      end if
2182      if ((vale_mesh%mesh_type/=vloc_mesh%mesh_type).or.&
2183 &        (vale_mesh%rstep    /=vloc_mesh%rstep)    .or.&
2184 &        (vale_mesh%lstep    /=vloc_mesh%lstep)) then
2185        LIBPAW_ALLOCATE(work1,(vale_mesh%mesh_size))
2186        msz1=msz;if (vale_mesh%rmax<vloc_mesh%rmax) msz1=pawrad_ifromr(vloc_mesh,vale_mesh%rmax)
2187        call bound_deriv(tnvale(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn)
2188        call paw_spline(vale_mesh%rad,tnvale,vale_mesh%mesh_size,yp1,ypn,work1)
2189        call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,tnvale,work1,msz1,vloc_mesh%rad(1:msz1),nwk(1:msz1))
2190        call bound_deriv(nhat(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn)
2191        call paw_spline(vale_mesh%rad,nhat,vale_mesh%mesh_size,yp1,ypn,work1)
2192        call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,nhat,work1,msz1,vloc_mesh%rad(1:msz1),nhatwk(1:msz1))
2193        LIBPAW_DEALLOCATE(work1)
2194      else
2195        msz1=min(vale_mesh%mesh_size,msz)
2196        nwk   (1:msz1)=tnvale(1:msz1)
2197        nhatwk(1:msz1)=nhat  (1:msz1)
2198      end if
2199      nwk=nwk-nhatwk
2200      nwk=sqrt(four_pi)*nwk;nhatwk=sqrt(four_pi)*nhatwk ! 0th-order moment of densities
2201 
2202 !    Compute Vxc without nhat (vxc1) and with nhat (vxc2)
2203      nspden=1
2204 #if defined LIBPAW_HAVE_LIBXC
2205      if (ixc<0) nspden=libxc_functionals_nspin()
2206 #endif
2207      if (ixc<0) then
2208        LIBPAW_ALLOCATE(vxc1,(msz*nspden))
2209        LIBPAW_ALLOCATE(vxc2,(msz*nspden))
2210        LIBPAW_ALLOCATE(work1,(msz))
2211        LIBPAW_ALLOCATE(work1b,(msz))
2212        LIBPAW_ALLOCATE(work2,(msz*nspden))
2213        LIBPAW_ALLOCATE(work3,(msz*nspden))
2214        work2(1:msz)=nwk
2215        work3(1:msz)=nhatwk
2216        if (nspden==2) work2(msz+1:2*msz)=half*nwk
2217        if (nspden==2) work3(msz+1:2*msz)=half*nhatwk
2218        if (pawxcdev/=0) then
2219          call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,work3,0,non_magnetic_xc,msz,nspden,5,&
2220 &         pawang_tmp,vloc_mesh,pawxcdev,work2,pawtab%usetcore,0,vxc1,xclevel,xc_denpos)
2221          call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,work3,0,non_magnetic_xc,msz,nspden,5,&
2222 &         pawang_tmp,vloc_mesh,pawxcdev,work2,pawtab%usetcore,2,vxc2,xclevel,xc_denpos)
2223          vxc1=vxc1/sqrt(four_pi);vxc2=vxc2/sqrt(four_pi) ! Deduce Vxc from its first moment
2224        else
2225          call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,work3,0,0,non_magnetic_xc,msz,nspden,5,&
2226 &         pawang_tmp,vloc_mesh,work2,pawtab%usetcore,0,vxc1,xclevel,xc_denpos)
2227          call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,work3,0,0,non_magnetic_xc,msz,nspden,5,&
2228 &         pawang_tmp,vloc_mesh,work2,pawtab%usetcore,2,vxc2,xclevel,xc_denpos)
2229        end if
2230        LIBPAW_DEALLOCATE(nwk)
2231        LIBPAW_DEALLOCATE(ncorwk)
2232        LIBPAW_DEALLOCATE(nhatwk)
2233        LIBPAW_DEALLOCATE(work1)
2234        LIBPAW_DEALLOCATE(work1b)
2235        LIBPAW_DEALLOCATE(work2)
2236        LIBPAW_DEALLOCATE(work3)
2237      else
2238        LIBPAW_ALLOCATE(vxc1,(msz))
2239        LIBPAW_ALLOCATE(vxc2,(msz))
2240        LIBPAW_ALLOCATE(work1,(msz))
2241        LIBPAW_ALLOCATE(work1b,(msz))
2242        if (pawxcdev/=0) then
2243          call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,nhatwk,0,non_magnetic_xc,msz,1,5,&
2244 &         pawang_tmp,vloc_mesh,pawxcdev,nwk,pawtab%usetcore,0,vxc1,xclevel,xc_denpos)
2245          call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,nhatwk,0,non_magnetic_xc,msz,1,5,&
2246 &         pawang_tmp,vloc_mesh,pawxcdev,nwk,pawtab%usetcore,2,vxc2,xclevel,xc_denpos)
2247          vxc1=vxc1/sqrt(four_pi);vxc2=vxc2/sqrt(four_pi) ! Deduce Vxc from its first moment
2248        else
2249          call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,nhatwk,0,0,non_magnetic_xc,msz,1,5,&
2250 &         pawang_tmp,vloc_mesh,nwk,pawtab%usetcore,0,vxc1,xclevel,xc_denpos)
2251          call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,nhatwk,0,0,non_magnetic_xc,msz,1,5,&
2252 &         pawang_tmp,vloc_mesh,nwk,pawtab%usetcore,2,vxc2,xclevel,xc_denpos)
2253        end if
2254        LIBPAW_DEALLOCATE(nwk)
2255        LIBPAW_DEALLOCATE(ncorwk)
2256        LIBPAW_DEALLOCATE(nhatwk)
2257        LIBPAW_DEALLOCATE(work1)
2258        LIBPAW_DEALLOCATE(work1b)
2259      endif
2260 !    Compute difference of XC potentials
2261      if (usexcnhat==0.and.pawtab%usexcnhat/=0)  vxc1(1:msz)=vxc2(1:msz)-vxc1(1:msz)
2262      if (usexcnhat/=0.and.pawtab%usexcnhat==0)  vxc1(1:msz)=vxc1(1:msz)-vxc2(1:msz)
2263 !    Modify VH(tnzc)
2264      vlocr(1:msz)=vlocr(1:msz)-vxc1(1:msz)
2265      if (has_dij0) then
2266 !      Modify  Dij0
2267        LIBPAW_ALLOCATE(work2,(pawtab%lmn2_size))
2268        call atompaw_kij(pawtab%indlmn,work2,pawtab%lmn_size,ncore,0,0,pawtab,pawrad,&
2269 &                       core_mesh,vloc_mesh,vxc1(1:msz),znucl)
2270        pawtab%dij0=work2
2271        LIBPAW_DEALLOCATE(work2)
2272      end if
2273      LIBPAW_DEALLOCATE(vxc1)
2274      LIBPAW_DEALLOCATE(vxc2)
2275    end if ! has_tvale/=0
2276  end if
2277  if (pawtab%usexcnhat==0) then
2278    write(msg,'(a)') &
2279 &   ' Compensation charge density is not taken into account in XC energy/potential'
2280    call wrtout(ab_out,msg,'COLL')
2281    call wrtout(std_out,  msg,'COLL')
2282  end if
2283  if (pawtab%usexcnhat==1) then
2284    write(msg,'(a)') &
2285 &   ' Compensation charge density is taken into account in XC energy/potential'
2286    call wrtout(ab_out,msg,'COLL')
2287    call wrtout(std_out,  msg,'COLL')
2288  end if
2289 
2290 !==========================================================
2291 ! Calculate the coefficient beta = \int { vH[nZc](r) - vloc(r) } 4pi r^2 dr
2292 !
2293  LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size))
2294  LIBPAW_ALLOCATE(nwk,(core_mesh%mesh_size))
2295 ! get vH[nZc]
2296  call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl)
2297 
2298 !Transpose vlocr mesh into core mesh
2299  nwk(:)=zero
2300  if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.&
2301 & (core_mesh%rstep    /=vloc_mesh%rstep)    .or.&
2302 & (core_mesh%lstep    /=vloc_mesh%lstep)) then
2303    msz=core_mesh%mesh_size;if (vloc_mesh%rmax<core_mesh%rmax) msz=pawrad_ifromr(core_mesh,vloc_mesh%rmax)
2304    call bound_deriv(vlocr(1:vloc_mesh%mesh_size),vloc_mesh,vloc_mesh%mesh_size,yp1,ypn)
2305    LIBPAW_ALLOCATE(work1,(vloc_mesh%mesh_size))
2306    LIBPAW_ALLOCATE(work2,(vloc_mesh%mesh_size))
2307    call paw_spline(vloc_mesh%rad,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1)
2308    call paw_splint(vloc_mesh%mesh_size,vloc_mesh%rad,vlocr,work1,msz,core_mesh%rad(1:msz),nwk)
2309    LIBPAW_DEALLOCATE(work1)
2310    LIBPAW_DEALLOCATE(work2)
2311  else
2312    msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size)
2313    nwk(1:msz)=vlocr(1:msz)
2314  end if
2315 
2316 !Difference
2317  nwk(1:msz)=vhnzc(1:msz)-nwk(1:msz)
2318  if (msz<core_mesh%mesh_size) nwk(msz+1:core_mesh%mesh_size)=zero
2319 
2320 !Perform the spherical integration
2321  nwk(1:msz)=nwk(1:msz)*four_pi*core_mesh%rad(1:msz)**2
2322 
2323  call simp_gen(pawtab%beta,nwk,core_mesh)
2324 
2325  LIBPAW_DEALLOCATE(vhnzc)
2326  LIBPAW_DEALLOCATE(nwk)
2327 
2328  write(msg,'(a,e18.6)') &
2329 &  ' beta integral value: ',pawtab%beta
2330  call wrtout(std_out,msg,'COLL')
2331 
2332 
2333 !==========================================================
2334 !Try to optimize CPU time:
2335 !If Vloc mesh size is big, spline Vloc into a smaller log. mesh
2336 
2337  reduced_vloc=(vloc_mesh%mesh_size>int(reduced_mshsz))
2338  if (reduced_vloc) then
2339    msz=vloc_mesh%mesh_size
2340    lstep_tmp=log(0.9999999_dp*vloc_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2)
2341    call pawrad_init(rvloc_mesh,mesh_size=reduced_mshsz,mesh_type=3,&
2342 &   rstep=reduced_rstep,lstep=lstep_tmp)
2343    LIBPAW_ALLOCATE(rvlocr,(reduced_mshsz))
2344    call bound_deriv(vlocr(1:msz),vloc_mesh,msz,yp1,ypn)
2345    LIBPAW_ALLOCATE(work1,(msz))
2346    LIBPAW_ALLOCATE(work2,(msz))
2347    LIBPAW_ALLOCATE(work3,(msz))
2348    work3(1:msz)=vloc_mesh%rad(1:msz)
2349    call paw_spline(work3,vlocr,msz,yp1,ypn,work1)
2350    call paw_splint(msz,work3,vlocr,work1,reduced_mshsz,rvloc_mesh%rad,rvlocr)
2351    LIBPAW_DEALLOCATE(work1)
2352    LIBPAW_DEALLOCATE(work2)
2353    LIBPAW_DEALLOCATE(work3)
2354  end if
2355 
2356 !Keep VH(tnZc) eventually in memory
2357  if (pawtab%has_vhtnzc==1) then
2358    LIBPAW_ALLOCATE(pawtab%vhtnzc,(pawtab%mesh_size))
2359    if ((reduced_vloc).and.(rvloc_mesh%mesh_type==pawrad%mesh_type)&
2360 &   .and.(rvloc_mesh%rstep==pawrad%rstep).and.(rvloc_mesh%lstep==pawrad%lstep)) then
2361      pawtab%vhtnzc(1:pawtab%mesh_size)=rvlocr(1:pawtab%mesh_size)
2362      pawtab%has_vhtnzc=2
2363    else if ((vloc_mesh%mesh_type==pawrad%mesh_type)&
2364 &     .and.(vloc_mesh%rstep==pawrad%rstep).and.(vloc_mesh%lstep==pawrad%lstep)) then
2365      pawtab%vhtnzc(1:pawtab%mesh_size)=vlocr(1:pawtab%mesh_size)
2366      pawtab%has_vhtnzc=2
2367    else
2368      msg = 'Vloc mesh is not right !'
2369      MSG_ERROR(msg)
2370    end if
2371  end if
2372 
2373 !==========================================================
2374 !Try to optimize CPU time:
2375 !If ncore mesh size is big, spline tncore into a smaller log. mesh
2376 
2377  reduced_ncor=(core_mesh%mesh_size>int(reduced_mshsz)).and.(pawtab%usetcore/=0)
2378  if (reduced_ncor) then
2379    msz=core_mesh%mesh_size
2380    lstep_tmp=log(0.9999999_dp*core_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2)
2381    call pawrad_init(rcore_mesh,mesh_size=reduced_mshsz,mesh_type=3,&
2382 &   rstep=reduced_rstep,lstep=lstep_tmp)
2383    LIBPAW_ALLOCATE(rtncor,(reduced_mshsz))
2384    call bound_deriv(tncore(1:msz),core_mesh,msz,yp1,ypn)
2385    LIBPAW_ALLOCATE(work1,(msz))
2386    LIBPAW_ALLOCATE(work2,(msz))
2387    LIBPAW_ALLOCATE(work3,(msz))
2388    work3(1:msz)=core_mesh%rad(1:msz)
2389    call paw_spline(work3,tncore,msz,yp1,ypn,work1)
2390    call paw_splint(msz,work3,tncore,work1,reduced_mshsz,rcore_mesh%rad,rtncor)
2391    LIBPAW_DEALLOCATE(work1)
2392    LIBPAW_DEALLOCATE(work2)
2393    LIBPAW_DEALLOCATE(work3)
2394  end if
2395 
2396 !==========================================================
2397 !Try to optimize CPU time:
2398 !If vale mesh size is big, spline tnvale into a smaller log. mesh
2399 
2400  if (pawtab%has_tvale==1) then
2401    reduced_nval=(vale_mesh%mesh_size>int(reduced_mshsz))
2402    if (reduced_nval) then
2403      msz=vale_mesh%mesh_size
2404      lstep_tmp=log(0.9999999_dp*vale_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2)
2405      call pawrad_init(rvale_mesh,mesh_size=reduced_mshsz,mesh_type=3,&
2406 &     rstep=reduced_rstep,lstep=lstep_tmp)
2407      LIBPAW_ALLOCATE(rtnval,(reduced_mshsz))
2408      call bound_deriv(tnvale(1:msz),vale_mesh,msz,yp1,ypn)
2409      LIBPAW_ALLOCATE(work1,(msz))
2410      LIBPAW_ALLOCATE(work2,(msz))
2411      LIBPAW_ALLOCATE(work3,(msz))
2412      work3(1:msz)=vale_mesh%rad(1:msz)
2413      call paw_spline(work3,tnvale,msz,yp1,ypn,work1)
2414      call paw_splint(msz,work3,tnvale,work1,reduced_mshsz,rvale_mesh%rad,rtnval)
2415      LIBPAW_DEALLOCATE(work1)
2416      LIBPAW_DEALLOCATE(work2)
2417      LIBPAW_DEALLOCATE(work3)
2418    end if
2419  else
2420    reduced_nval=.false.
2421  end if
2422 
2423 !==========================================================
2424 !Compute Vlspl(q) (and second derivative) from Vloc(r)
2425 
2426 !Compute Vlspl(q)=q^2.Vloc(q) from vloc(r)
2427  if(mqgrid_vl>0) then
2428    if (reduced_vloc) then
2429      call pawpsp_lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),rvloc_mesh,rvlocr,yp1,ypn,zion)
2430    else
2431      call pawpsp_lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),vloc_mesh,vlocr,yp1,ypn,zion)
2432    end if
2433 !  Compute second derivative of Vlspl(q)
2434    call paw_spline(qgrid_vl,vlspl(:,1),mqgrid_vl,yp1,ypn,vlspl(:,2))
2435  else
2436    ! Only to compute epsatm
2437    epsatm=zero
2438    if (reduced_vloc) then
2439      call pawpsp_lo(epsatm,1,tmp_qgrid,tmp_q2vq,rvloc_mesh,rvlocr,yp1,ypn,zion)
2440    else
2441      call pawpsp_lo(epsatm,1,tmp_qgrid,tmp_q2vq,vloc_mesh,vlocr,yp1,ypn,zion)
2442    end if
2443  end if
2444 
2445 !==========================================================
2446 !Compute tcorespl(q) (and second derivative) from tNcore(r)
2447 
2448  pawtab%mqgrid=mqgrid_vl
2449  xcccrc=core_mesh%rmax
2450  LIBPAW_ALLOCATE(pawtab%tcorespl,(pawtab%mqgrid,2))
2451 
2452  if(mqgrid_vl>0.and.pawtab%usetcore/=0) then
2453 !  Compute tcorespl(q)=tNc(q) from tNcore(r)
2454    if (reduced_ncor) then
2455      call pawpsp_cg(pawtab%dncdq0,pawtab%d2ncdq0,mqgrid_vl,qgrid_vl,pawtab%tcorespl(:,1),rcore_mesh,rtncor,yp1,ypn)
2456    else
2457    call pawpsp_cg(pawtab%dncdq0,pawtab%d2ncdq0,mqgrid_vl,qgrid_vl,pawtab%tcorespl(:,1),core_mesh,tncore,yp1,ypn)
2458    end if
2459 !  Compute second derivative of tcorespl(q)
2460    call paw_spline(qgrid_vl,pawtab%tcorespl(:,1),mqgrid_vl,yp1,ypn,pawtab%tcorespl(:,2))
2461  else
2462    pawtab%tcorespl=zero
2463    pawtab%dncdq0=zero
2464    pawtab%d2ncdq0=zero
2465  end if
2466 
2467 !==========================================================
2468 !Compute tvalespl(q) (and second derivative) from tNvale(r)
2469 
2470  if (pawtab%has_tvale/=0.and.mqgrid_vl>0) then
2471    LIBPAW_ALLOCATE(pawtab%tvalespl,(pawtab%mqgrid,2))
2472    if (reduced_nval) then
2473      call pawpsp_cg(pawtab%dnvdq0,d2nvdq0,mqgrid_vl,qgrid_vl,pawtab%tvalespl(:,1),rvale_mesh,rtnval,yp1,ypn)
2474      pawtab%tnvale_mesh_size=rvale_mesh%mesh_size
2475    else
2476      call pawpsp_cg(pawtab%dnvdq0,d2nvdq0,mqgrid_vl,qgrid_vl,pawtab%tvalespl(:,1),vale_mesh,tnvale,yp1,ypn)
2477      pawtab%tnvale_mesh_size=vale_mesh%mesh_size
2478    end if
2479 !  Compute second derivative of tvalespl(q)
2480    call paw_spline(qgrid_vl,pawtab%tvalespl(:,1),mqgrid_vl,yp1,ypn,pawtab%tvalespl(:,2))
2481  else
2482    pawtab%dnvdq0=zero
2483    pawtab%tnvale_mesh_size=0
2484  end if
2485 
2486 !==================================================
2487 !Compute Ex-correlation energy for the core density
2488 
2489  nspden=1
2490 #if defined LIBPAW_HAVE_LIBXC
2491  if (ixc<0) nspden=libxc_functionals_nspin()
2492 #endif
2493 
2494  LIBPAW_ALLOCATE(work1,(core_mesh%mesh_size*nspden))
2495  LIBPAW_ALLOCATE(work1b,(core_mesh%mesh_size*nspden))
2496  LIBPAW_ALLOCATE(work2,(core_mesh%mesh_size*nspden))
2497  LIBPAW_ALLOCATE(work3,(1))
2498  LIBPAW_ALLOCATE(work4,(core_mesh%mesh_size))
2499  work1(:)=zero
2500  if (pawxcdev/=0) then
2501    call pawxcm(ncore,pawtab%exccore,yp1,0,ixc,work4,1,tmp_lmselect,work3,0,non_magnetic_xc,core_mesh%mesh_size,&
2502 &   nspden,4,pawang_tmp,core_mesh,pawxcdev,work1,1,0,work2,xclevel,xc_denpos)
2503  else
2504    call pawxc(ncore,pawtab%exccore,yp1,ixc,work4,work1b,1,tmp_lmselect,work3,0,0,non_magnetic_xc,core_mesh%mesh_size,&
2505 &   nspden,4,pawang_tmp,core_mesh,work1,1,0,work2,xclevel,xc_denpos)
2506  end if
2507  LIBPAW_DEALLOCATE(work1)
2508  LIBPAW_DEALLOCATE(work1b)
2509  LIBPAW_DEALLOCATE(work2)
2510  LIBPAW_DEALLOCATE(work3)
2511  LIBPAW_DEALLOCATE(work4)
2512 
2513 !==================================================
2514 !Compute atomic contribution to Dij (Dij0)
2515 !if not already in memory
2516  if ((.not.has_dij0).and.(pawtab%has_kij==2.or.pawtab%has_kij==-1)) then
2517    LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size))
2518    if (reduced_vloc) then
2519      call atompaw_dij0(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,&
2520 &                      rvloc_mesh,rvlocr,znucl)
2521    else
2522      call atompaw_dij0(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,&
2523 &                      vloc_mesh,vlocr,znucl)
2524    end if
2525    has_dij0=.true.
2526  end if
2527 !==================================================
2528 !Compute kinetic operator contribution to Dij
2529 
2530  if (pawtab%has_kij==1.and.has_dij0) then
2531    LIBPAW_ALLOCATE(pawtab%kij,(pawtab%lmn2_size))
2532    call atompaw_kij(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,1,pawtab,pawrad,core_mesh,&
2533 &                   vloc_mesh,vlocr,znucl)
2534    pawtab%has_kij=2
2535  end if
2536 
2537 !pawtab%has_kij=-1 means that kij does not have to be kept in memory
2538  if (pawtab%has_kij==-1) then
2539    LIBPAW_DEALLOCATE(pawtab%kij)
2540    pawtab%has_kij=0
2541  end if
2542 
2543 !==========================================================
2544 !If projectors have to be kept in memory, we need
2545 !them on the main radial mesh (so, spline them if necessary)
2546 
2547  if (pawtab%has_tproj>0) then
2548    if ((tproj_mesh%mesh_type/=pawrad%mesh_type).or.&
2549 &      (tproj_mesh%rstep    /=pawrad%rstep).or.&
2550 &      (tproj_mesh%lstep    /=pawrad%lstep)) then
2551      ir=pawrad_ifromr(pawrad,tproj_mesh%rmax)
2552      call pawrad_init(tproj_mesh_new,mesh_size=ir,mesh_type=pawrad%mesh_type,&
2553 &                     rstep=pawrad%rstep,lstep=pawrad%lstep)
2554      LIBPAW_ALLOCATE(pawtab%tproj,(tproj_mesh_new%mesh_size,pawtab%basis_size))
2555      LIBPAW_ALLOCATE(work1,(tproj_mesh%mesh_size))
2556      do ib=1,pawtab%basis_size
2557        call bound_deriv(tproj(:,ib),tproj_mesh,tproj_mesh%mesh_size,yp1,ypn)
2558        call paw_spline(tproj_mesh%rad,tproj(:,ib),tproj_mesh%mesh_size,yp1,ypn,work1)
2559        call paw_splint(tproj_mesh%mesh_size,tproj_mesh%rad,tproj(:,ib),work1,&
2560 &           tproj_mesh_new%mesh_size,tproj_mesh_new%rad,pawtab%tproj(:,ib))
2561      end do
2562      LIBPAW_DEALLOCATE(work1)
2563      call pawrad_free(tproj_mesh_new)
2564    else
2565      LIBPAW_ALLOCATE(pawtab%tproj,(tproj_mesh%mesh_size,pawtab%basis_size))
2566      pawtab%tproj(:,:)=tproj(:,:)
2567    end if
2568    pawtab%has_tproj=2
2569  end if
2570 
2571 !==========================================================
2572 !Free temporary allocated space
2573 
2574  if (pawtab%has_tvale==1)  then
2575    LIBPAW_DEALLOCATE(nhat)
2576  end if
2577  if (reduced_vloc) then
2578    call pawrad_free(rvloc_mesh)
2579    LIBPAW_DEALLOCATE(rvlocr)
2580  end if
2581  if (reduced_ncor)  then
2582    call pawrad_free(rcore_mesh)
2583    LIBPAW_DEALLOCATE(rtncor)
2584  end if
2585  if (reduced_nval)  then
2586    call pawrad_free(rvale_mesh)
2587    LIBPAW_DEALLOCATE(rtnval)
2588  end if
2589  if (pawxcdev==0)  then
2590    LIBPAW_DEALLOCATE(pawang_tmp%angwgth)
2591    LIBPAW_DEALLOCATE(pawang_tmp%anginit)
2592    LIBPAW_DEALLOCATE(pawang_tmp%ylmr)
2593    LIBPAW_DEALLOCATE(pawang_tmp%ylmrgr)
2594  end if
2595 
2596 end subroutine pawpsp_calc

m_pawpsp/pawpsp_calc_d5 [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_calc_d5

FUNCTION

  Compute the first to the 5th derivatives of
  a given function in a pawrad mesh

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

2624 subroutine pawpsp_calc_d5(mesh,mesh_size,tcoredens)
2625 
2626 
2627 !This section has been created automatically by the script Abilint (TD).
2628 !Do not modify the following lines by hand.
2629 #undef ABI_FUNC
2630 #define ABI_FUNC 'pawpsp_calc_d5'
2631 !End of the abilint section
2632 
2633  implicit none
2634 
2635 !Arguments ------------------------------------
2636  integer,intent(in) :: mesh_size
2637  type(pawrad_type),intent(in) :: mesh
2638  real(dp),intent(inout) :: tcoredens(mesh_size,6)
2639 
2640 !Local variables-------------------------------
2641  integer,parameter :: it=1 !number of steps for smoothing function
2642  logical,parameter :: use_smooth=.true.
2643 
2644 ! *************************************************************************
2645 
2646 !calculate first derivative from density,
2647 !and store it
2648  call nderiv_gen(tcoredens(:,2),tcoredens(:,1),mesh)
2649 
2650 !get second derivative from density, and store it
2651  call paw_spline(mesh%rad,tcoredens(:,1),mesh_size,&
2652 &                zero,zero,tcoredens(:,3))
2653 
2654 !smooth functions, to avoid numerical instabilities
2655  if(use_smooth) then
2656    call paw_smooth(tcoredens(:,2),mesh_size,it)
2657    call paw_smooth(tcoredens(:,3),mesh_size,it)
2658  end if
2659 
2660 !get third derivative from first derivative:
2661  call paw_spline(mesh%rad,tcoredens(:,2),mesh_size,&
2662 &                zero,zero,tcoredens(:,4))
2663 
2664 !get fourth derivative from second derivative:
2665  call paw_spline(mesh%rad,tcoredens(:,3),mesh_size,&
2666 &                zero,zero,tcoredens(:,5))
2667 
2668 !smooth 3rd and 4th order derivatives
2669  if(use_smooth) then
2670    call paw_smooth(tcoredens(:,4),mesh_size,it)
2671    call paw_smooth(tcoredens(:,5),mesh_size,it)
2672  end if
2673 
2674 !get fifth derivative from third derivative:
2675  call paw_spline(mesh%rad,tcoredens(:,4),mesh_size,&
2676 &                zero,zero,tcoredens(:,6))
2677 
2678 !smooth 5th order derivative
2679  if(use_smooth) then
2680    call paw_smooth(tcoredens(:,6),mesh_size,it)
2681  end if
2682 
2683 end subroutine pawpsp_calc_d5

m_pawpsp/pawpsp_cg [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_cg

FUNCTION

 Compute sine transform to transform from n(r) to n(q).
 Computes integrals on (generalized) grid using corrected trapezoidal integration.

INPUTS

  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=q grid values (bohr**-1).
  radmesh <type(pawrad_type)>=data containing radial grid information
  nr(:)=n(r) on radial grid.

OUTPUT

  dnqdq0= 1/q dn(q)/dq for q=0
  d2nqdq0 = Gives contribution of d2(tNcore(q))/d2q for q=0
            compute \int{(16/15)*pi^5*n(r)*r^6* dr}
  nq(mqgrid)= n(q)
            = 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 n(r))dr].
  yp1,ypn=derivatives of n(q) wrt q at q=0 and q=qmax (needed for spline fitter).

PARENTS

      dfpt_eltfrxc,m_pawpsp,m_psps

CHILDREN

SOURCE

492 subroutine pawpsp_cg(dnqdq0,d2nqdq0,mqgrid,qgrid,nq,radmesh,nr,yp1,ypn)
493 
494 
495 !This section has been created automatically by the script Abilint (TD).
496 !Do not modify the following lines by hand.
497 #undef ABI_FUNC
498 #define ABI_FUNC 'pawpsp_cg'
499 !End of the abilint section
500 
501  implicit none
502 
503 !Arguments----------------------------------------------------------
504 !scalars
505  integer,intent(in) :: mqgrid
506  real(dp),intent(out) :: dnqdq0,d2nqdq0,yp1,ypn
507  type(pawrad_type),intent(in) :: radmesh
508 !arrays
509  real(dp),intent(in) :: nr(:)
510  real(dp),intent(in) :: qgrid(mqgrid)
511  real(dp),intent(out) :: nq(mqgrid)
512 
513 !Local variables-------------------------------
514 !scalars
515  integer :: iq,ir,mesh_size
516  real(dp) :: aexp,arg,bexp,dn,r0tor1,r1torm,rm,rmtoin
517  logical :: begin_r0
518  !character(len=500) :: msg
519 !arrays
520  real(dp),allocatable :: ff(:),rnr(:)
521 
522 ! *************************************************************************
523 
524  mesh_size=min(size(nr),radmesh%mesh_size)
525  LIBPAW_ALLOCATE(ff,(mesh_size))
526  LIBPAW_ALLOCATE(rnr,(mesh_size))
527  ff=zero;rnr=zero
528 
529  do ir=1,mesh_size
530    rnr(ir)=radmesh%rad(ir)*nr(ir)
531  end do
532 
533 !Is mesh beginning with r=0 ?
534  begin_r0=(radmesh%rad(1)<1.d-20)
535 
536 !Adjustment of an exponentional at r_max (n_exp(r)=aexp*Exp[-bexp*r])
537  rm=radmesh%rad(mesh_size)
538  dn=one/(12._dp*radmesh%stepint*radmesh%radfact(mesh_size)) &
539 & *( 3._dp*nr(mesh_size-4) &
540 &  -16._dp*nr(mesh_size-3) &
541 &  +36._dp*nr(mesh_size-2) &
542 &  -48._dp*nr(mesh_size-1) &
543 &  +25._dp*nr(mesh_size))
544  if (dn<0._dp.and. &
545 & abs(radmesh%rad(mesh_size)*nr(mesh_size))>1.d-20) then
546    bexp=-dn/nr(mesh_size)
547    if (bexp * rm > 50._dp) then
548      ! This solves the problem with the weird core charge used in v4[62] in which bexp x rm ~= 10^3
549      !write(msg,"(a,es16.8)")"Tooooo large bexp * rm: ", bexp*rm, ", setting aexp to 0"
550      !MSG_WARNING(msg)
551      bexp=0.001_dp;aexp=zero
552    else
553      aexp=nr(mesh_size)*exp(bexp*rm)
554      if (abs(aexp)<1.d-20) then
555        bexp=0.001_dp;aexp=zero
556      end if
557    end if
558  else
559    bexp=0.001_dp;aexp=zero
560  end if
561 
562 !===========================================
563 !=== Compute n(q) for q=0 separately
564 !===========================================
565 
566 !Integral from 0 to r1 (only if r1<>0)
567  r0tor1=zero
568  if (.not.begin_r0) r0tor1=(rnr(1)*radmesh%rad(1)**2)/3.d0
569 
570 !Integral from r1 to rmax
571  do ir=1,mesh_size
572    if (abs(rnr(ir))>1.d-20) ff(ir)=rnr(ir)*radmesh%rad(ir)
573  end do
574  call simp_gen(r1torm,ff,radmesh)
575 
576 !Integral from rmax to infinity
577 !This part is approximated using an exponential density aexp*Exp[-bexp*r]
578 !(formulae obtained with mathematica)
579  rmtoin=aexp*exp(-bexp*rm)/bexp**3*(two+two*bexp*rm+bexp*bexp*rm*rm)
580 
581 !Some of the three parts
582  nq(1)=four_pi*(r0tor1+r1torm+rmtoin)
583 
584 !===========================================
585 !=== Compute n(q) for other q''s
586 !===========================================
587 
588 !Loop over q values
589  do iq=2,mqgrid
590    arg=two_pi*qgrid(iq)
591 
592 !  Integral from 0 to r1 (only if r1<>0)
593    r0tor1=zero;if (.not.begin_r0) &
594 &   r0tor1=nr(1)*(sin(arg*radmesh%rad(1))/arg/arg&
595 &   -radmesh%rad(1)*cos(arg*radmesh%rad(1))/arg)
596 
597 !  Integral from r1 to rmax
598    do ir=1,mesh_size
599      if (abs(rnr(ir))>1.d-20) ff(ir)=sin(arg*radmesh%rad(ir))*rnr(ir)
600    end do
601    call simp_gen(r1torm,ff,radmesh)
602 
603 !  Integral from rmax to infinity
604 !  This part is approximated using an exponential density aexp*Exp[-bexp*r]
605 !  (formulae obtained with mathematica)
606    rmtoin=aexp*exp(-bexp*rm)/(arg**2+bexp**2)**2 &
607 &   *(arg*(two*bexp+arg**2*rm+bexp**2*rm)*cos(arg*rm) &
608 &   +(arg**2*(bexp*rm-one)+bexp**2*(bexp*rm+one))*sin(arg*rm))
609 
610 !  Store q^2 v(q)
611    nq(iq)=two/qgrid(iq)*(r0tor1+r1torm+rmtoin)
612  end do
613 
614 !===========================================
615 !=== Compute derivatives of n(q)
616 !=== at ends of interval
617 !===========================================
618 
619 !yp(0)=zero
620  yp1=zero
621 
622 !yp(qmax)=$ 2\int_0^\infty[(-\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r) r n(r) dr]$
623  arg=two_pi*qgrid(mqgrid)
624 
625 !Integral from 0 to r1 (only if r1<>0)
626  r0tor1=zero;if (.not.begin_r0) &
627 & r0tor1=two_pi*nr(1)*(3.d0*radmesh%rad(1)/arg /arg*cos(arg*radmesh%rad(1))+ &
628 & (radmesh%rad(1)**2/arg-3.0d0/arg**3)*sin(arg*radmesh%rad(1)))
629 
630 !Integral from r1 to rmax
631  do ir=1,mesh_size
632    if (abs(rnr(ir))>1.d-20) ff(ir)=(two_pi*radmesh%rad(ir)*cos(arg*radmesh%rad(ir)) &
633 &   - sin(arg*radmesh%rad(ir))/qgrid(mqgrid)) *rnr(ir)
634  end do
635  call simp_gen(r1torm,ff,radmesh)
636 
637 !Integral from rmax to infinity
638 !This part is approximated using an exponential density aexp*Exp[-bexp*r]
639 !(formulae obtained with mathematica)
640  rmtoin=-one/(qgrid(mqgrid)*(arg**2+bexp**2)**3) &
641 & *aexp*exp(-bexp*rm) &
642 & *((arg**5*rm-two_pi*arg**4*qgrid(mqgrid)*rm*(bexp*rm-two) &
643 & +two*arg**3*bexp*(bexp*rm+one)+arg*bexp**3*(bexp*rm+two) &
644 & -four_pi*arg**2*bexp*qgrid(mqgrid)*(bexp**2*rm**2-three) &
645 & -two_pi*bexp**3*qgrid(mqgrid)*(bexp**2*rm**2+two*bexp*rm+two))*cos(arg*rm) &
646 & +(two*arg**2*bexp**3*rm+two_pi*arg**5*qgrid(mqgrid)*rm**2 &
647 & +arg**4*(bexp*rm-one)+bexp**4*(bexp*rm+one) &
648 & +four_pi*arg**3*qgrid(mqgrid)*(bexp**2*rm**2+two*bexp*rm-one) &
649 & +two_pi*arg*bexp**2*qgrid(mqgrid)*(bexp**2*rm**2+four*bexp*rm+6._dp))*sin(arg*rm))
650 
651 !Some of the three parts
652  ypn=two/qgrid(mqgrid)*(r0tor1+r1torm+rmtoin)
653 
654 !===========================================
655 !=== Compute 1/q dn(q)/dq at q=0
656 !===========================================
657 
658 !Integral from 0 to r1 (only if r1<>0)
659  r0tor1=zero
660  if (.not.begin_r0) r0tor1=(rnr(1)*radmesh%rad(1)**4)/5.d0
661 
662 !Integral from r1 to rmax
663  do ir=1,mesh_size
664    if (abs(rnr(ir))>1.d-20) ff(ir)=rnr(ir)*radmesh%rad(ir)**3
665  end do
666  call simp_gen(r1torm,ff,radmesh)
667 
668 !Integral from rmax to infinity
669 !This part is approximated using an exponential density aexp*Exp[-bexp*r]
670 !(formulae obtained with mathematica)
671  rmtoin=aexp*exp(-bexp*rm)/bexp**5 &
672 & *(24._dp+24._dp*bexp*rm+12._dp*bexp**2*rm**2+four*bexp**3*rm**3+bexp**4*rm**4)
673 
674 !Some of the three parts
675  dnqdq0=-(2.d0/3.d0)*two_pi**3*(r0tor1+r1torm+rmtoin)
676 
677  LIBPAW_DEALLOCATE(ff)
678  LIBPAW_DEALLOCATE(rnr)
679 
680  d2nqdq0 = 1_dp
681 
682 end subroutine pawpsp_cg

m_pawpsp/pawpsp_header_type [ Types ]

[ Top ] [ m_pawpsp ] [ Types ]

NAME

 pawpsp_header_type

FUNCTION

 For PAW, header related data

SOURCE

 89  type, public :: pawpsp_header_type
 90 
 91 !Integer scalars
 92   integer :: basis_size    ! Number of elements of the wf basis ((l,n) quantum numbers)
 93   integer :: l_size        ! Maximum value of l+1 leading to a non zero Gaunt coefficient
 94   integer :: lmn_size      ! Number of elements of the paw basis
 95   integer :: mesh_size     ! Dimension of (main) radial mesh
 96   integer :: pawver        ! Version number of paw psp format
 97   integer :: shape_type    ! Type of shape function
 98   real(dp) :: rpaw         ! Radius for paw spheres
 99   real(dp) :: rshp         ! Cut-off radius of shape function
100 
101  end type pawpsp_header_type

m_pawpsp/pawpsp_lo [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_lo

FUNCTION

 Compute sine transform to transform from V(r) to q^2 V(q).
 Computes integrals on (generalized) grid using corrected trapezoidal integration.

INPUTS

  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=q grid values (bohr**-1).
  radmesh <type(pawrad_type)>=data containing radial grid information
  vloc(:)=V(r) on radial grid.
  zion=nominal valence charge of atom.

OUTPUT

  epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
  q2vq(mqgrid)
   =q^2 V(q)
   = -\frac{Zv}{\pi}
     + q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
  yp1,ypn=derivatives of q^2 V(q) wrt q at q=0 and q=qmax (needed for spline fitter).

PARENTS

      m_pawpsp

CHILDREN

SOURCE

315 subroutine pawpsp_lo(epsatm,mqgrid,qgrid,q2vq,radmesh,vloc,yp1,ypn,zion)
316 
317 
318 !This section has been created automatically by the script Abilint (TD).
319 !Do not modify the following lines by hand.
320 #undef ABI_FUNC
321 #define ABI_FUNC 'pawpsp_lo'
322 !End of the abilint section
323 
324  implicit none
325 
326 !Arguments----------------------------------------------------------
327 !scalars
328  integer,intent(in) :: mqgrid
329  real(dp),intent(in) :: zion
330  real(dp),intent(out) :: epsatm,yp1,ypn
331  type(pawrad_type),intent(in) :: radmesh
332 !arrays
333  real(dp),intent(in) :: qgrid(mqgrid)
334  real(dp),intent(in) :: vloc(:)
335  real(dp),intent(out) :: q2vq(mqgrid)
336 
337 !Local variables ------------------------------
338 !scalars
339  integer :: iq,ir,irmax,mesh_size
340  real(dp) :: arg,r0tor1,r1torm,rmtoin
341  logical :: begin_r0
342 !arrays
343  real(dp),allocatable :: ff(:),rvpz(:)
344 
345 !************************************************************************
346 
347  mesh_size=size(vloc)
348  irmax=pawrad_ifromr(radmesh,min(20._dp,radmesh%rmax))
349  irmax=min(irmax,mesh_size)
350 
351 !Particular case of a zero potential
352  if (maxval(abs(vloc(1:irmax)))<=1.e-20_dp) then
353    q2vq=zero;yp1=zero;ypn=zero;epsatm=zero
354    return
355  end if
356 
357  LIBPAW_ALLOCATE(ff,(mesh_size))
358  LIBPAW_ALLOCATE(rvpz,(mesh_size))
359  ff=zero;rvpz=zero
360 
361 !Is mesh beginning with r=0 ?
362  begin_r0=(radmesh%rad(1)<1.e-20_dp)
363 
364 !Store r.V+Z
365  do ir=1,irmax
366    rvpz(ir)=radmesh%rad(ir)*vloc(ir)+zion
367  end do
368 
369 !===========================================
370 !=== Compute q^2 v(q) for q=0 separately
371 !===========================================
372 
373 !Integral from 0 to r1 (only if r1<>0)
374  r0tor1=zero;if (.not.begin_r0) &
375 & r0tor1=(zion*0.5_dp+radmesh%rad(1)*vloc(1)/3._dp)*radmesh%rad(1)**2
376 
377 !Integral from r1 to rmax
378  do ir=1,irmax
379    if (abs(rvpz(ir))>1.e-20_dp) then
380      ff(ir)=radmesh%rad(ir)*rvpz(ir)
381    end if
382  end do
383 
384  call simp_gen(r1torm,ff,radmesh)
385 
386 !Integral from rmax to infinity
387 !This part is neglected... might be improved.
388  rmtoin=zero
389 
390 !Some of the three parts
391  epsatm=four_pi*(r0tor1+r1torm+rmtoin)
392 
393  q2vq(1)=-zion/pi
394 
395 !===========================================
396 !=== Compute q^2 v(q) for other q''s
397 !===========================================
398 
399 !Loop over q values
400  do iq=2,mqgrid
401    arg=two_pi*qgrid(iq)
402 
403 !  Integral from 0 to r1 (only if r1<>0)
404    r0tor1=zero;if (.not.begin_r0) &
405 &   r0tor1=( vloc(1)/arg*sin(arg*radmesh%rad(1)) &
406 &   -rvpz(1)    *cos(arg*radmesh%rad(1)) +zion )/pi
407 
408 !  Integral from r1 to rmax
409    do ir=1,irmax
410      if (abs(rvpz(ir))>1.e-20_dp) ff(ir)=sin(arg*radmesh%rad(ir))*rvpz(ir)
411    end do
412    call simp_gen(r1torm,ff,radmesh)
413 
414 !  Integral from rmax to infinity
415 !  This part is neglected... might be improved.
416    rmtoin=zero
417 
418 !  Store q^2 v(q)
419    q2vq(iq)=-zion/pi + two*qgrid(iq)*(r0tor1+r1torm+rmtoin)
420  end do
421 
422 !===========================================
423 !=== Compute derivatives of q^2 v(q)
424 !=== at ends of interval
425 !===========================================
426 
427 !yp(0)=zero
428  yp1=zero
429 
430 !yp(qmax)=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$
431  arg=two_pi*qgrid(mqgrid)
432 
433 !Integral from 0 to r1 (only if r1<>0)
434  r0tor1=zero;if (.not.begin_r0) &
435 & r0tor1=zion*radmesh%rad(1)                  *sin(arg*radmesh%rad(1)) &
436 & +three*radmesh%rad(1)*vloc(1)/arg         *cos(arg*radmesh%rad(1)) &
437 & +(radmesh%rad(1)**2-one/arg**2)*vloc(1)*sin(arg*radmesh%rad(1))
438 
439 !Integral from r1 to rmax
440  do ir=1,irmax
441    if (abs(rvpz(ir))>1.e-20_dp) ff(ir)=( arg*radmesh%rad(ir)*cos(arg*radmesh%rad(ir)) &
442 &   +                    sin(arg*radmesh%rad(ir))) *rvpz(ir)
443  end do
444  call simp_gen(r1torm,ff,radmesh)
445 
446 !Integral from rmax to infinity
447 !This part is neglected... might be improved.
448  rmtoin=zero
449 
450 !Some of the three parts
451  ypn=two*(r0tor1+r1torm+rmtoin)
452 
453  LIBPAW_DEALLOCATE(ff)
454  LIBPAW_DEALLOCATE(rvpz)
455 
456 end subroutine pawpsp_lo

m_pawpsp/pawpsp_main [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_main

FUNCTION

 Reads a PAW dataset (atomic data)

INPUTS

  filpsp=name of the file containing the PAW dataset
  usewvl=1 if we use a wavelet basis, 0 other wise (plane waves)
  icoulomb=1 if we use a Poisson routine with wavelets, 0 otherwise
  ixc=index of the XC correlation functional
  xclevel=type of XC functional (1=LDA, 2=GGA, ...)
  pawxcdev=order of the developement of the PAW on-site terms
           (0: full calculation, 1: order 1, 2:order 2)
  usexcnhat=flag controlling the use of compensation charge (nhat) in XC potential
  qgrid_ff=size of the mesh for the sin FFT transform of the non-local projectors (form factors)
           (plane waves only, 0 otherwise)
  qgrid_vl=size of the mesh for the sin FFT transform of the local potential
           (plane waves only, 0 otherwise)
  ffspl=sin FFT transform of the non-local projectors (form factors) (plane waves only)
  vlspl=sin FFT transform of the local potential (plane waves only)
  epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
  xcccrc=XC core correction cutoff radius (bohr)
  zionpsp=valence of atom as specified in input file
  znuclpsp=atomic number of atom as specified in input file
  ===== Optional arguments for wvl =====
    wvl_ngauss
  ===== Other optional arguments =====
    psxml=datastructure containing a XMP PAW dataset
    comm_mpi=MPI communicator
    xc_denpos=tolerance on density/potential for the calculation of XC potential
              (if density<xc_denpos, density=zero)

OUTPUT

  pawrad <type(pawrad_type)>=data containing PAW radial grid information
  pawtab <type(pawtab_type)>=data containing the PAW dataset (partial waves...)

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

4785 subroutine pawpsp_main( &
4786 & pawrad,pawtab,&
4787 & filpsp,usewvl,icoulomb,ixc,xclevel,pawxcdev,usexcnhat,&
4788 & qgrid_ff,qgrid_vl,ffspl,vlspl,epsatm,xcccrc,zionpsp,znuclpsp,&
4789 & wvl_ngauss,psxml,comm_mpi,xc_denpos)
4790 
4791 
4792 !This section has been created automatically by the script Abilint (TD).
4793 !Do not modify the following lines by hand.
4794 #undef ABI_FUNC
4795 #define ABI_FUNC 'pawpsp_main'
4796 !End of the abilint section
4797 
4798  implicit none
4799 
4800 !Arguments ------------------------------------
4801 !scalars
4802  integer,intent(in) :: icoulomb,ixc
4803  integer,intent(in) :: pawxcdev,usewvl,usexcnhat,xclevel
4804  integer,optional,intent(in) :: comm_mpi
4805  real(dp),intent(in):: zionpsp,znuclpsp
4806  real(dp),optional,intent(in) :: xc_denpos
4807  real(dp),intent(out) :: epsatm,xcccrc
4808  character(len=fnlen),intent(in):: filpsp   ! name of the psp file
4809  type(pawrad_type),intent(inout) :: pawrad
4810  type(pawtab_type),intent(inout) :: pawtab
4811  type(paw_setup_t),optional,intent(in) :: psxml
4812 !arrays
4813  integer,optional,intent(in) :: wvl_ngauss(2)
4814  real(dp),intent(in) :: qgrid_ff(:),qgrid_vl(:)
4815  real(dp),intent(inout) :: ffspl(:,:,:)
4816  real(dp),intent(out) :: vlspl(:,:)
4817 
4818 !Local variables-------------------------------
4819  integer :: has_tproj,has_wvl,ipsp,lmax,lloc,lnmax,mmax,me,mqgrid_ff,mqgrid_vl
4820  integer :: pspcod,pspxc,usexml
4821  real(dp),parameter :: xc_denpos_default=tol14
4822  real(dp) :: my_xc_denpos,r2well,zion,znucl
4823  character(len=500) :: msg
4824  type(pawpsp_header_type) :: pawpsp_header
4825 !arrays
4826 
4827 ! *************************************************************************
4828 
4829 !Check consistency of parameters
4830  if (icoulomb/= 0.or.usewvl==1) then
4831    if (.not.present(wvl_ngauss)) then
4832      msg='usewvl==1 or icoulomb/=0: a mandatory argument is missing!'
4833      MSG_BUG(msg)
4834    end if
4835  end if
4836 
4837  mqgrid_ff=size(qgrid_ff)
4838  mqgrid_vl=size(qgrid_vl)
4839  lnmax=size(ffspl,3)
4840  if (size(ffspl,1)/=mqgrid_ff.or.size(ffspl,2)/=2) then
4841    msg='invalid sizes for ffspl!'
4842    MSG_BUG(msg)
4843  end if
4844  if (size(vlspl,1)/=mqgrid_vl.or.size(vlspl,2)/=2) then
4845    msg='invalid sizes for vlspl!'
4846    MSG_BUG(msg)
4847  end if
4848 
4849  my_xc_denpos=xc_denpos_default;if (present(xc_denpos)) my_xc_denpos=xc_denpos
4850  pawtab%usexcnhat=usexcnhat
4851  me=0;if (present(comm_mpi))me=xmpi_comm_rank(comm_mpi)
4852 
4853  has_wvl=0; if (usewvl==1.or.icoulomb/=0) has_wvl=1
4854  has_tproj=0; if (usewvl==1) has_tproj=1
4855  call pawtab_set_flags(pawtab,has_tvale=1,has_wvl=has_wvl,has_tproj=has_tproj)
4856 
4857  if(me==0) then
4858    write(msg, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(filpsp)
4859    call wrtout(ab_out,  msg,'COLL')
4860    call wrtout(std_out,  msg,'COLL')
4861 
4862 !  This checks if file is xml or UPF
4863 !  It sets usexml as well
4864    call pawpsp_check_xml_upf(filpsp)
4865 
4866 !  ----------------------------------------------------------------------------
4867    if (usexml /= 1) then
4868 !    Open the atomic data file, and read the three first lines
4869      open (unit=tmp_unit,file=filpsp,form='formatted',status='old')
4870      rewind (unit=tmp_unit)
4871 !    Read first 3 lines of psp file:
4872      call pawpsp_read_header(tmp_unit,lloc,lmax,mmax,pspcod,&
4873 &     pspxc,r2well,zion,znucl)
4874 
4875    else if (usexml == 1 .and. present(psxml)) then
4876      write(msg,'(a,a)')  &
4877 &     '- pawpsp : Reading pseudopotential header in XML form from ', trim(filpsp)
4878      call wrtout(ab_out,msg,'COLL')
4879      call wrtout(std_out,  msg,'COLL')
4880 
4881 !    Return header information
4882      call pawpsp_read_header_xml(lloc,lmax,pspcod,&
4883 &     pspxc,psxml,r2well,zion,znucl)
4884 !    Fill in pawpsp_header object:
4885      call pawpsp_read_pawheader(pawpsp_header%basis_size,&
4886 &   lmax,pawpsp_header%lmn_size,&
4887 &   pawpsp_header%l_size,pawpsp_header%mesh_size,&
4888 &   pawpsp_header%pawver,psxml,&
4889 &   pawpsp_header%rpaw,pawpsp_header%rshp,pawpsp_header%shape_type)
4890    end if
4891 
4892 !  Check data for consistency against main routine input
4893    call pawpsp_consistency()
4894 
4895 !  Read rest of the PSP file
4896    if (pspcod==7) then
4897 !    ABINIT proprietary format
4898      call pawpsp_7in(epsatm,ffspl,icoulomb,ixc,&
4899 &     lmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,&
4900 &     pawrad,pawtab,pawxcdev,qgrid_ff,qgrid_vl,&
4901 &     usewvl,usexcnhat,vlspl,xcccrc,xclevel,my_xc_denpos,zion,znucl)
4902 
4903    else if (pspcod==17)then
4904 !    XML format
4905      ipsp=1
4906      call pawpsp_17in(epsatm,ffspl,icoulomb,ipsp,ixc,lmax,&
4907 &     lnmax,mmax,mqgrid_ff,mqgrid_vl,pawpsp_header,pawrad,pawtab,&
4908 &     pawxcdev,qgrid_ff,qgrid_vl,usewvl,usexcnhat,vlspl,xcccrc,&
4909 &     xclevel,my_xc_denpos,zion,znucl)
4910 
4911    end if
4912  end if!me==0
4913 
4914  close(unit=tmp_unit)
4915 
4916  write(msg,'(3a)') ' pawpsp: atomic psp has been read ',&
4917 & ' and splines computed',ch10
4918  call wrtout(ab_out,msg,'COLL')
4919  call wrtout(std_out,  msg,'COLL')
4920 
4921 !Communicate PAW objects
4922  if(present(comm_mpi)) then
4923    if(xmpi_comm_size(comm_mpi)>1) then
4924      call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc)
4925    end if
4926  end if
4927 
4928 !WVL+PAW:
4929  if(icoulomb/=0.or.usewvl==1) then
4930    if(present(comm_mpi))then
4931     call pawpsp_wvl(filpsp,pawrad,pawtab,usewvl,wvl_ngauss,comm_mpi)
4932    else
4933     call pawpsp_wvl(filpsp,pawrad,pawtab,usewvl,wvl_ngauss)
4934    end if
4935  end if
4936 
4937 contains

m_pawpsp/pawpsp_nl [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_nl

FUNCTION

 Make paw projector form factors f_l(q) for each l

INPUTS

  indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn
  lmnmax=max number of (l,m,n) components
  lnmax=max number of (l,n) components
  mqgrid=number of grid points for q grid
  qgrid(mqgrid)=values at which form factors are returned
  radmesh <type(pawrad_type)>=data containing radial grid information
  wfll(:,lnmax)=paw projector on radial grid

OUTPUT

  ffspl(mqgrid,2,lnmax)= form factor f_l(q) and second derivative

NOTES

  u_l(r) is the paw projector (input as wfll);
  j_l(q) is a spherical Bessel function;
  f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r)  r dr]$

PARENTS

      m_pawpsp,pawinit

CHILDREN

SOURCE

141 subroutine pawpsp_nl(ffspl,indlmn,lmnmax,lnmax,mqgrid,qgrid,radmesh,wfll)
142 
143 
144 !This section has been created automatically by the script Abilint (TD).
145 !Do not modify the following lines by hand.
146 #undef ABI_FUNC
147 #define ABI_FUNC 'pawpsp_nl'
148 !End of the abilint section
149 
150  implicit none
151 
152 !Arguments ------------------------------------
153 !scalars
154  integer,intent(in) :: lmnmax,lnmax,mqgrid
155  type(pawrad_type),intent(in) :: radmesh
156 !arrays
157  integer,intent(in) :: indlmn(6,lmnmax)
158  real(dp),intent(in) :: qgrid(mqgrid)
159  real(dp),intent(in) ::  wfll(:,:)
160  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax)
161 
162 !Local variables-------------------------------
163 !scalars
164  integer :: ilmn,iln,iln0,iq,ir,ll,meshsz,mmax
165  real(dp),parameter :: eps=tol14**4,TOLJ=0.001_dp
166  real(dp) :: arg,argn,bes
167  real(dp) :: besp,qr
168  real(dp) :: yp1,ypn
169  character(len=100) :: msg
170  type(pawrad_type) :: tmpmesh
171 !arrays
172  real(dp),allocatable :: ff(:),gg(:),rr(:),rr2(:),rr2wf(:),rrwf(:),work(:)
173 
174 !*************************************************************************
175 
176 !Is mesh beginning with r=0 ?
177  if (radmesh%rad(1)>tol10) then
178    msg='Radial mesh cannot begin with r<>0!'
179    MSG_BUG(msg)
180  end if
181 
182  meshsz=size(wfll,1)
183  if (meshsz>radmesh%mesh_size) then
184    msg='wrong size for wfll!'
185    MSG_BUG(msg)
186  end if
187 
188 !Init. temporary arrays and variables
189  LIBPAW_ALLOCATE(ff,(meshsz))
190  LIBPAW_ALLOCATE(gg,(meshsz))
191  LIBPAW_ALLOCATE(rr,(meshsz))
192  LIBPAW_ALLOCATE(rr2,(meshsz))
193  LIBPAW_ALLOCATE(rrwf,(meshsz))
194  LIBPAW_ALLOCATE(rr2wf,(meshsz))
195  LIBPAW_ALLOCATE(work,(mqgrid))
196  rr(1:meshsz) =radmesh%rad(1:meshsz)
197  rr2(1:meshsz)=two_pi*rr(1:meshsz)*rr(1:meshsz)
198  argn=two_pi*qgrid(mqgrid)
199  mmax=meshsz
200 
201 !Loop on (l,n) projectors
202  iln0=0
203  do ilmn=1,lmnmax
204    iln=indlmn(5,ilmn)
205    if(iln>iln0) then
206      iln0=iln;ll=indlmn(1,ilmn)
207 
208      ir=meshsz
209      do while (abs(wfll(ir,iln))<eps)
210        ir=ir-1
211      end do
212      mmax=min(ir+1,meshsz)
213      if (mmax/=radmesh%int_meshsz) then
214        call pawrad_init(tmpmesh,mesh_size=meshsz,mesh_type=radmesh%mesh_type, &
215 &       rstep=radmesh%rstep,lstep=radmesh%lstep,r_for_intg=rr(mmax))
216      else
217        call pawrad_copy(radmesh,tmpmesh)
218      end if
219 
220      rrwf(:) =rr (:)*wfll(:,iln)
221      rr2wf(:)=rr2(:)*wfll(:,iln)
222 
223 !    1-Compute f_l(0<q<qmax)
224      if (mqgrid>2) then
225        do iq=2,mqgrid-1
226          arg=two_pi*qgrid(iq)
227          do ir=1,mmax
228            qr=arg*rr(ir)
229            call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ)
230            ff(ir)=bes*rrwf(ir)
231          end do
232          call simp_gen(ffspl(iq,1,iln),ff,tmpmesh)
233        end do
234      end if
235 
236 !    2-Compute f_l(q=0) and first derivative
237      ffspl(1,1,iln)=zero;yp1=zero
238      if (ll==0) then
239        call simp_gen(ffspl(1,1,iln),rrwf,tmpmesh)
240      end if
241      if (ll==1) then
242        call simp_gen(yp1,rr2wf,tmpmesh)
243        yp1=yp1*third
244      end if
245 
246 !    3-Compute f_l(q=qmax) and first derivative
247      if (mqgrid>1) then
248 !      if (ll==0.or.ll==1) then
249        do ir=1,mmax
250          qr=argn*rr(ir)
251          call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ)
252          ff(ir)=bes*rrwf(ir)
253          gg(ir)=besp*rr2wf(ir)
254        end do
255        call simp_gen(ffspl(mqgrid,1,iln),ff,tmpmesh)
256        call simp_gen(ypn,gg,tmpmesh)
257      else
258        ypn=yp1
259      end if
260 
261 !    4-Compute second derivative of f_l(q)
262      call paw_spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
263 
264      call pawrad_free(tmpmesh)
265 
266 !    End loop on (l,n) projectors
267    end if
268  end do
269 
270  LIBPAW_DEALLOCATE(ff)
271  LIBPAW_DEALLOCATE(gg)
272  LIBPAW_DEALLOCATE(rr)
273  LIBPAW_DEALLOCATE(rr2)
274  LIBPAW_DEALLOCATE(rrwf)
275  LIBPAW_DEALLOCATE(rr2wf)
276  LIBPAW_DEALLOCATE(work)
277 
278 end subroutine pawpsp_nl

m_pawpsp/pawpsp_read [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_read

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

   File format of formatted PAW psp input (the 3 first lines
   have already been read in calling -pspatm- routine) :
   (1) title (character) line
   (2) psps%znuclpsp(ipsp), 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 core_mesh_size)
   (23) comment(character)
   (24) radial mesh index for pseudo_core_density
   (25) tcore_density (for ir=1 to core_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), 2=VH(tnzc) without nhat in XC)
   (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 shape_mesh_size)
   (36) comment(character)
   (37) radial mesh index for pseudo_valence_density
   (38) tvale(r) (for ir=1 to vale_mesh_size)

   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 ; gl(r)=numeric (read from psp file)
   shape_type= 1 ; gl(r)=k(r).r^l; k(r)=exp[-(r/sigma)**lambda]
   shape_type= 2 ; gl(r)=k(r).r^l; k(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2 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

PARENTS

      m_pawpsp

CHILDREN

SOURCE

 772 subroutine pawpsp_read(core_mesh,funit,imainmesh,lmax,&
 773 & ncore,nmesh,pawrad,pawtab,pspversion,radmesh,save_core_msz,&
 774 & tncore,tnvale,tproj,tproj_mesh,usexcnhat_in,usexcnhat_out,vale_mesh,&
 775 & vlocopt,vlocr,vloc_mesh,znucl)
 776 
 777 
 778 !This section has been created automatically by the script Abilint (TD).
 779 !Do not modify the following lines by hand.
 780 #undef ABI_FUNC
 781 #define ABI_FUNC 'pawpsp_read'
 782 !End of the abilint section
 783 
 784  implicit none
 785 
 786 !Arguments ------------------------------------
 787  integer,intent(in):: funit,lmax,usexcnhat_in
 788  integer,intent(out) :: imainmesh,pspversion,usexcnhat_out,vlocopt
 789  logical,intent(in) :: save_core_msz
 790  real(dp),intent(in):: znucl
 791 !arrays
 792  real(dp),pointer :: ncore(:),tncore(:),tnvale(:),tproj(:,:),vlocr(:)
 793  type(pawrad_type),intent(inout) :: pawrad
 794  type(pawrad_type),intent(out)::core_mesh,tproj_mesh,vale_mesh,vloc_mesh
 795  type(pawrad_type),pointer :: radmesh(:)
 796  type(pawtab_type),intent(inout) :: pawtab
 797  integer,intent(out)::nmesh
 798 
 799 !Local variables-------------------------------
 800  integer :: creatorid,imsh
 801  integer :: icoremesh,ishpfmesh,ivalemesh,ivlocmesh
 802  integer :: ib,il,ilm,ilmn,iln,iprojmesh
 803  integer :: ii,ir,iread1,iread2,jj
 804  integer :: msz,pngau_,ptotgau_
 805  real(dp):: rc,rread1,rread2
 806  real(dp) :: yp1,ypn
 807 !arrays
 808  integer,allocatable :: nprj(:)
 809  real(dp),allocatable :: shpf(:,:),val(:),vhnzc(:)
 810  real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:)
 811  character :: blank=' ',numb=' '
 812  character(len=80) :: pspline
 813  character(len=500) :: msg,submsg
 814  logical :: read_gauss=.false.
 815  type(pawrad_type)::shpf_mesh
 816 
 817 ! *************************************************************************
 818 
 819 !==========================================================
 820 !Read lines 4 to 11 of the header
 821 
 822 !This is important for BigDFT in standalone mode
 823  call pawpsp_read_header_2(funit,pspversion,pawtab%basis_size,pawtab%lmn_size)
 824 
 825 !Check pspversion for wvl-paw
 826  if(pspversion<4 .and. pawtab%has_wvl>0)then
 827    write(msg, '(a,i2,a,a)' )&
 828 &   'In reading atomic psp file, finds pspversion=',pspversion,ch10,&
 829 &   'For WVL-PAW, pspversion >= 4 is required.'
 830    MSG_BUG(msg)
 831  end if
 832 
 833 
 834 !Have to maintain compatibility with Abinit v4.2.x
 835  if (pspversion==1) then
 836    LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size))
 837    read(funit,*) (pawtab%orbitals(ib), ib=1,pawtab%basis_size)
 838    pawtab%l_size=2*maxval(pawtab%orbitals)+1
 839    nmesh=3
 840    LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh))
 841    read(funit,'(a80)') pspline
 842    radmesh(1)%lstep=zero
 843    read(unit=pspline,fmt=*,err=10,end=10) radmesh(1)%mesh_type,&
 844 &   radmesh(1)%rstep,radmesh(1)%lstep
 845    10 read(funit,*) pawtab%rpaw
 846    read(funit,*) radmesh(1)%mesh_size,radmesh(2)%mesh_size,&
 847 &   radmesh(3)%mesh_size
 848    read(funit,'(a80)') pspline
 849    pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99
 850    read(unit=pspline,fmt=*,err=11,end=11) pawtab%shape_type,&
 851 &   pawtab%shape_lambda,pawtab%shape_sigma
 852    11 read(funit,*) creatorid
 853    if (pawtab%shape_type==3) pawtab%shape_type=-1
 854    radmesh(2)%mesh_type=radmesh(1)%mesh_type
 855    radmesh(3)%mesh_type=radmesh(1)%mesh_type
 856    radmesh(2)%rstep=radmesh(1)%rstep
 857    radmesh(3)%rstep=radmesh(1)%rstep
 858    radmesh(2)%lstep=radmesh(1)%lstep
 859    radmesh(3)%lstep=radmesh(1)%lstep
 860  else
 861 
 862 !  Here psp file for Abinit 4.3+
 863    LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size))
 864    read(funit,*) (pawtab%orbitals(ib), ib=1,pawtab%basis_size)
 865    pawtab%l_size=2*maxval(pawtab%orbitals)+1
 866    read(funit,*) nmesh
 867    LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh))
 868    do imsh=1,nmesh
 869      rread2=zero
 870      read(funit,'(a80)') pspline
 871      read(unit=pspline,fmt=*,err=20,end=20) ii,iread1,iread2,rread1,rread2
 872      20 continue
 873      if (ii<=nmesh) then
 874        radmesh(ii)%mesh_type=iread1
 875        radmesh(ii)%mesh_size=iread2
 876        radmesh(ii)%rstep=rread1
 877        radmesh(ii)%lstep=rread2
 878      else
 879        write(msg, '(3a)' )&
 880 &       'Index of mesh out of range !',ch10,&
 881 &       'Action : check your pseudopotential file.'
 882        MSG_ERROR(msg)
 883      end if
 884    end do
 885    read(funit,*) pawtab%rpaw
 886    read(funit,'(a80)') pspline
 887    read(unit=pspline,fmt=*) pawtab%shape_type
 888    pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99
 889  end if
 890 
 891 !Initialize radial meshes
 892  do imsh=1,nmesh
 893    call pawrad_init(radmesh(imsh))
 894  end do
 895 
 896 !==========================================================
 897 !Initialize various dims and indexes
 898 
 899  pawtab%l_size=2*maxval(pawtab%orbitals)+1
 900  pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2
 901  pawtab%ij_size=pawtab%basis_size*(pawtab%basis_size+1)/2
 902  pawtab%usexcnhat=usexcnhat_in
 903 
 904 !indlmn calculation (indices for (l,m,n) basis)
 905  if (allocated(pawtab%indlmn)) then
 906    LIBPAW_DEALLOCATE(pawtab%indlmn)
 907  end if
 908  LIBPAW_ALLOCATE(pawtab%indlmn,(6,pawtab%lmn_size))
 909  LIBPAW_BOUND1_ALLOCATE(nprj,BOUNDS(0,maxval(pawtab%orbitals)))
 910  pawtab%indlmn(:,:)=0
 911  ilmn=0;iln=0;nprj=0
 912  do ib=1,pawtab%basis_size
 913    il=pawtab%orbitals(ib)
 914    nprj(il)=nprj(il)+1
 915    iln=iln+1
 916    do ilm=1,2*il+1
 917      pawtab%indlmn(1,ilmn+ilm)=il
 918      pawtab%indlmn(2,ilmn+ilm)=ilm-(il+1)
 919      pawtab%indlmn(3,ilmn+ilm)=nprj(il)
 920      pawtab%indlmn(4,ilmn+ilm)=il*il+ilm
 921      pawtab%indlmn(5,ilmn+ilm)=iln
 922      pawtab%indlmn(6,ilmn+ilm)=1
 923    end do
 924    ilmn=ilmn+2*il+1
 925  end do
 926  LIBPAW_DEALLOCATE(nprj)
 927 !Are ilmn (found here) and pawtab%lmn_size compatibles ?
 928  if (ilmn/=pawtab%lmn_size) then
 929    write(msg, '(a,a,a,a,a)' )&
 930 &   'Calculated lmn size differs from',ch10,&
 931 &   'lmn_size read from pseudo !',ch10,&
 932 &   'Action: check your pseudopotential file.'
 933    MSG_ERROR(msg)
 934  end if
 935 
 936 !==========================================================
 937 !Here reading shapefunction parameters
 938 
 939 !Shapefunction parameters for Abinit 4.3...4.5
 940  if (pspversion==2) then
 941    if (pawtab%shape_type==1) read(unit=pspline,fmt=*) ii,pawtab%shape_lambda,pawtab%shape_sigma
 942    if (pawtab%shape_type==3) pawtab%shape_type=-1
 943    pawtab%rshp=zero
 944 
 945 !Shapefunction parameters for Abinit 4.6+
 946  else if (pspversion>=3) then
 947    pawtab%rshp=zero
 948    if (pawtab%shape_type==-1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp
 949    if (pawtab%shape_type== 1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp, &
 950 &   pawtab%shape_lambda,pawtab%shape_sigma
 951    if (pawtab%shape_type== 2) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp
 952    if (pawtab%shape_type== 3) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp
 953  end if
 954  21 continue
 955 !If shapefunction type is gaussian, check exponent
 956  if (pawtab%shape_type==1) then
 957    if (pawtab%shape_lambda<2) then
 958      write(msg, '(3a)' )&
 959 &     'For a gaussian shape function, exponent lambda must be >1 !',ch10,&
 960 &     'Action: check your psp file.'
 961      MSG_ERROR(msg)
 962    end if
 963  end if
 964 !If shapefunction type is Bessel, deduce here its parameters from rc
 965  if (pawtab%shape_type==3) then
 966    LIBPAW_ALLOCATE(pawtab%shape_alpha,(2,pawtab%l_size))
 967    LIBPAW_ALLOCATE(pawtab%shape_q,(2,pawtab%l_size))
 968    rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw
 969    do il=1,pawtab%l_size
 970      call atompaw_shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc)
 971    end do
 972  end if
 973 
 974 !==========================================================
 975 !Mirror pseudopotential parameters to the output and log files
 976 
 977  write(msg,'(a,i1)')' Pseudopotential format is: paw',pspversion
 978  call wrtout(ab_out,msg,'COLL')
 979  call wrtout(std_out,  msg,'COLL')
 980  write(msg,'(2(a,i3),a,64i4)') &
 981 & ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',&
 982 & pawtab%lmn_size,'), orbitals=',pawtab%orbitals(1:pawtab%basis_size)
 983  call wrtout(ab_out,msg,'COLL')
 984  call wrtout(std_out,  msg,'COLL')
 985  write(msg,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw
 986  call wrtout(ab_out,msg,'COLL')
 987  call wrtout(std_out,  msg,'COLL')
 988  write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:'
 989  call wrtout(ab_out,msg,'COLL')
 990  call wrtout(std_out,  msg,'COLL')
 991  do imsh=1,nmesh
 992    if (radmesh(imsh)%mesh_type==1) &
 993 &   write(msg,'(a,i1,a,i4,a,g12.5)') &
 994 &   '  - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,&
 995 &   ' , step=',radmesh(imsh)%rstep
 996    if (radmesh(imsh)%mesh_type==2) &
 997 &   write(msg,'(a,i1,a,i4,2(a,g12.5))') &
 998 &   '  - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,&
 999 &   ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep
1000    if (radmesh(imsh)%mesh_type==3) &
1001 &   write(msg,'(a,i1,a,i4,2(a,g12.5))') &
1002 &   '  - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,&
1003 &   ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep
1004    if (radmesh(imsh)%mesh_type==4) &
1005 &   write(msg,'(a,i1,a,i4,a,g12.5)') &
1006 &   '  - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,&
1007 &   ' , AA=',radmesh(imsh)%rstep
1008    call wrtout(ab_out,msg,'COLL')
1009    call wrtout(std_out,  msg,'COLL')
1010  end do
1011  if (pawtab%shape_type==-1) then
1012    write(msg,'(a)')&
1013    ' Shapefunction is NUMERIC type: directly read from atomic data file'
1014    call wrtout(ab_out,msg,'COLL')
1015    call wrtout(std_out,  msg,'COLL')
1016  end if
1017  if (pawtab%shape_type==1) then
1018    write(msg,'(2a,a,f6.3,a,i3)')&
1019 &   ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,&
1020 &   '                            with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda
1021    call wrtout(ab_out,msg,'COLL')
1022    call wrtout(std_out,  msg,'COLL')
1023  end if
1024  if (pawtab%shape_type==2) then
1025    write(msg,'(a)')&
1026    ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2'
1027    call wrtout(ab_out,msg,'COLL')
1028    call wrtout(std_out,  msg,'COLL')
1029  end if
1030  if (pawtab%shape_type==3) then
1031    write(msg,'(a)')&
1032 &   ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)'
1033    call wrtout(ab_out,msg,'COLL')
1034    call wrtout(std_out,  msg,'COLL')
1035  end if
1036  if (pawtab%rshp<1.d-8) then
1037    write(msg,'(a)') ' Radius for shape functions = sphere core radius'
1038  else
1039    write(msg,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp
1040  end if
1041  call wrtout(ab_out,msg,'COLL')
1042  call wrtout(std_out,  msg,'COLL')
1043 
1044 !==========================================================
1045 !Perfom tests
1046 
1047 !Are lmax and orbitals compatibles ?
1048  if (lmax/=maxval(pawtab%orbitals)) then
1049    write(msg, '(a,a,a)' )&
1050 &   'lmax /= MAX(orbitals) !',ch10,&
1051 &   'Action: check your pseudopotential file.'
1052    MSG_ERROR(msg)
1053  end if
1054 
1055 !Only mesh_type=1,2, 3 or 4 allowed
1056  do imsh=1,nmesh
1057    if (radmesh(imsh)%mesh_type>4) then
1058      write(msg, '(a,a,a)' )&
1059 &     'Only mesh types 1,2, 3 or 4 allowed !',ch10,&
1060 &     'Action : check your pseudopotential or input file.'
1061      MSG_ERROR(msg)
1062    end if
1063  end do
1064 
1065 !==========================================================
1066 !Read tabulated atomic data
1067 
1068 !---------------------------------
1069 !Read wave-functions (phi)
1070  do ib=1,pawtab%basis_size
1071    read (funit,*)
1072    if (pspversion==1) iread1=1
1073    if (pspversion>1) read (funit,*) iread1
1074    if (ib==1) then
1075      call pawrad_free(pawrad)
1076      call pawrad_init(pawrad,mesh_size=radmesh(iread1)%mesh_size,mesh_type=radmesh(iread1)%mesh_type,&
1077 &     rstep=radmesh(iread1)%rstep,lstep=radmesh(iread1)%lstep,r_for_intg=pawtab%rpaw)
1078      pawtab%partialwave_mesh_size=pawrad%mesh_size
1079      pawtab%mesh_size=pawrad_ifromr(pawrad,pawtab%rpaw)+5
1080      pawtab%mesh_size=min(pawtab%mesh_size,pawrad%mesh_size)
1081      if (pawtab%mesh_size>pawrad%mesh_size-2) pawtab%mesh_size=pawrad%mesh_size
1082      imainmesh=iread1
1083      LIBPAW_ALLOCATE(pawtab%phi,(pawtab%partialwave_mesh_size,pawtab%basis_size))
1084    else if (iread1/=imainmesh) then
1085      write(msg, '(a,a,a)' )&
1086 &     'All Phi and tPhi must be given on the same radial mesh !',ch10,&
1087 &     'Action: check your pseudopotential file.'
1088      MSG_ERROR(msg)
1089    end if
1090    read (funit,*) (pawtab%phi(ir,ib),ir=1,pawtab%partialwave_mesh_size)
1091  end do
1092 
1093 !---------------------------------
1094 !Read pseudo wave-functions (tphi)
1095  LIBPAW_ALLOCATE(pawtab%tphi,(pawtab%partialwave_mesh_size,pawtab%basis_size))
1096  do ib=1,pawtab%basis_size
1097    read (funit,*)
1098    if (pspversion==1) iread1=1
1099    if (pspversion>1) read (funit,*) iread1
1100    if (iread1/=imainmesh) then
1101      write(msg, '(a,a,a)' )&
1102 &     'All Phi and tPhi must be given on the same radial mesh !',ch10,&
1103 &     'Action: check your pseudopotential file.'
1104      MSG_ERROR(msg)
1105    end if
1106    read (funit,*) (pawtab%tphi(ir,ib),ir=1,pawtab%partialwave_mesh_size)
1107  end do
1108  write(msg,'(a,i1)') &
1109 & ' Radial grid used for partial waves is grid ',imainmesh
1110  call wrtout(ab_out,msg,'COLL')
1111  call wrtout(std_out,  msg,'COLL')
1112 
1113 !---------------------------------
1114 !Read projectors (tproj)
1115  do ib=1,pawtab%basis_size
1116    read (funit,*)
1117    if (pspversion==1) iread1=2
1118    if (pspversion>1) read (funit,*) iread1
1119    if (ib==1) then
1120      iprojmesh=iread1
1121      call pawrad_copy(radmesh(iprojmesh),tproj_mesh)
1122      LIBPAW_POINTER_ALLOCATE(tproj,(tproj_mesh%mesh_size,pawtab%basis_size))
1123    else if (iread1/=iprojmesh) then
1124      write(msg, '(a,a,a)' )&
1125 &     'All tprojectors must be given on the same radial mesh !',ch10,&
1126 &     'Action: check your pseudopotential file.'
1127      MSG_ERROR(msg)
1128    end if
1129 !  read projectors from a mesh
1130    read (funit,*) (tproj(ir,ib),ir=1,tproj_mesh%mesh_size)
1131  end do
1132  write(msg,'(a,i2)') &
1133 & ' Radial grid used for projectors is grid ',iprojmesh
1134  call wrtout(ab_out,msg,'COLL')
1135  call wrtout(std_out,  msg,'COLL')
1136 
1137 !---------------------------------
1138 !Read gaussian projectors for wavelets
1139 !  -- only if pawtab%has_wvl flag is on
1140 !  -- if not, we skip the lines
1141  read(funit,'(a80)') pspline
1142  if(index(trim(pspline),'GAUSSIAN')/=0) read_gauss=.true.
1143  if (read_gauss) then
1144    if (pawtab%has_wvl>0) then
1145      call wvlpaw_allocate(pawtab%wvl)
1146      jj=0
1147      do ib=1,pawtab%basis_size
1148        if(ib/=1) read(funit,*) pspline
1149 !      read Gaussian coefficients
1150        read(funit,*) pngau_, ptotgau_ !total number of gaussians
1151        if(ib==1) then
1152          pawtab%wvl%ptotgau=ptotgau_
1153          LIBPAW_ALLOCATE(pawtab%wvl%pngau,(pawtab%basis_size))
1154          LIBPAW_ALLOCATE(pawtab%wvl%parg,(2,pawtab%wvl%ptotgau))
1155          LIBPAW_ALLOCATE(pawtab%wvl%pfac,(2,pawtab%wvl%ptotgau))
1156        else
1157          if(pawtab%wvl%ptotgau/=ptotgau_) then
1158            write(msg,'(3a)')&
1159 &           'Total number of gaussians, should be the same for all projectors !',ch10,&
1160 &           'Action: check your pseudopotential file.'
1161            MSG_ERROR(msg)
1162          end if
1163        end if !ib==1
1164        read(funit,*)(pawtab%wvl%parg(:,ii),ii=jj+1,jj+pngau_)
1165        read(funit,*)(pawtab%wvl%pfac(:,ii),ii=jj+1,jj+pngau_)
1166        pawtab%wvl%pngau(ib)=pngau_
1167        jj=jj+pngau_
1168      end do
1169      pawtab%has_wvl=2
1170    else
1171 !    If pawtab%has_wvl=0, we skip the lines
1172      do ib=1,pawtab%basis_size
1173        if(ib/=1) read(funit,*)
1174        read(funit,*) pngau_, ptotgau_
1175        LIBPAW_ALLOCATE(val, (pngau_  *2))
1176        read(funit,*) val
1177        read(funit,*) val
1178        LIBPAW_DEALLOCATE(val)
1179      end do
1180    end if
1181  end if
1182 
1183 !---------------------------------
1184 !Read core density (coredens)
1185  if(read_gauss) read (funit,*) !if not read_gauss, this line was already read
1186  if (pspversion==1) iread1=1
1187  if (pspversion>1) read (funit,*) iread1
1188  icoremesh=iread1
1189  call pawrad_copy(radmesh(icoremesh),core_mesh)
1190  if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.&
1191 & (radmesh(icoremesh)%rstep    /=pawrad%rstep)    .or.&
1192 & (radmesh(icoremesh)%lstep    /=pawrad%lstep)) then
1193    write(msg, '(a,a,a,a,a)' )&
1194 &   'Ncore must be given on a radial mesh with the same',ch10,&
1195 &   'type and step(s) than the main radial mesh (mesh for Phi) !',ch10,&
1196 &   'Action: check your pseudopotential file.'
1197    MSG_ERROR(msg)
1198  end if
1199  LIBPAW_POINTER_ALLOCATE(ncore,(core_mesh%mesh_size))
1200  read (funit,*) (ncore(ir),ir=1,core_mesh%mesh_size)
1201 
1202 !Construct and save VH[z_NC] if requested
1203  if (pawtab%has_vhnzc==1) then
1204    LIBPAW_ALLOCATE(pawtab%VHnZC,(pawtab%mesh_size))
1205    LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size))
1206    call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl)
1207    pawtab%VHnZC(1:pawtab%mesh_size)=vhnzc(1:pawtab%mesh_size)
1208    pawtab%has_vhnzc=2
1209    LIBPAW_DEALLOCATE(vhnzc)
1210  end if
1211 
1212  pawtab%core_mesh_size=pawrad%mesh_size
1213  if(save_core_msz) pawtab%core_mesh_size=core_mesh%mesh_size
1214  LIBPAW_ALLOCATE(pawtab%coredens,(pawtab%core_mesh_size))
1215  pawtab%rcore=core_mesh%rad(pawtab%core_mesh_size)
1216  pawtab%coredens(1:pawtab%core_mesh_size)=ncore(1:pawtab%core_mesh_size)
1217 
1218 !---------------------------------
1219 !Read pseudo core density (tcoredens)
1220  if(save_core_msz)  then
1221    LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,6))
1222  else
1223    LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,1))
1224  end if
1225  pawtab%tcoredens=zero
1226  read (funit,*)
1227  if (pspversion==1) iread1=1
1228  if (pspversion>1) read (funit,*) iread1
1229  if (iread1/=icoremesh) then
1230    write(msg, '(a,a,a,a,a,a,a,a)' )&
1231 &   'Pseudized core density (tNcore) must be given',ch10,&
1232 &   'on the same radial mesh as core density (Ncore) !',ch10,&
1233 &   'Action: check your pseudopotential file.'
1234    MSG_ERROR(msg)
1235  end if
1236  LIBPAW_POINTER_ALLOCATE(tncore,(core_mesh%mesh_size))
1237  read (funit,*) (tncore(ir),ir=1,core_mesh%mesh_size)
1238  if (maxval(abs(tncore(:)))<tol6) then
1239    pawtab%usetcore=0
1240  else
1241    pawtab%usetcore=1
1242    pawtab%tcoredens(1:pawtab%core_mesh_size,1)=tncore(1:pawtab%core_mesh_size)
1243  end if
1244  write(msg,'(a,i1)') &
1245 & ' Radial grid used for (t)core density is grid ',icoremesh
1246  call wrtout(ab_out,msg,'COLL')
1247  call wrtout(std_out,  msg,'COLL')
1248 
1249 !---------------------------------
1250 !Read frozen part of Dij terms (dij0)
1251  LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size))
1252  read (funit,*)
1253  read (funit,*) (pawtab%dij0(ib),ib=1,pawtab%lmn2_size)
1254 
1255 !---------------------------------
1256 !Read initial guess of rhoij (rhoij0)
1257  LIBPAW_ALLOCATE(pawtab%rhoij0,(pawtab%lmn2_size))
1258  read (funit,*)
1259  read (funit,*) (pawtab%rhoij0(ib),ib=1,pawtab%lmn2_size)
1260 
1261 !---------------------------------
1262 !Read local pseudopotential=Vh(tn_zc) or Vbare
1263  read (funit,*)
1264  if (pspversion==1) ivlocmesh=3
1265  vlocopt=1
1266  if (pspversion==2) then
1267    read (funit,*) ivlocmesh
1268  else if (pspversion>2) then
1269 !  read (funit,fmt=*,err=30,end=30) ivlocmesh,vlocopt
1270    msg=blank
1271    read (funit,fmt='(a)') msg
1272    read (msg,fmt=*) ivlocmesh
1273    write(numb,'(i1)')ivlocmesh
1274    ii=index(msg,numb)
1275    if(len_trim(trim(msg(ii+1:)))/=0)then
1276      submsg=trim(msg(ii+1:))
1277      if(len_trim(submsg)/=0)then
1278        do ii=1,len_trim(submsg)
1279          numb=submsg(ii:ii)
1280          if(numb==blank)cycle
1281          jj=index('0123456789',numb)
1282          if(jj<1 .or. jj>10)exit
1283          vlocopt=jj-1
1284        end do
1285      end if
1286    end if
1287  end if
1288  usexcnhat_out=0;if (vlocopt==1) usexcnhat_out=1
1289  call pawrad_copy(radmesh(ivlocmesh),vloc_mesh)
1290  LIBPAW_POINTER_ALLOCATE(vlocr,(vloc_mesh%mesh_size))
1291  read (funit,*) (vlocr(ir),ir=1,vloc_mesh%mesh_size)
1292  write(msg,'(a,i1)') &
1293 & ' Radial grid used for Vloc is grid ',ivlocmesh
1294  call wrtout(ab_out,msg,'COLL')
1295  call wrtout(std_out,  msg,'COLL')
1296 
1297 !---------------------------------
1298 !Eventually read "numeric" shapefunctions (if shape_type=-1)
1299  if (pawtab%shape_type==-1) then
1300    LIBPAW_ALLOCATE(pawtab%shapefunc,(pawtab%mesh_size,pawtab%l_size))
1301    do il=1,pawtab%l_size
1302      read (funit,*)
1303      if (pspversion==1) iread1=1
1304      if (pspversion>1) read (funit,*) iread1
1305      if (il==1) then
1306        call pawrad_copy(radmesh(iread1),shpf_mesh)
1307        ishpfmesh=iread1
1308        LIBPAW_ALLOCATE(shpf,(shpf_mesh%mesh_size,pawtab%l_size))
1309      else if (iread1/=ishpfmesh) then
1310        write(msg, '(a,a,a)' )&
1311 &       'All shape functions must be given on the same radial mesh !',ch10,&
1312 &       'Action: check your pseudopotential file.'
1313        MSG_ERROR(msg)
1314      end if
1315      read (funit,*) (shpf(ir,il),ir=1,shpf_mesh%mesh_size)
1316    end do
1317    write(msg,'(a,i1)') &
1318 &   ' Radial grid used for shape functions is grid ',iread1
1319    call wrtout(ab_out,msg,'COLL')
1320    call wrtout(std_out,  msg,'COLL')
1321 
1322 !  Has to spline shape functions if mesh is not the "main" mesh
1323    if (ishpfmesh/=imainmesh) then
1324      msz=shpf_mesh%mesh_size
1325      LIBPAW_ALLOCATE(work1,(msz))
1326      LIBPAW_ALLOCATE(work2,(msz))
1327      LIBPAW_ALLOCATE(work3,(msz))
1328      LIBPAW_ALLOCATE(work4,(pawtab%mesh_size))
1329      work3(1:pawtab%mesh_size)=shpf_mesh%rad(1:pawtab%mesh_size)
1330      work4(1:pawtab%mesh_size)=pawrad%rad(1:pawtab%mesh_size)
1331      do il=1,pawtab%l_size
1332        call bound_deriv(shpf(1:msz,il),shpf_mesh,msz,yp1,ypn)
1333        call paw_spline(work3,shpf(:,il),msz,yp1,ypn,work1)
1334        call paw_splint(msz,work3,shpf(:,il),work1,pawtab%mesh_size,work4,pawtab%shapefunc(:,il))
1335      end do
1336      LIBPAW_DEALLOCATE(work1)
1337      LIBPAW_DEALLOCATE(work2)
1338      LIBPAW_DEALLOCATE(work3)
1339      LIBPAW_DEALLOCATE(work4)
1340    else
1341      pawtab%shapefunc(:,:)=shpf(:,:)
1342    end if
1343    LIBPAW_DEALLOCATE(shpf)
1344    call pawrad_free(shpf_mesh)
1345  end if
1346 
1347 !---------------------------------
1348 !Read pseudo valence density (if psp version >=4)
1349  if (pspversion>=4) then
1350    read (funit,*)
1351    read (funit,*) iread1
1352    ivalemesh=iread1
1353    call pawrad_copy(radmesh(iread1),vale_mesh)
1354    LIBPAW_POINTER_ALLOCATE(tnvale,(vale_mesh%mesh_size))
1355    read (funit,*) (tnvale(ir),ir=1,vale_mesh%mesh_size)
1356    pawtab%has_tvale=1
1357    write(msg,'(a,i1)') &
1358 &   ' Radial grid used for pseudo valence density is grid ',ivalemesh
1359    call wrtout(ab_out,msg,'COLL')
1360    call wrtout(std_out,  msg,'COLL')
1361  else
1362    pawtab%has_tvale=0
1363    LIBPAW_POINTER_ALLOCATE(tnvale,(0))
1364  end if
1365 
1366 end subroutine pawpsp_read

m_pawpsp/pawpsp_read_corewf [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_read_corewf

FUNCTION

INPUTS

 [filename]= (optional) core WF file name

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      optics_paw_core,posdoppler

CHILDREN

SOURCE

1392 subroutine pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,radmesh,phi_cor,&
1393 &                             filename) ! optional argument
1394 
1395 
1396 !This section has been created automatically by the script Abilint (TD).
1397 !Do not modify the following lines by hand.
1398 #undef ABI_FUNC
1399 #define ABI_FUNC 'pawpsp_read_corewf'
1400 !End of the abilint section
1401 
1402  implicit none
1403 
1404 !Arguments ------------------------------------
1405  integer,intent(out) :: lmncmax,nphicor
1406  character(len=*),optional :: filename
1407 !arrays
1408  integer,allocatable,intent(inout) :: indlmn_core(:,:),lcor(:),ncor(:)
1409  real(dp),allocatable,intent(inout) :: phi_cor(:,:),energy_cor(:)
1410  type(pawrad_type),intent(in) :: radmesh
1411 
1412 !Local variables-------------------------------
1413  integer :: ib,i1,i2,il,ilm,ilmn,iln,ios,jln,nmesh,npts,unt
1414  real(dp) :: noccor,r1,r2
1415  logical :: ex,oldformat,usexml
1416  character(len=8) :: dum,dum1,dum2,dum3,dum4
1417  character(len=80) :: fline
1418  character(len=500) :: filename_,msg
1419 !arrays
1420  integer,allocatable :: meshtp(:),meshsz(:)
1421  real(dp),allocatable :: rad(:),radstp(:),work(:)
1422  real(dp),allocatable :: logstp(:),phitmp(:)
1423  type(pawrad_type) :: tmpmesh
1424 
1425 ! ************************************************************************
1426 
1427 !Check for core WF file existence and XML format
1428  usexml=.false.;oldformat=.false.
1429  if (present(filename)) then
1430 !  Core WF file given as optional argument
1431    filename_=filename;ex=.false.
1432    inquire(file=trim(filename_),iostat=ios,exist=ex)
1433    if (ios/=0) then
1434      write(msg,'(2a)') 'INQUIRE returns an error for file ',trim(filename_)
1435      MSG_ERROR(msg)
1436    end if
1437    if (.not.ex) then
1438      write(msg,'(3a)') 'This file does not exist: ',trim(filename_),'!'
1439      MSG_ERROR(msg)
1440    end if
1441    unt = libpaw_get_free_unit()
1442    open(unit=unt,file=trim(filename_),form='formatted',status='old', action="read")
1443    read(unt,*) fline
1444    close(unt)
1445    usexml=(fline(1:5)=='<?xml')
1446  else
1447 !  Core WF file: new format
1448    filename_='corewf.abinit';ex=.false.
1449    inquire(file=trim(filename_),iostat=ios,exist=ex)
1450    if (ios/=0) then
1451      write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!'
1452      MSG_ERROR(msg)
1453    end if
1454    if (.not.ex) then
1455 !    Core WF file: new format XML
1456      filename_='corewf.abinit.xml';ex=.false.
1457      inquire(file=trim(filename_),iostat=ios,exist=ex)
1458      if (ios/=0) then
1459        write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!'
1460        MSG_ERROR(msg)
1461      end if
1462      usexml=ex
1463      if (.not.ex) then
1464 !      Core WF file: old format
1465        filename_='corewf.dat';ex=.false.
1466        inquire(file=trim(filename_),iostat=ios,exist=ex)
1467        if (ios/=0) then
1468          write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!'
1469          MSG_ERROR(msg)
1470        end if
1471        oldformat=ex
1472        if (.not.ex) then
1473 !        No core WF file found
1474          write(msg, '(3a)' )&
1475 &         'Checks for existence of files corewf.abinit[.xml] or corewf.dat',ch10,&
1476 &         'but INQUIRE finds file does not exist!'
1477          MSG_ERROR(msg)
1478        end if
1479      end if
1480    end if
1481  end if
1482 
1483 !Core WF file is in new XML format
1484  if ((.not.oldformat).and.(usexml)) then
1485    call rdpawpsxml_core(energy_cor,trim(filename_),lcor,ncor,nphicor,radmesh,phi_cor)
1486  endif
1487 
1488 !Core WF file is in new (proprietary) format
1489  if ((.not.oldformat).and.(.not.usexml)) then
1490    unt = libpaw_get_free_unit()
1491    open(unt,file=trim(filename_),form='formatted',action="read")
1492    read(unt,*) ! skip title
1493    read(unt,*) ! skip method,nspinor,nsppol
1494    read(unt,*) ! skip zatom,zcore,pspdat
1495    read(unt,*) ! skip pspcod,pspxc,lmax
1496    read(unt,*) ! skip pspfmt,creatorID
1497    read(unt,*) nphicor
1498    read(unt,*) ! skip orbitals
1499    read(unt,*) nmesh
1500    LIBPAW_ALLOCATE(meshsz,(nmesh))
1501    LIBPAW_ALLOCATE(meshtp,(nmesh))
1502    LIBPAW_ALLOCATE(radstp,(nmesh))
1503    LIBPAW_ALLOCATE(logstp,(nmesh))
1504    do iln=1,nmesh
1505      r2=zero;read(unt,'(a80)') fline
1506      read(unit=fline,fmt=*,err=20,end=20) ib,i1,i2,r1,r2
1507      20 continue
1508      if (ib<=nmesh) then
1509        meshtp(ib)=i1;meshsz(ib)=i2
1510        radstp(ib)=r1;logstp(ib)=r2
1511      end if
1512    end do
1513    read(unt,*) ! skip rmax(core)
1514    LIBPAW_ALLOCATE(ncor,(nphicor))
1515    LIBPAW_ALLOCATE(lcor,(nphicor))
1516    LIBPAW_ALLOCATE(energy_cor,(nphicor))
1517    LIBPAW_ALLOCATE(phi_cor,(radmesh%mesh_size,nphicor))
1518    do iln=1,nphicor
1519      read(unt,*) ! skip comment
1520      read(unt,*) i1
1521      read(unt,*) ncor(iln),lcor(iln)
1522      read(unt,*) energy_cor(iln)
1523      LIBPAW_ALLOCATE(phitmp,(meshsz(i1)))
1524      read(unt,*) phitmp
1525      if ((radmesh%mesh_type/=meshtp(i1)) &
1526 &     .or.(radmesh%rstep/=radstp(i1)) &
1527 &     .or.(radmesh%lstep/=logstp(i1))) then
1528        call pawrad_init(tmpmesh,mesh_size=meshsz(i1),mesh_type=meshtp(i1),rstep=radstp(i1),lstep=logstp(i1))
1529        npts=radmesh%mesh_size
1530        if (tmpmesh%rmax<radmesh%rmax+tol8) npts=pawrad_ifromr(radmesh,tmpmesh%rmax)-1
1531        LIBPAW_ALLOCATE(work,(meshsz(i1)))
1532        call bound_deriv(phitmp,tmpmesh,meshsz(i1),r1,r2)
1533        call paw_spline(tmpmesh%rad,phitmp,meshsz(i1),r1,r2,work)
1534        call paw_splint(meshsz(i1),tmpmesh%rad,phitmp,work,npts,radmesh%rad(1:npts),phi_cor(1:npts,iln))
1535        if (npts<radmesh%mesh_size) phi_cor(npts+1:radmesh%mesh_size,iln)=zero
1536        LIBPAW_DEALLOCATE(work)
1537        call pawrad_free(tmpmesh)
1538      else
1539        npts=min(meshsz(i1),radmesh%mesh_size)
1540        phi_cor(1:npts,iln)=phitmp(1:npts)
1541        if (npts<radmesh%mesh_size) phi_cor(npts+1:radmesh%mesh_size,iln)=zero
1542      end if
1543      LIBPAW_DEALLOCATE(phitmp)
1544    end do
1545    LIBPAW_DEALLOCATE(meshsz)
1546    LIBPAW_DEALLOCATE(meshtp)
1547    LIBPAW_DEALLOCATE(radstp)
1548    LIBPAW_DEALLOCATE(logstp)
1549    close(unt)
1550  end if
1551 
1552 !Core WF file is in old (proprietary) format
1553  if ((oldformat).and.(.not.usexml)) then
1554    unt = libpaw_get_free_unit()
1555    open(unt,file=trim(filename_),form='formatted',action="read")
1556    do while (dum/='atompaw ')
1557      read(unt,'(a8)') dum
1558    end do
1559    read(unt,'(2i4)') npts,nphicor
1560    LIBPAW_ALLOCATE(ncor,(nphicor))
1561    LIBPAW_ALLOCATE(lcor,(nphicor))
1562    LIBPAW_ALLOCATE(energy_cor,(nphicor))
1563    LIBPAW_ALLOCATE(phi_cor,(npts,nphicor))
1564    LIBPAW_ALLOCATE(rad,(npts))
1565    do iln=1,nphicor
1566      read(unt,'(a4,i4,a3,i4,a6,f15.7,a8,f15.7)') &
1567 &     dum1,ncor(iln),dum2,lcor(iln),dum3,noccor,dum4,energy_cor(iln)
1568 
1569      do jln=1,npts
1570        read(unt,*) rad(jln),phi_cor(jln,iln)
1571      end do
1572      read(unt,*)
1573    end do
1574    LIBPAW_DEALLOCATE(rad)
1575    close(unt)
1576  end if
1577 
1578 !Set an array 'a la' indlmn
1579  lmncmax=0
1580  do ib=1,nphicor
1581    il=lcor(ib)
1582    lmncmax=lmncmax+2*il+1
1583  end do
1584  LIBPAW_ALLOCATE(indlmn_core,(6,lmncmax))
1585  indlmn_core=0;ilmn=0;iln=0
1586  do ib=1,nphicor
1587    il=lcor(ib)
1588    iln=iln+1
1589    do ilm=1,2*il+1
1590      indlmn_core(1,ilmn+ilm)=il
1591      indlmn_core(2,ilmn+ilm)=ilm-(il+1)
1592      indlmn_core(3,ilmn+ilm)=1
1593      indlmn_core(4,ilmn+ilm)=il*il+ilm
1594      indlmn_core(5,ilmn+ilm)=iln
1595      indlmn_core(6,ilmn+ilm)=1
1596    end do
1597    ilmn=ilmn+2*il+1
1598  end do
1599 
1600 end subroutine pawpsp_read_corewf

m_pawpsp/pawpsp_read_header [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_read_header

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

4089 subroutine pawpsp_read_header(funit,lloc,lmax,mmax,pspcod,pspxc,r2well,zion,znucl)
4090 
4091 
4092 !This section has been created automatically by the script Abilint (TD).
4093 !Do not modify the following lines by hand.
4094 #undef ABI_FUNC
4095 #define ABI_FUNC 'pawpsp_read_header'
4096 !End of the abilint section
4097 
4098 implicit none
4099 !Arguments ------------------------------------
4100 !scalars
4101  integer,intent(in):: funit
4102  integer,intent(out):: lloc,lmax,mmax,pspcod,pspxc
4103  real(dp),intent(out):: r2well,zion,znucl
4104 !Local variables-------------------------------
4105  integer:: pspdat
4106  character(len=fnlen):: title
4107  character(len=500) :: msg
4108 
4109 ! *************************************************************************
4110 
4111 !Read and write some description of file from first line (character data)
4112  read (funit,'(a)') title
4113  write(msg, '(a,a)' ) '- ',trim(title)
4114  call wrtout(ab_out,msg,'COLL')
4115  call wrtout(std_out,  msg,'COLL')
4116 
4117 !Read and write more data describing psp parameters
4118  read (funit,*) znucl,zion,pspdat
4119  write(msg, '(a,f9.5,f10.5,2x,i8,t47,a)' ) &
4120 & '-',znucl,zion,pspdat,'znucl, zion, pspdat'
4121  call wrtout(ab_out,msg,'COLL')
4122  call wrtout(std_out,  msg,'COLL')
4123 
4124  read (funit,*) pspcod,pspxc,lmax,lloc,mmax,r2well
4125  if(pspxc<0) then
4126    write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) &
4127 &   pspcod,pspxc,lmax,lloc,mmax,r2well,&
4128 &   'pspcod,pspxc,lmax,lloc,mmax,r2well'
4129  else
4130    write(msg, '(4i5,i10,f10.5,t47,a)' ) &
4131 &   pspcod,pspxc,lmax,lloc,mmax,r2well,&
4132 &   'pspcod,pspxc,lmax,lloc,mmax,r2well'
4133  end if
4134  call wrtout(ab_out,msg,'COLL')
4135  call wrtout(std_out,  msg,'COLL')
4136 
4137 end subroutine pawpsp_read_header

m_pawpsp/pawpsp_read_header_2 [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_read_header_2

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 Reads pspversion, basis_size and lmn_size

PARENTS

      m_pawpsp

CHILDREN

SOURCE

4166 subroutine pawpsp_read_header_2(funit,pspversion,basis_size,lmn_size)
4167 
4168 
4169 !This section has been created automatically by the script Abilint (TD).
4170 !Do not modify the following lines by hand.
4171 #undef ABI_FUNC
4172 #define ABI_FUNC 'pawpsp_read_header_2'
4173 !End of the abilint section
4174 
4175 implicit none
4176 
4177 !Arguments ------------------------------------
4178 !scalars
4179  integer,intent(in):: funit
4180  integer,intent(out) :: pspversion,basis_size,lmn_size
4181 
4182 !Local variables-------------------------------
4183  integer :: creatorid
4184  character(len=80) :: pspline
4185  character(len=500) :: msg
4186 
4187 ! *************************************************************************
4188 
4189 !Read psp version in line 4 of the header
4190  pspversion=1
4191  read (funit,'(a80)') pspline;pspline=adjustl(pspline)
4192  if (pspline(1:3)=="paw".or.pspline(1:3)=="PAW") &
4193 & read(unit=pspline(4:80),fmt=*) pspversion
4194  if (pspversion<1.or.pspversion>5) then
4195    write(msg, '(a,i2,a,a,a)' )&
4196 &   'This version of PAW psp file (',pspversion,') is not compatible with',ch10,&
4197 &   'current version of Abinit.'
4198    MSG_ERROR(msg)
4199  end if
4200 
4201  if (pspversion==1) then
4202    read (unit=pspline,fmt=*) basis_size,lmn_size
4203  else
4204 !  Here psp file for Abinit 4.3+
4205    read (unit=pspline(5:80),fmt=*) creatorid
4206    read (funit,*) basis_size,lmn_size
4207  end if
4208 
4209 end subroutine pawpsp_read_header_2

m_pawpsp/pawpsp_read_header_xml [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_read_header_xml

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

 This is done instead of: call pawpsxml2ab( psxml, pspheads,1)
 since pspheads does not exist in PAW library.
 should we include it to avoid the following code replica?
 check pspheads commented out in pawpsp_17in, and routine pawpsp_read_xml_2

PARENTS

      m_pawpsp,pawpsxml2ab,pspatm

CHILDREN

SOURCE

4362 subroutine pawpsp_read_header_xml(lloc,lmax,pspcod,pspxc,&
4363 & psxml,r2well,zion,znucl)
4364 
4365 
4366 !This section has been created automatically by the script Abilint (TD).
4367 !Do not modify the following lines by hand.
4368 #undef ABI_FUNC
4369 #define ABI_FUNC 'pawpsp_read_header_xml'
4370 !End of the abilint section
4371 
4372 implicit none
4373 !Arguments ------------------------------------
4374 !scalars
4375  type(paw_setup_t),intent(in) :: psxml
4376  integer,intent(out):: lloc,lmax,pspcod,pspxc
4377  real(dp),intent(out):: r2well,zion,znucl
4378 !Local variables-------------------------------
4379  integer :: il
4380 #if defined LIBPAW_HAVE_LIBXC
4381  integer :: ii
4382 #endif
4383  character(len=100) :: xclibxc
4384  character(len=500) :: msg
4385 !arrays
4386 
4387 ! *************************************************************************
4388 
4389  lloc   = 0
4390  r2well = 0
4391  pspcod=17
4392  znucl=psxml%atom%znucl
4393  zion =psxml%atom%zval
4394 
4395 !lmax:
4396  lmax = 0
4397  do il=1,psxml%valence_states%nval
4398    if(psxml%valence_states%state(il)%ll>lmax) lmax=psxml%valence_states%state(il)%ll
4399  end do
4400 !pspxc
4401  select case(trim(psxml%xc_functional%name))
4402    case('PZ')
4403      pspxc = 2
4404 #if defined LIBPAW_HAVE_LIBXC
4405      pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 &
4406 &             +libxc_functionals_getid('XC_LDA_C_PZ'))
4407 #endif
4408    case('W')
4409      pspxc = 4
4410 #if defined LIBPAW_HAVE_LIBXC
4411      pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 &
4412 &             +libxc_functionals_getid('XC_LDA_C_WIGNER'))
4413 #endif
4414    case('HL')
4415      pspxc = 5
4416 #if defined LIBPAW_HAVE_LIBXC
4417      pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 &
4418 &             +libxc_functionals_getid('XC_LDA_C_HL'))
4419 #endif
4420    case('GL')
4421 #if defined LIBPAW_HAVE_LIBXC
4422      pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 &
4423 &             +libxc_functionals_getid('XC_LDA_C_GL'))
4424 #else
4425      write(msg, '(7a)' )&
4426 &     'The exchange and correlation functional by Gunnarson-Lundqvist', ch10,&
4427 &     'is not implemented in Abinit.',ch10,&
4428 &     'Action : choose another XC functional in the pseudopotential',ch10, &
4429 &     '         generation or compile ABINIT with the libXC library.'
4430      MSG_ERROR(msg)
4431 #endif
4432    case('VWN')
4433 #if defined LIBPAW_HAVE_LIBXC
4434      pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 &
4435 &             +libxc_functionals_getid('XC_LDA_C_VWN'))
4436 #else
4437      write(msg, '(7a)' )&
4438 &     'The exchange and correlation functional by Vosko,Wilk and Nusair', ch10,&
4439 &     'is not implemented in Abinit.',ch10,&
4440 &     'Action : choose another XC functional in the pseudopotential',ch10, &
4441 &     '         generation or compile ABINIT with the libXC library.'
4442      MSG_ERROR(msg)
4443 #endif
4444    case('PW')
4445      pspxc = 7
4446 #if defined LIBPAW_HAVE_LIBXC
4447      pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 &
4448 &             +libxc_functionals_getid('XC_LDA_C_PW'))
4449 #endif
4450    case('PBE')
4451      pspxc = 11
4452 #if defined LIBPAW_HAVE_LIBXC
4453      pspxc = -(libxc_functionals_getid('XC_GGA_X_PBE')*1000 &
4454 &             +libxc_functionals_getid('XC_GGA_C_PBE'))
4455 #endif
4456    case('revPBE')
4457      pspxc = 14
4458 #if defined LIBPAW_HAVE_LIBXC
4459      pspxc = -(libxc_functionals_getid('XC_GGA_X_PBE_R')*1000 &
4460 &             +libxc_functionals_getid('XC_GGA_C_PBE'))
4461 #endif
4462    case('RPBE')
4463      pspxc = 15
4464 #if defined LIBPAW_HAVE_LIBXC
4465      pspxc = -(libxc_functionals_getid('XC_GGA_X_RPBE')*1000 &
4466 &             +libxc_functionals_getid('XC_GGA_C_PBE'))
4467 #endif
4468    case('PW91')
4469 #if defined LIBPAW_HAVE_LIBXC
4470      pspxc = -(libxc_functionals_getid('XC_GGA_X_PW91')*1000 &
4471 &             +libxc_functionals_getid('XC_GGA_C_PW91'))
4472 #else
4473      write(msg, '(7a)' )&
4474 &     'The exchange and correlation functional by Perdew and Wang 91', ch10,&
4475 &     'is not implemented in Abinit.',ch10,&
4476 &     'Action : choose another XC functional in the pseudopotential',ch10, &
4477 &     '         generation or compile ABINIT with the libXC library.'
4478      MSG_ERROR(msg)
4479 #endif
4480    case('BLYP')
4481 #if defined LIBPAW_HAVE_LIBXC
4482      pspxc = -(libxc_functionals_getid('XC_GGA_X_B88')*1000 &
4483 &             +libxc_functionals_getid('XC_GGA_C_LYP'))
4484 #else
4485      write(msg, '(7a)' )&
4486 &     'The exchange and correlation functional BLYP', ch10,&
4487 &     'is not implemented in Abinit.',ch10,&
4488 &     'Action : choose another XC functional in the pseudopotential',ch10, &
4489 &     '         generation or compile ABINIT with the libXC library.'
4490      MSG_ERROR(msg)
4491 #endif
4492    case DEFAULT
4493      xclibxc=trim(psxml%xc_functional%name)
4494      if (xclibxc(1:3)=='XC_'  .or.xclibxc(1:3)=='xc_'  .or. &
4495 &        xclibxc(1:5)=='LDA_X'.or.xclibxc(1:5)=='LDA_C'.or. &
4496 &        xclibxc(1:5)=='lda_x'.or.xclibxc(1:5)=='lda_c'.or. &
4497 &        xclibxc(1:5)=='GGA_X'.or.xclibxc(1:5)=='GGA_C'.or. &
4498 &        xclibxc(1:5)=='gga_x'.or.xclibxc(1:5)=='gga_c') then
4499 #if defined LIBPAW_HAVE_LIBXC
4500        ii=index(xclibxc,'+')
4501        if (ii>0) then
4502          pspxc=-(libxc_functionals_getid(xclibxc(1:ii-1))*1000 &
4503 &               +libxc_functionals_getid(xclibxc(ii+1:)))
4504        else
4505          pspxc=-libxc_functionals_getid(xclibxc)
4506        end if
4507 #else
4508        msg='Cannot use LibXC functional because ABINIT is not compiled with LibXC !'
4509        MSG_ERROR(msg)
4510 #endif
4511 !      To be eliminated later (temporary)
4512      else if(trim(psxml%xc_functional%functionaltype)=='LIBXC')then
4513 #if defined LIBPAW_HAVE_LIBXC
4514        xclibxc=trim(psxml%xc_functional%name)
4515        read(unit=xclibxc,fmt=*) pspxc
4516        pspxc=-pspxc
4517 #else
4518        msg='Cannot use LibXC functional because ABINIT is not compiled with LibXC !'
4519        MSG_ERROR(msg)
4520 #endif
4521      else
4522        write(msg, '(3a)') 'Unknown XC functional in psp file: ',trim(xclibxc),' !'
4523        MSG_ERROR(msg)
4524      end if
4525  end select
4526 
4527 end subroutine pawpsp_read_header_xml

m_pawpsp/pawpsp_read_pawheader [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_read_pawheader

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp,pawpsxml2ab,pspatm

CHILDREN

SOURCE

4553 subroutine pawpsp_read_pawheader(basis_size,lmax,lmn_size,&
4554 & l_size,mesh_size,pspversion,psxml,rpaw,rshp,shape_type)
4555 
4556 
4557 !This section has been created automatically by the script Abilint (TD).
4558 !Do not modify the following lines by hand.
4559 #undef ABI_FUNC
4560 #define ABI_FUNC 'pawpsp_read_pawheader'
4561 !End of the abilint section
4562 
4563 implicit none
4564 !Arguments ------------------------------------
4565 !scalars
4566  integer,intent(in):: lmax
4567  integer,intent(out):: basis_size,mesh_size,lmn_size,l_size
4568  integer,intent(out):: pspversion,shape_type
4569  real(dp),intent(out)::rpaw,rshp
4570  type(paw_setup_t),intent(in) :: psxml
4571 !Local variables-------------------------------
4572  integer::il
4573 
4574 ! *************************************************************************
4575 
4576 !All of this was moved from pawpsxml2ab,
4577 !basis_size
4578  basis_size=psxml%valence_states%nval
4579 !mesh_size
4580  do il=1,psxml%ngrid
4581    if(psxml%radial_grid(il)%id==psxml%idgrid) &
4582 &   mesh_size=psxml%radial_grid(il)%iend-psxml%radial_grid(il)%istart+1
4583  end do
4584 !lmn_size:
4585  lmn_size=0
4586  do il=1,psxml%valence_states%nval
4587    lmn_size=lmn_size+2*psxml%valence_states%state(il)%ll+1
4588  end do
4589 !lsize
4590  l_size=2*lmax+1
4591 !pspversion
4592  pspversion=10
4593 !rpaw:
4594  rpaw=0.d0
4595  if (psxml%rpaw<0.d0) then
4596    do il=1,psxml%valence_states%nval
4597      if(psxml%valence_states%state(il)%rc>rpaw) rpaw=psxml%valence_states%state(il)%rc
4598    end do
4599  else
4600    rpaw=psxml%rpaw
4601  end if
4602 !shape_type, rshp:
4603  select case(trim(psxml%shape_function%gtype))
4604    case('gauss')
4605      shape_type=1
4606      rshp=rpaw
4607    case('bessel')
4608      shape_type=3
4609      rshp=psxml%shape_function%rc
4610    case('sinc')
4611      shape_type=2
4612      rshp=psxml%shape_function%rc
4613    case('exp')
4614      shape_type=1
4615      rshp=rpaw
4616    case('num')
4617      shape_type=-1
4618      rshp=rpaw
4619  end select
4620 
4621 end subroutine pawpsp_read_pawheader

m_pawpsp/pawpsp_rw_atompaw [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_rw_atompaw

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

1627 subroutine pawpsp_rw_atompaw(basis_size,filpsp,wvl)
1628 
1629 
1630 !This section has been created automatically by the script Abilint (TD).
1631 !Do not modify the following lines by hand.
1632 #undef ABI_FUNC
1633 #define ABI_FUNC 'pawpsp_rw_atompaw'
1634 !End of the abilint section
1635 
1636  implicit none
1637 
1638 !Arguments ------------------------------------
1639  integer,intent(in):: basis_size
1640  type(wvlpaw_type),intent(in)::wvl
1641  character(len=fnlen),intent(in)::filpsp
1642 !arrays
1643  character(strlen) :: pspline
1644  character(len=fnlen)::fname
1645 
1646 !Local variables-------------------------------
1647  integer :: ib,ii,ios,jj,step,iunt,ount
1648 !arrays
1649 
1650 ! *************************************************************************
1651  iunt = libpaw_get_free_unit()
1652  ount = libpaw_get_free_unit()
1653 
1654  step=0
1655 ! Open psp file for reading
1656   open(unit=iunt,file=trim(filpsp),form='formatted',status='old',action="read")
1657 ! Open the file for writing
1658   write(fname,'(2a)') libpaw_basename(trim(filpsp)),".wvl"
1659   open(unit=ount,file=fname,form='formatted',status='unknown',action="write")
1660 
1661   read_loop: do
1662     if(step==0) then
1663       read(iunt,'(a)',IOSTAT=ios) pspline
1664       if ( ios /= 0 ) exit read_loop
1665       if(index(trim(pspline),'CORE_DENSITY')/=0 .and. &
1666 &        index(trim(pspline),'PSEUDO_CORE_DENSITY')==0 .and. &
1667 &        index(trim(pspline),'TCORE_DENSITY')==0 ) then
1668         step=1
1669       else
1670         write(ount,'(a)') trim(pspline)
1671       end if
1672     elseif(step==1) then
1673 !Write Gaussian projectors:
1674       jj=0
1675       do ib=1,basis_size
1676         write(ount,'(a,i1,a)') "===== GAUSSIAN_TPROJECTOR ",ib,&
1677 &       " =====   "
1678         write(ount,'(i5,1x,i5,1x,a)')wvl%pngau(ib),wvl%ptotgau, ":ngauss, total ngauss"
1679         write(ount,'(3(1x,es23.16))')(wvl%parg(:,ii),&
1680 &        ii=jj+1,jj+wvl%pngau(ib))
1681         write(ount,'(3(1x,es23.16))')(wvl%pfac(:,ii),&
1682 &        ii=jj+1,jj+wvl%pngau(ib))
1683         jj=jj+wvl%pngau(ib)
1684       end do
1685       write(ount,'(a)') trim(pspline)
1686       step=0
1687     end if
1688   end do read_loop
1689 
1690   close(iunt)
1691   close(ount)
1692 
1693 end subroutine pawpsp_rw_atompaw

m_pawpsp/pawpsp_vhar2rho [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_vhar2rho

FUNCTION

  gets rho(r) from v(r), solving the Poisson equation
  \lap v(r) =  4 \pi rho(r)

INPUTS

 radmesh = radial grid (datastructure)
 vv(:)= potential

OUTPUT

  rho(:)= density

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

2714 subroutine pawpsp_vhar2rho(radmesh,rho,vv)
2715 
2716 
2717 !This section has been created automatically by the script Abilint (TD).
2718 !Do not modify the following lines by hand.
2719 #undef ABI_FUNC
2720 #define ABI_FUNC 'pawpsp_vhar2rho'
2721 !End of the abilint section
2722 
2723  implicit none
2724 
2725 !Arguments ------------------------------------
2726  type(pawrad_type),intent(in) :: radmesh
2727  real(dp), intent(in) :: vv(:)
2728  real(dp), intent(out):: rho(:)
2729 
2730 !Local variables-------------------------------
2731  integer :: nr
2732  real(dp) :: dfdr(radmesh%mesh_size),d2fdr(radmesh%mesh_size)
2733 
2734 ! *************************************************************************
2735 
2736  nr=size(vv)
2737  if (nr/=size(rho)) then
2738    MSG_BUG('wrong sizes!')
2739  end if
2740 
2741 !Laplacian =
2742 !\frac{\partial^2}{\partial r^2} + 2/r \frac{\partial}{\partial r}
2743 
2744 !Calculate derivatives
2745  call nderiv_gen(dfdr(1:nr),vv,radmesh,der2=d2fdr(1:nr))
2746 
2747  rho(2:nr)=d2fdr(2:nr) + 2._dp*dfdr(2:nr)/radmesh%rad(2:nr)
2748  call pawrad_deducer0(rho,nr,radmesh)
2749 
2750  rho(1:nr)=-rho(1:nr)/(4._dp*pi)
2751 
2752 end subroutine pawpsp_vhar2rho

m_pawpsp/pawpsp_wvl [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_wvl

FUNCTION

 WVL+PAW related operations

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp,pspatm

CHILDREN

SOURCE

4237 subroutine pawpsp_wvl(filpsp,pawrad, pawtab,usewvl, wvl_ngauss, comm_mpi)
4238 
4239 
4240 !This section has been created automatically by the script Abilint (TD).
4241 !Do not modify the following lines by hand.
4242 #undef ABI_FUNC
4243 #define ABI_FUNC 'pawpsp_wvl'
4244 !End of the abilint section
4245 
4246 implicit none
4247 
4248 !Arguments------------------------------------
4249 !scalars
4250  integer, optional,intent(in):: comm_mpi
4251  integer, intent(in):: usewvl, wvl_ngauss(2)
4252  character(len=fnlen),intent(in)::filpsp
4253  type(pawrad_type),intent(in) :: pawrad
4254  type(pawtab_type),intent(inout):: pawtab
4255 !arrays
4256 
4257 !Local variables-------------------------------
4258 !scalars
4259  integer:: ii, me, mparam, nterm_bounds(2)
4260  type(pawrad_type)::tproj_mesh
4261  character(len=500) :: msg
4262 !arrays
4263  integer,allocatable:: ngauss_param(:)
4264  real(dp),allocatable:: gauss_param(:,:)
4265 
4266 ! *************************************************************************
4267 
4268  me=0; if (present(comm_mpi))me=xmpi_comm_rank(comm_mpi)
4269 
4270 !If usewvl flag is on, we must have the pawtab%wvl pointer allocated
4271  if (usewvl==1.and.pawtab%has_wvl==0) then
4272    call wvlpaw_allocate(pawtab%wvl)
4273    pawtab%has_wvl=1
4274  end if
4275 
4276 !Fit projectors to a sum of Gaussians:
4277  if (usewvl ==1 .and. pawtab%wvl%ptotgau==0 ) then
4278 
4279    if (pawtab%has_tproj==0) then
4280      msg='pawtab%tproj must be allocated'
4281      MSG_BUG(msg)
4282    end if
4283 
4284 !  1) fit projectors to gaussians
4285    write(msg,'(a,a)')ch10,'Fitting tproj to Gaussians'
4286    call wrtout(std_out,msg,'COLL')
4287 
4288 !  See fit_gen (option==4):
4289    do ii=1,2
4290      nterm_bounds(ii)=ceiling(wvl_ngauss(ii)/2.0)
4291    end do
4292    mparam=nterm_bounds(2)*4
4293    LIBPAW_ALLOCATE(gauss_param,(mparam,pawtab%basis_size))
4294    LIBPAW_ALLOCATE(ngauss_param,(pawtab%basis_size))
4295 !  compute tproj_mesh
4296    call pawrad_init(tproj_mesh,mesh_size=size(pawtab%tproj,1),&
4297 &    mesh_type=pawrad%mesh_type,rstep=pawrad%rstep, lstep=pawrad%lstep)
4298 
4299    if(present(comm_mpi)) then
4300      call gaussfit_projector(pawtab%basis_size,mparam,&
4301 &     ngauss_param,nterm_bounds,pawtab%orbitals,&
4302 &     gauss_param,tproj_mesh,&
4303 &     pawtab%rpaw,pawtab%tproj,comm_mpi)
4304    else
4305      call gaussfit_projector(pawtab%basis_size,mparam,&
4306 &     ngauss_param,nterm_bounds,pawtab%orbitals,&
4307 &     gauss_param,tproj_mesh,&
4308 &     pawtab%rpaw,pawtab%tproj)
4309    end if
4310 !  tproj is now as a sum of sin+cos functions,
4311 !  convert it to a sum of complex gaussians and fill %wvl object:
4312    call pawpsp_wvl_sin2gauss(pawtab%basis_size,mparam,&
4313 &   ngauss_param,gauss_param,pawtab%wvl)
4314    LIBPAW_DEALLOCATE(gauss_param)
4315    LIBPAW_DEALLOCATE(ngauss_param)
4316 
4317    if(me==0) then
4318      call pawpsp_rw_atompaw(pawtab%basis_size,filpsp,pawtab%wvl)
4319    end if
4320 
4321    pawtab%has_wvl=2
4322 
4323  end if
4324 
4325 !Projectors in real space are no more needed
4326  call pawrad_free(tproj_mesh)
4327  if(allocated(pawtab%tproj)) then
4328    LIBPAW_DEALLOCATE(pawtab%tproj)
4329    pawtab%has_tproj=0
4330  end if
4331 
4332 end subroutine pawpsp_wvl

m_pawpsp/pawpsp_wvl_calc [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

 pawpsp_wvl_calc

FUNCTION

 Performs tests and compute data related to pspcod=7 or 17 ("PAW pseudopotentials")

INPUTS

  tnvale(vale_mesh%mesh_size)= pseudo valence density (+ nhat in output)
  usewvl= flag for wavelets method
  vale_mesh<type(pawrad_type)>= radial mesh for the valence density
  vloc_mesh<type(pawrad_type)>= radial mesh for the local potential
  vlocr(vloc_mesh%mesh_size)= local potential according to vlocopt.

OUTPUT

  Sets pawtab%rholoc

SIDE EFFECTS

  pawtab <type(pawtab_type)>= objects are modified

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

2786 subroutine pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr)
2787 
2788 
2789 !This section has been created automatically by the script Abilint (TD).
2790 !Do not modify the following lines by hand.
2791 #undef ABI_FUNC
2792 #define ABI_FUNC 'pawpsp_wvl_calc'
2793 !End of the abilint section
2794 
2795  implicit none
2796 
2797 !Arguments ------------------------------------
2798 !scalars
2799  integer,intent(in)::usewvl
2800  type(pawrad_type),intent(in) :: vale_mesh
2801  type(pawtab_type),intent(inout) :: pawtab
2802  type(pawrad_type),intent(in) ::vloc_mesh
2803 
2804 !arrays
2805  real(dp),intent(in) :: tnvale(vale_mesh%mesh_size*pawtab%has_tvale)
2806  real(dp),intent(in) :: vlocr(vloc_mesh%mesh_size)
2807 
2808 
2809 !Local variables ------------------------------
2810 !scalars
2811  integer :: msz
2812  character(len=500) :: msg
2813 !arrays
2814 
2815 ! *************************************************************************
2816 
2817 !If usewvl flag is on, we must have the pawtab%wvl pointer allocated
2818  if (pawtab%has_wvl==0) then
2819    msg='pawtab%has_wvl flag should be on o entry'
2820    MSG_BUG(msg)
2821  end if
2822  call wvlpaw_allocate(pawtab%wvl)
2823 
2824 !==========================================================
2825 !Change mesh_size of tvalespl
2826 !Compute second derivative from tNvale(r)
2827 
2828  if (pawtab%has_tvale/=0) then
2829    if(usewvl==1) then
2830      if(allocated(pawtab%tvalespl)) then
2831        LIBPAW_DEALLOCATE(pawtab%tvalespl)
2832      end if
2833      LIBPAW_ALLOCATE(pawtab%tvalespl,(vale_mesh%mesh_size,2))
2834      pawtab%tnvale_mesh_size=vale_mesh%mesh_size
2835      pawtab%tvalespl(:,1)=tnvale
2836 !    Compute second derivative of tvalespl(r)
2837      call paw_spline(vale_mesh%rad,pawtab%tvalespl(:,1),vale_mesh%mesh_size,zero,zero,pawtab%tvalespl(:,2))
2838    end if
2839  else
2840    pawtab%dnvdq0=zero
2841    pawtab%tnvale_mesh_size=0
2842  end if
2843 
2844 !==========================================================
2845 !Save rholoc:
2846 !Get local density from local potential
2847 !use the poisson eq.
2848  msz=vloc_mesh%mesh_size
2849  call wvlpaw_rholoc_free(pawtab%wvl%rholoc)
2850  LIBPAW_ALLOCATE(pawtab%wvl%rholoc%d,(msz,4))
2851  LIBPAW_ALLOCATE(pawtab%wvl%rholoc%rad,(msz))
2852  pawtab%wvl%rholoc%msz=msz
2853  pawtab%wvl%rholoc%rad(1:msz)=vloc_mesh%rad(1:msz)
2854 
2855 !get rho from v:
2856  call pawpsp_vhar2rho(vloc_mesh,pawtab%wvl%rholoc%d(:,1),vlocr)
2857 !
2858 !get second derivative, and store it
2859  call paw_spline(pawtab%wvl%rholoc%rad,pawtab%wvl%rholoc%d(:,1),pawtab%wvl%rholoc%msz,&
2860 & zero,zero,pawtab%wvl%rholoc%d(:,2))
2861 
2862 !save also vlocr:
2863  pawtab%wvl%rholoc%d(:,3)=vlocr
2864 
2865 !get second derivative, and store it
2866  call paw_spline(pawtab%wvl%rholoc%rad,vlocr,pawtab%wvl%rholoc%msz,&
2867 & zero,zero,pawtab%wvl%rholoc%d(:,4))
2868 
2869 !Test
2870 !do ii=1,pawtab%wvl%rholoc%msz
2871 !write(503,'(3(f16.10,x))')pawtab%wvl%rholoc%rad(ii),pawtab%wvl%rholoc%d(ii,1),pawtab%wvl%rholoc%d(ii,3)
2872 !end do
2873 !
2874 !Do splint
2875 !
2876 !nmesh=4000
2877 !rread1= (9.9979999d0/real(nmesh-1,dp)) ! 0.0025001d0  !step
2878 !allocate(raux1(nmesh),raux2(nmesh))
2879 !do ii=1,nmesh
2880 !raux1(ii)=rread1*real(ii-1,dp)  !mesh
2881 !end do
2882 !call splint(pawtab%wvl%rholoc%msz,pawtab%wvl%rholoc%rad,pawtab%wvl%rholoc%d(:,1),pawtab%wvl%rholoc%d(:,2),&
2883 !&  nmesh,raux1,raux2,ierr)
2884 !do ii=1,nmesh
2885 !write(401,'(10(f20.7,x))')raux1(ii),raux2(ii),raux2(ii)*raux1(ii)**2
2886 !end do
2887 !deallocate(raux1,raux2)
2888 
2889 end subroutine pawpsp_wvl_calc

m_pawpsp/pawpsp_wvl_sin2gauss [ Functions ]

[ Top ] [ m_pawpsp ] [ Functions ]

NAME

  pawpsp_wvl_sin2gauss

FUNCTION

  Converts a f(x)=sum_i^N_i a_i sin(b_i x)+ c_i cos( d_i x) to
    f(x)=sum_j e_j exp(f_j x), where e and f are complex numbers.

INPUTS

 basis_size =  size of the lmn basis
 mparam = number of terms in the summatory (N_i, see the expression above)
 nparam = Array containing the parameters (a_i, b_i,c_i,d_i)
 wvl = wavelets data type

OUTPUT

SIDE EFFECTS

 On output wvl%pfac and wvl%parg are filled with complex parameters (e_i, f_i)

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

3928  subroutine pawpsp_wvl_sin2gauss(basis_size,mparam,nparam,&
3929 & param,wvl)
3930 
3931 
3932 !This section has been created automatically by the script Abilint (TD).
3933 !Do not modify the following lines by hand.
3934 #undef ABI_FUNC
3935 #define ABI_FUNC 'pawpsp_wvl_sin2gauss'
3936 !End of the abilint section
3937 
3938   implicit none
3939   !
3940 !Arguments ------------------------------------
3941   integer,intent(in) :: mparam,basis_size
3942   integer,intent(in) :: nparam(basis_size)
3943   real(dp),intent(in) :: param(mparam,basis_size)
3944   type(wvlpaw_type),intent(inout):: wvl
3945 
3946 !Local variables ------------------------------
3947   integer :: i,ii,ib,ngauss,nterm
3948   real(dp) :: sep
3949   real(dp) :: a1(mparam),a2(mparam),a3(mparam),a4(mparam),a5(mparam)
3950   real(dp) :: b1r(mparam),b2r(mparam),b1i(mparam),b2i(mparam)
3951  character(len=500) :: message
3952   !
3953   !extra variables, use to debug
3954   !
3955   !integer::igau,nr,unitp
3956   !real(dp)::step,rmax
3957   !real(dp),allocatable::r(:), y(:)
3958   !complex::fac,arg
3959   !complex(dp),allocatable::f(:)
3960 !************************************************************************
3961 
3962 !  Convert from \sum(sin+cos) expressions to sums of complex gaussians
3963 !  (only works for option=4, see fit_gen)
3964 
3965 !  get number of coefficients:
3966    ii=0
3967    do ib=1,basis_size
3968      nterm=nparam(ib)/4  !option=4, there are 4 parameters for each term
3969      ii=ii+nterm*2 !two gaussians for each term
3970    end do
3971 !
3972 !  Allocate objects
3973 !
3974    ngauss=ii
3975    wvl%ptotgau=ngauss !total number of complex gaussians
3976    LIBPAW_ALLOCATE(wvl%pfac,(2,ngauss))
3977    LIBPAW_ALLOCATE(wvl%parg,(2,ngauss))
3978    LIBPAW_ALLOCATE(wvl%pngau,(basis_size))
3979    wvl%pngau(1:basis_size)=nparam(1:basis_size)/2 !option=4
3980 !
3981    ii=0
3982    do ib=1,basis_size
3983 !
3984 !    Get parameters in sin+cos expansion:
3985 !    Option4: \sum a1 exp(-a2 x^2) ( a3 sin(k x^2) + a4 cos(k x^2))
3986 !
3987      nterm=nparam(ib)/4  !option=4
3988 !
3989      a1(1:nterm)=param(1:nterm,ib)
3990      a2(1:nterm)=param(nterm+1:nterm*2,ib)
3991      a3(1:nterm)=param(nterm*2+1:nterm*3,ib)
3992      a4(1:nterm)=param(nterm*3+1:nterm*4,ib)
3993      sep=1.1d0
3994      do i=1,nterm
3995        a5(i)=sep**(i)
3996      end do
3997 
3998 !    First check that "a2" is a positive number (it is multiplied by -1, so
3999 !    that gaussians decay to zero:
4000      if( any(a2(1:nterm) < tol12) ) then
4001        message = 'Real part of Gaussians should be a negative number (they should go to zero at infty)'
4002        MSG_ERROR(message)
4003      end if
4004 
4005 !
4006 !    Now translate them to a sum of complex gaussians:
4007 !    pngau(ib)=nterm*2
4008 !    Two gaussians by term:
4009 !
4010 !    First gaussian
4011      b1r(1:nterm)= a1(1:nterm)*a4(1:nterm)/2.d0 !coefficient, real
4012      b1i(1:nterm)=-a1(1:nterm)*a3(1:nterm)/2.d0 !coefficient, imag
4013      b2r(1:nterm)=-a2(1:nterm)   !exponential, real
4014      b2i(1:nterm)= a5(1:nterm)   !exponential, imag
4015 !
4016      wvl%pfac(1,ii+1:ii+nterm)=b1r(1:nterm)
4017      wvl%pfac(2,ii+1:ii+nterm)=b1i(1:nterm)
4018      wvl%parg(1,ii+1:ii+nterm)=b2r(1:nterm)
4019      wvl%parg(2,ii+1:ii+nterm)=b2i(1:nterm)
4020 !    Second gaussian
4021      wvl%pfac(1,ii+nterm+1:ii+nterm*2)= b1r(1:nterm)
4022      wvl%pfac(2,ii+nterm+1:ii+nterm*2)=-b1i(1:nterm)
4023      wvl%parg(1,ii+nterm+1:ii+nterm*2)= b2r(1:nterm)
4024      wvl%parg(2,ii+nterm+1:ii+nterm*2)=-b2i(1:nterm)
4025 !
4026      ii=ii+nterm*2
4027    end do
4028 
4029 !  begin debug
4030 !  write(*,*)'pawpsp_wvl_sin2gauss, comment me'
4031 !  nr=3000
4032 !  rmax=10.d0
4033 !  LIBPAW_ALLOCATE(r,(nr))
4034 !  LIBPAW_ALLOCATE(f,(nr))
4035 !  LIBPAW_ALLOCATE(y,(nr))
4036 !  step=rmax/real(nr-1,dp)
4037 !  do ir=1,nr
4038 !  r(ir)=real(ir-1,dp)*step
4039 !  end do
4040 !  !
4041 !  ii=0
4042 !  do ib=1,basis_size
4043 !  unitp=500+ib
4044 !  f(:)=czero
4045 !  !
4046 !  do igau=1,wvl%pngau(ib)
4047 !  ii=ii+1
4048 !  arg=cmplx(wvl%parg(1,ii),wvl%parg(2,ii))
4049 !  fac=cmplx(wvl%pfac(1,ii),wvl%pfac(2,ii))
4050 !  f(:)=f(:)+fac*exp(arg*r(:)**2)
4051 !  end do
4052 !  do ir=1,nr
4053 !  write(unitp,'(3f16.7)')r(ir),real(f(ir))!,y(ir)
4054 !  end do
4055 !  end do
4056 !  LIBPAW_DEALLOCATE(r)
4057 !  LIBPAW_DEALLOCATE(f)
4058 !  LIBPAW_DEALLOCATE(y)
4059 !  end debug
4060 
4061  end subroutine pawpsp_wvl_sin2gauss

pawpsp_main/pawpsp_check_xml_upf [ Functions ]

[ Top ] [ pawpsp_main ] [ Functions ]

NAME

  pawpsp_main_checks

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

4963 subroutine pawpsp_check_xml_upf(filpsp)
4964 
4965 
4966 !This section has been created automatically by the script Abilint (TD).
4967 !Do not modify the following lines by hand.
4968 #undef ABI_FUNC
4969 #define ABI_FUNC 'pawpsp_check_xml_upf'
4970 !End of the abilint section
4971 
4972 implicit none
4973 
4974 !Arguments ------------------------------------
4975 !scalars
4976  character(len=fnlen),intent(in):: filpsp   ! name of the psp file
4977 
4978 !Local variables-------------------------------
4979  integer :: unt
4980  character(len=70):: testxml
4981 
4982 ! *************************************************************************
4983 
4984 !  Check if the file pseudopotential file is written in XML
4985    usexml = 0
4986    unt = libpaw_get_free_unit()
4987    open (unit=unt,file=filpsp,form='formatted',status='old',action="read")
4988    rewind (unit=unt)
4989    read(unt,*) testxml
4990    if(testxml(1:5)=='<?xml')then
4991      usexml = 1
4992      read(unt,*) testxml
4993      if(testxml(1:4)/='<paw')then
4994        msg='Reading a NC pseudopotential for a PAW calculation?'
4995        MSG_BUG(msg)
4996      end if
4997    else
4998      usexml = 0
4999    end if
5000    close (unit=unt)
5001 
5002 !  Check if pseudopotential file is a Q-espresso UPF file
5003    unt = libpaw_get_free_unit()
5004    open (unit=unt,file=filpsp,form='formatted',status='old',action="read")
5005    rewind (unit=unt)
5006    read(unt,*) testxml ! just a string, no relation to xml.
5007    if(testxml(1:9)=='<PP_INFO>')then
5008      msg='UPF format not allowed with PAW (USPP part not read yet)!'
5009      MSG_ERROR(msg)
5010    end if
5011    close (unit=unt)
5012 
5013 end subroutine pawpsp_check_xml_upf

pawpsp_main/pawpsp_consistency [ Functions ]

[ Top ] [ pawpsp_main ] [ Functions ]

NAME

  pawpsp_consistency

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

      m_pawpsp

CHILDREN

SOURCE

5040 subroutine pawpsp_consistency()
5041 
5042 
5043 !This section has been created automatically by the script Abilint (TD).
5044 !Do not modify the following lines by hand.
5045 #undef ABI_FUNC
5046 #define ABI_FUNC 'pawpsp_consistency'
5047 !End of the abilint section
5048 
5049 implicit none
5050 
5051 !Local variables-------------------------------
5052 
5053 ! *************************************************************************
5054 
5055 !Check pspcod=7 or 17
5056  if(pspcod/=7 .and. pspcod/=17)then
5057    write(msg, '(a,i2,a,a)' )&
5058 &   'In reading atomic psp file, finds pspcod=',pspcod,ch10,&
5059 &   'This is not an allowed value within PAW.'
5060    MSG_BUG(msg)
5061  end if
5062 
5063 !Does nuclear charge znuclpsp agree with psp input znucl
5064  if (abs(znuclpsp-znucl)>tol8) then
5065    write(msg, '(a,f10.5,2a,f10.5,5a)' )&
5066 &   'Pseudopotential file znucl=',znucl,ch10,&
5067 &   'does not equal input znuclpsp=',znuclpsp,' better than 1e-08 .',ch10,&
5068 &   'znucl is read from the psp file in pspatm_abinit, while',ch10,&
5069 &   'znuclpsp is read in iofn2.'
5070    MSG_BUG(msg)
5071  end if
5072 
5073 !Does nuclear charge zionpsp agree with psp input zion
5074  if (abs(zionpsp-zion)>tol8) then
5075    write(msg, '(a,f10.5,2a,f10.5,5a)' )&
5076 &   'Pseudopotential file zion=',zion,ch10,&
5077 &   'does not equal input zionpsp=',zionpsp,' better than 1e-08 .',ch10,&
5078 &   'zion is read from the psp file in pawpsp_main, while',ch10,&
5079 &   'zionpsp is read in iofn2.'
5080    MSG_BUG(msg)
5081  end if
5082 
5083 !Check several choices for ixc against pspxc
5084 !ixc is from ABINIT code; pspxc is from atomic psp file
5085  if (ixc==0) then
5086    msg='Note that input ixc=0 => no xc is being used.'
5087    MSG_WARNING(msg)
5088  else if(ixc/=pspxc) then
5089    write(msg, '(a,i8,a,a,a,i8,a,a,a,a,a,a,a,a,a,a)' ) &
5090 &   'Pseudopotential file pspxc=',pspxc,',',ch10,&
5091 &   'not equal to input ixc=',ixc,'.',ch10,&
5092 &   'These parameters must agree to get the same xc ',ch10,&
5093 &   'in ABINIT code as in psp construction.',ch10,&
5094 &   'Action : check psp design or input file.',ch10,&
5095 &   'Assume experienced user. Execution will continue.',ch10
5096    MSG_WARNING(msg)
5097  end if
5098 
5099  if (lloc>lmax ) then
5100    write(msg, '(a,2i12,a,a,a,a)' )&
5101 &   'lloc,lmax=',lloc,lmax,ch10,&
5102 &   'chosen l of local psp exceeds range from input data.',ch10,&
5103 &   'Action : check pseudopotential input file.'
5104    MSG_ERROR(msg)
5105  end if
5106 
5107 end subroutine pawpsp_consistency