TABLE OF CONTENTS
- ABINIT/m_pawpsp
- m_pawpsp/pawpsp_17in
- m_pawpsp/pawpsp_7in
- m_pawpsp/pawpsp_bcast
- m_pawpsp/pawpsp_calc
- m_pawpsp/pawpsp_calc_d5
- m_pawpsp/pawpsp_cg
- m_pawpsp/pawpsp_header_type
- m_pawpsp/pawpsp_lo
- m_pawpsp/pawpsp_main
- m_pawpsp/pawpsp_nl
- m_pawpsp/pawpsp_read
- m_pawpsp/pawpsp_read_corewf
- m_pawpsp/pawpsp_read_header
- m_pawpsp/pawpsp_read_header_2
- m_pawpsp/pawpsp_read_header_xml
- m_pawpsp/pawpsp_read_pawheader
- m_pawpsp/pawpsp_rw_atompaw
- m_pawpsp/pawpsp_vhar2rho
- m_pawpsp/pawpsp_wvl
- m_pawpsp/pawpsp_wvl_calc
- m_pawpsp/pawpsp_wvl_sin2gauss
- pawpsp_main/pawpsp_check_xml_upf
- pawpsp_main/pawpsp_consistency
ABINIT/m_pawpsp [ Modules ]
NAME
m_pawpsp
FUNCTION
Module to read PAW atomic data
COPYRIGHT
Copyright (C) 2012-2024 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
20 #include "libpaw.h" 21 22 module m_pawpsp 23 24 USE_DEFS 25 USE_MSG_HANDLING 26 USE_MPI_WRAPPERS 27 USE_MEMORY_PROFILING 28 29 use m_libpaw_libxc 30 #if defined LIBPAW_HAVE_FOX 31 use fox_sax 32 #endif 33 34 use m_libpaw_tools, only : libpaw_basename, libpaw_get_free_unit 35 36 use m_pawang, only: pawang_type 37 use m_pawtab, only: pawtab_type, wvlpaw_type, wvlpaw_allocate, wvlpaw_rholoc_free, & 38 & pawtab_free, wvlpaw_free, wvlpaw_rholoc_nullify, pawtab_bcast, & 39 & pawtab_set_flags, wvlpaw_allocate, wvlpaw_free, wvlpaw_rholoc_nullify, & 40 & wvlpaw_rholoc_free 41 use m_pawxmlps, only: rdpawpsxml_core, paw_setup_t, paw_setuploc, paw_setup_free 42 use m_pawrad, only: pawrad_type, pawrad_init, pawrad_free, pawrad_copy, & 43 & pawrad_bcast, pawrad_ifromr, simp_gen, nderiv_gen, bound_deriv, pawrad_deducer0, poisson 44 use m_paw_numeric, only: paw_splint, paw_spline, paw_smooth, paw_jbessel_4spline 45 use m_paw_atom, only: atompaw_shapebes, atompaw_vhnzc, atompaw_shpfun, & 46 & atompaw_dij0, atompaw_kij 47 use m_pawxc, only: pawxc, pawxcm, pawxc_get_usekden 48 use m_paw_gaussfit, only: gaussfit_projector 49 50 implicit none 51 52 private 53 54 public:: pawpsp_calc_d5 !calculate up to the 5th derivative 55 public:: pawpsp_main !main routine to read psp 56 public:: pawpsp_nl !make paw projector form factors f_l(q) 57 public:: pawpsp_read !read psp from file 58 public:: pawpsp_read_header !read header of psp file 59 public:: pawpsp_read_corewf !read core wavefunction 60 public:: pawpsp_read_header_2 !reads pspversion, basis_size and lmn_size 61 public:: pawpsp_rw_atompaw !read and writes ATOMPAW psp with gaussian |p> 62 public:: pawpsp_wvl !wavelet and icoulomb>0 related operations 63 public:: pawpsp_wvl_calc !wavelet related operations 64 public:: pawpsp_7in !reads non-XML atomic data 65 public:: pawpsp_17in !reads XML atomic data 66 public:: pawpsp_calc !calculates atomic quantities from psp info 67 public:: pawpsp_read_header_xml !read header of psp file for XML 68 public:: pawpsp_read_pawheader !read header variables from XML objects 69 public:: pawpsp_bcast ! broadcast PAW psp data 70 public:: pawpsp_cg !compute sin FFT transform of a density 71 public:: pawpsp_lo !compute sin FFT transform of local potential 72 73 ! Private procedures 74 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 xc_denpos= lowest allowed density (usually for the computation of the XC functionals) [xc_taupos]= lowest allowed kinetic energy density (for mGGA XC functionals) 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
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
SOURCE
3000 subroutine pawpsp_17in(epsatm,ffspl,icoulomb,ipsp,hyb_mixing,ixc,lmax,& 3001 & lnmax,mmax,mqgrid_ff,mqgrid_vl,pawpsp_header,pawrad,pawtab,& 3002 & pawxcdev, qgrid_ff,qgrid_vl,usewvl,usexcnhat_in,vlspl,xcccrc,& 3003 & xclevel,xc_denpos,zion,znucl,& 3004 & xc_taupos) ! Optional argument 3005 3006 !Arguments ------------------------------------ 3007 !scalars 3008 integer,intent(in) :: ipsp,ixc,lmax,lnmax,mqgrid_ff,mqgrid_vl,pawxcdev,usexcnhat_in 3009 integer,intent(inout) ::mmax 3010 integer,intent(in) :: xclevel,icoulomb,usewvl 3011 real(dp),intent(in) :: hyb_mixing,xc_denpos,zion,znucl 3012 real(dp),intent(in),optional :: xc_taupos 3013 real(dp),intent(out) :: epsatm,xcccrc 3014 type(pawpsp_header_type),intent(in) :: pawpsp_header 3015 type(pawrad_type),intent(inout) :: pawrad 3016 type(pawtab_type),intent(inout) :: pawtab 3017 !arrays 3018 real(dp),intent(in) :: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl) 3019 real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax) 3020 real(dp),intent(out) :: vlspl(mqgrid_vl,2) 3021 3022 !Local variables ------------------------------ 3023 !scalars 3024 integer :: has_v_minushalf,ib,icoremesh,icoretaumesh,il,ilm,ilmn,ilmn0,iln,imainmesh,imsh,iprojmesh 3025 integer :: ir,iread1,ishpfmesh,ivalemesh,ivlocmesh,j0lmn,jlm,pngau 3026 integer :: jlmn,jln,klmn,msz,nmesh,nval,pspversion,shft,sz10,usexcnhat,vlocopt 3027 real(dp), parameter :: rmax_vloc=10.0_dp 3028 real(dp) :: fourpi,my_xc_taupos,occ,rc,yp1,ypn 3029 logical :: save_core_msz 3030 character(len=500) :: msg 3031 type(pawrad_type) :: core_mesh,coretau_mesh,shpf_mesh,tproj_mesh,vale_mesh,vloc_mesh 3032 !arrays 3033 integer,allocatable :: mesh_shift(:),nprj(:) 3034 real(dp),allocatable :: kij(:),ncore(:),shpf(:,:),tncore(:),coretau(:) 3035 real(dp),allocatable :: tcoretau(:),tnvale(:),tproj(:,:),vhnzc(:),vlocr(:) 3036 real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:) 3037 type(pawrad_type),allocatable :: radmesh(:) 3038 3039 !************************************************************************ 3040 3041 if (.False.) write(std_out,*) ipsp 3042 3043 !========================================================== 3044 !Destroy everything in pawtab but optional flags 3045 call pawtab_free(pawtab) 3046 !Destroy everything in pawrad 3047 call pawrad_free(pawrad) 3048 3049 !========================================================== 3050 !Initialize useful data 3051 3052 pawtab%usexcnhat=usexcnhat_in 3053 fourpi=4*acos(-1.d0) 3054 pspversion=pawpsp_header%pawver 3055 save_core_msz=(usewvl==1 .or. icoulomb .ne. 0) 3056 imainmesh=-1;icoremesh=-1;icoretaumesh=-1;iprojmesh=-1 3057 ishpfmesh=-1;ivalemesh=-1;ivlocmesh=-1 3058 my_xc_taupos=xc_denpos;if(present(xc_taupos)) my_xc_taupos=xc_taupos 3059 3060 !========================================================== 3061 !Initialize partial waves quantum numbers 3062 3063 pawtab%basis_size=pawpsp_header%basis_size 3064 LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size)) 3065 do ib=1,pawtab%basis_size 3066 pawtab%orbitals(ib)=paw_setuploc%valence_states%state(ib)%ll 3067 end do 3068 3069 !========================================================== 3070 !Initialize various dims and indexes 3071 3072 pawtab%lmn_size=pawpsp_header%lmn_size 3073 pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2 3074 pawtab%l_size=2*maxval(pawtab%orbitals)+1 3075 pawtab%ij_size=pawtab%basis_size*(pawtab%basis_size+1)/2 3076 3077 !indlmn calculation (indices for (l,m,n) basis) 3078 if (allocated(pawtab%indlmn)) then 3079 LIBPAW_DEALLOCATE(pawtab%indlmn) 3080 end if 3081 LIBPAW_ALLOCATE(pawtab%indlmn,(6,pawtab%lmn_size)) 3082 pawtab%indlmn(:,:)=0 3083 LIBPAW_BOUND1_ALLOCATE(nprj,BOUNDS(0,maxval(pawtab%orbitals))) 3084 ilmn=0;iln=0;nprj=0 3085 do ib=1,pawtab%basis_size 3086 il=pawtab%orbitals(ib) 3087 nprj(il)=nprj(il)+1 3088 iln=iln+1 3089 do ilm=1,2*il+1 3090 pawtab%indlmn(1,ilmn+ilm)=il 3091 pawtab%indlmn(2,ilmn+ilm)=ilm-(il+1) 3092 pawtab%indlmn(3,ilmn+ilm)=nprj(il) 3093 pawtab%indlmn(4,ilmn+ilm)=il*il+ilm 3094 pawtab%indlmn(5,ilmn+ilm)=iln 3095 pawtab%indlmn(6,ilmn+ilm)=1 3096 end do 3097 ilmn=ilmn+2*il+1 3098 end do 3099 LIBPAW_DEALLOCATE(nprj) 3100 !Are ilmn (found here) and pawtab%lmn_size compatibles ? 3101 if (ilmn/=pawtab%lmn_size) then 3102 write(msg, '(a,a,a,a,a)' )& 3103 & 'Calculated lmn size differs from',ch10,& 3104 & 'lmn_size read from pseudo !',ch10,& 3105 & 'Action: check your pseudopotential file.' 3106 LIBPAW_ERROR(msg) 3107 end if 3108 3109 !========================================================== 3110 !Read and initialize radial meshes 3111 3112 nmesh=paw_setuploc%ngrid 3113 LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh)) 3114 LIBPAW_ALLOCATE(mesh_shift,(nmesh)) 3115 do imsh=1,nmesh 3116 radmesh(imsh)%mesh_type=-1 3117 radmesh(imsh)%rstep=zero 3118 radmesh(imsh)%lstep=zero 3119 mesh_shift(imsh)=0 3120 select case(trim(paw_setuploc%radial_grid(imsh)%eq)) 3121 case("r=a*exp(d*i)") 3122 mesh_shift(imsh)=1 3123 radmesh(imsh)%mesh_type=3 3124 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3125 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3126 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa 3127 radmesh(imsh)%lstep=paw_setuploc%radial_grid(imsh)%dd 3128 case("r=a*i/(1-b*i)") 3129 write(msg, '(3a)' )& 3130 & 'The grid r=a*i/(1-b*i) is not implemented in ABINIT !',ch10,& 3131 & 'Action: check your psp file.' 3132 LIBPAW_ERROR(msg) 3133 case("r=a*i/(n-i)") 3134 mesh_shift(imsh)=0 3135 radmesh(imsh)%mesh_type=5 3136 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3137 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3138 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa 3139 radmesh(imsh)%lstep=dble(paw_setuploc%radial_grid(imsh)%nn) 3140 case("r=a*(exp(d*i)-1)") 3141 mesh_shift(imsh)=0 3142 radmesh(imsh)%mesh_type=2 3143 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3144 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3145 if(paw_setuploc%radial_grid(imsh)%istart==1)radmesh(imsh)%mesh_size=radmesh(imsh)%mesh_size+1 3146 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa 3147 radmesh(imsh)%lstep=paw_setuploc%radial_grid(imsh)%dd 3148 case("r=d*i") 3149 mesh_shift(imsh)=0 3150 radmesh(imsh)%mesh_type=1 3151 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3152 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3153 if(paw_setuploc%radial_grid(imsh)%istart==1)radmesh(imsh)%mesh_size=radmesh(imsh)%mesh_size+1 3154 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%dd 3155 case("r=(i/n+a)^5/a-a^4") 3156 write(msg, '(3a)' )& 3157 & 'The grid r=(i/n+a)^5/a-a^4 is not implemented in ABINIT !',ch10,& 3158 & 'Action: check your psp file.' 3159 LIBPAW_ERROR(msg) 3160 end select 3161 end do 3162 3163 !Initialize radial meshes 3164 do imsh=1,nmesh 3165 call pawrad_init(radmesh(imsh)) 3166 end do 3167 3168 pawtab%rpaw=pawpsp_header%rpaw 3169 3170 !========================================================== 3171 !Here reading shapefunction parameters 3172 3173 pawtab%shape_type=pawpsp_header%shape_type 3174 pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99 3175 pawtab%rshp=pawpsp_header%rshp 3176 pawtab%shape_lambda=paw_setuploc%shape_function%lamb 3177 if(trim(paw_setuploc%shape_function%gtype)=="gauss")pawtab%shape_lambda=2 3178 pawtab%shape_sigma=paw_setuploc%shape_function%rc 3179 !If shapefunction type is gaussian, check exponent 3180 if (pawtab%shape_type==1) then 3181 if (pawtab%shape_lambda<2) then 3182 write(msg, '(3a)' )& 3183 & 'For a gaussian shape function, exponent lambda must be >1 !',ch10,& 3184 & 'Action: check your psp file.' 3185 LIBPAW_ERROR(msg) 3186 end if 3187 end if 3188 3189 !If shapefunction type is Bessel, deduce here its parameters from rc 3190 if (pawtab%shape_type==3) then 3191 LIBPAW_ALLOCATE(pawtab%shape_alpha,(2,pawtab%l_size)) 3192 LIBPAW_ALLOCATE(pawtab%shape_q,(2,pawtab%l_size)) 3193 rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw 3194 do il=1,pawtab%l_size 3195 call atompaw_shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc) 3196 end do 3197 end if 3198 3199 !========================================================== 3200 !Mirror pseudopotential parameters to the output and log files 3201 3202 write(msg,'(a,i2)')' Pseudopotential format is: paw',pspversion 3203 call wrtout(ab_out,msg,'COLL') 3204 call wrtout(std_out, msg,'COLL') 3205 write(msg,'(2(a,i3),a,64i4)') & 3206 & ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',& 3207 & pawtab%lmn_size,'), orbitals=',pawtab%orbitals(1:pawtab%basis_size) 3208 call wrtout(ab_out,msg,'COLL') 3209 call wrtout(std_out, msg,'COLL') 3210 write(msg,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw 3211 call wrtout(ab_out,msg,'COLL') 3212 call wrtout(std_out, msg,'COLL') 3213 write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:' 3214 call wrtout(ab_out,msg,'COLL') 3215 call wrtout(std_out, msg,'COLL') 3216 3217 do imsh=1,nmesh 3218 if (radmesh(imsh)%mesh_type==1) & 3219 & write(msg,'(a,i1,a,i4,a,g12.5)') & 3220 & ' - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,& 3221 & ' , step=',radmesh(imsh)%rstep 3222 if (radmesh(imsh)%mesh_type==2) & 3223 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 3224 & ' - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,& 3225 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 3226 if (radmesh(imsh)%mesh_type==3) & 3227 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 3228 & ' - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,& 3229 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 3230 if (radmesh(imsh)%mesh_type==4) & 3231 & write(msg,'(a,i1,a,i4,a,g12.5)') & 3232 & ' - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,& 3233 & ' , AA=',radmesh(imsh)%rstep 3234 if (radmesh(imsh)%mesh_type==5) & 3235 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 3236 & ' - mesh ',imsh,': r(i)=-AA*i/(NN-i)), n=size=',radmesh(imsh)%mesh_size,& 3237 & ' , AA=',radmesh(imsh)%rstep,' NN=',radmesh(imsh)%lstep 3238 call wrtout(ab_out,msg,'COLL') 3239 call wrtout(std_out, msg,'COLL') 3240 end do 3241 if (pawtab%shape_type==-1) then 3242 write(msg,'(a)')& 3243 ' Shapefunction is NUMERIC type: directly read from atomic data file' 3244 call wrtout(ab_out,msg,'COLL') 3245 call wrtout(std_out, msg,'COLL') 3246 end if 3247 if (pawtab%shape_type==1) then 3248 write(msg,'(2a,a,f6.3,a,i3)')& 3249 & ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,& 3250 & ' with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda 3251 call wrtout(ab_out,msg,'COLL') 3252 call wrtout(std_out, msg,'COLL') 3253 end if 3254 if (pawtab%shape_type==2) then 3255 write(msg,'(a)')& 3256 ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2' 3257 call wrtout(ab_out,msg,'COLL') 3258 call wrtout(std_out, msg,'COLL') 3259 end if 3260 if (pawtab%shape_type==3) then 3261 write(msg,'(a)')& 3262 & ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)' 3263 call wrtout(ab_out,msg,'COLL') 3264 call wrtout(std_out, msg,'COLL') 3265 end if 3266 if (pawtab%rshp<1.d-8) then 3267 write(msg,'(a)') ' Radius for shape functions = sphere core radius' 3268 else 3269 write(msg,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp 3270 end if 3271 call wrtout(ab_out,msg,'COLL') 3272 call wrtout(std_out, msg,'COLL') 3273 3274 !========================================================== 3275 !Perfom tests 3276 3277 !Are lmax and orbitals compatibles ? 3278 if (lmax/=maxval(pawtab%orbitals)) then 3279 write(msg, '(a,a,a)' )& 3280 & 'lmax /= MAX(orbitals) !',ch10,& 3281 & 'Action: check your pseudopotential file.' 3282 LIBPAW_ERROR(msg) 3283 end if 3284 3285 !Only mesh_type=1,2, 3 or 5 allowed 3286 do imsh=1,nmesh 3287 if (radmesh(imsh)%mesh_type>5) then 3288 write(msg, '(a,a,a)' )& 3289 & 'Only mesh types 1,2,3 or 5 allowed !',ch10,& 3290 & 'Action : check your pseudopotential or input file.' 3291 LIBPAW_ERROR(msg) 3292 end if 3293 end do 3294 3295 !========================================================== 3296 !Read tabulated atomic data 3297 3298 !--------------------------------- 3299 !Read wave-functions (phi) 3300 3301 do ib=1,pawtab%basis_size 3302 if (ib==1) then 3303 do imsh=1,nmesh 3304 if(trim(paw_setuploc%ae_partial_wave(1)%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3305 mmax=radmesh(imsh)%mesh_size 3306 call pawrad_init(pawrad,mesh_size=mmax,mesh_type=radmesh(imsh)%mesh_type, & 3307 & rstep=radmesh(imsh)%rstep,lstep=radmesh(imsh)%lstep,r_for_intg=pawtab%rpaw) 3308 pawtab%partialwave_mesh_size=pawrad%mesh_size 3309 pawtab%mesh_size=pawrad_ifromr(pawrad,pawtab%rpaw)+5 3310 pawtab%mesh_size=min(pawtab%mesh_size,pawrad%mesh_size) 3311 if (pawtab%mesh_size>pawrad%mesh_size-2) pawtab%mesh_size=pawrad%mesh_size 3312 imainmesh=imsh 3313 exit 3314 end if 3315 end do 3316 LIBPAW_ALLOCATE(pawtab%phi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 3317 else if (trim(paw_setuploc%ae_partial_wave(ib)%grid)/=trim(paw_setuploc%radial_grid(imainmesh)%id)) then 3318 write(msg, '(a,a,a)' )& 3319 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 3320 & 'Action: check your pseudopotential file.' 3321 LIBPAW_ERROR(msg) 3322 end if 3323 shft=mesh_shift(imainmesh) 3324 pawtab%phi(1+shft:pawtab%partialwave_mesh_size,ib)= & 3325 & paw_setuploc%ae_partial_wave(ib)%data(1:pawtab%partialwave_mesh_size-shft) & 3326 & *pawrad%rad(1+shft:pawtab%partialwave_mesh_size) 3327 if (shft==1) pawtab%phi(1,ib)=zero 3328 end do 3329 write(msg,'(a,i4)') ' mmax= ',mmax 3330 call wrtout(ab_out,msg,'COLL') 3331 call wrtout(std_out,msg,'COLL') 3332 3333 !--------------------------------- 3334 !Read pseudo wave-functions (tphi) 3335 3336 LIBPAW_ALLOCATE(pawtab%tphi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 3337 do ib=1,pawtab%basis_size 3338 3339 if(trim(paw_setuploc%pseudo_partial_wave(ib)%grid)/=trim(paw_setuploc%radial_grid(imainmesh)%id)) then 3340 write(msg, '(a,a,a)' )& 3341 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 3342 & 'Action: check your pseudopotential file.' 3343 LIBPAW_ERROR(msg) 3344 end if 3345 shft=mesh_shift(imainmesh) 3346 pawtab%tphi(1+shft:pawtab%partialwave_mesh_size,ib)=& 3347 & paw_setuploc%pseudo_partial_wave(ib)%data(1:pawtab%partialwave_mesh_size-shft) & 3348 & *pawrad%rad(1+shft:pawtab%partialwave_mesh_size) 3349 if (shft==1) pawtab%tphi(1,ib)=zero 3350 end do 3351 write(msg,'(a,i1)') & 3352 & ' Radial grid used for partial waves is grid ',imainmesh 3353 call wrtout(ab_out,msg,'COLL') 3354 call wrtout(std_out, msg,'COLL') 3355 3356 !--------------------------------- 3357 !Read projectors (tproj) 3358 3359 if (allocated(paw_setuploc%projector_fit)) then 3360 call wvlpaw_allocate(pawtab%wvl) 3361 LIBPAW_ALLOCATE(pawtab%wvl%pngau,(pawtab%basis_size)) 3362 do ib=1,pawtab%basis_size 3363 pawtab%wvl%pngau(ib) = paw_setuploc%projector_fit(ib)%ngauss 3364 end do 3365 pawtab%wvl%ptotgau = sum(pawtab%wvl%pngau) * 2 3366 LIBPAW_ALLOCATE(pawtab%wvl%parg,(2,pawtab%wvl%ptotgau)) 3367 LIBPAW_ALLOCATE(pawtab%wvl%pfac,(2,pawtab%wvl%ptotgau)) 3368 pngau = 1 3369 do ib=1,pawtab%basis_size 3370 ! Complex gaussian 3371 pawtab%wvl%parg(:,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3372 & paw_setuploc%projector_fit(ib)%expos(:,1:pawtab%wvl%pngau(ib)) 3373 pawtab%wvl%pfac(:,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3374 & paw_setuploc%projector_fit(ib)%factors(:,1:pawtab%wvl%pngau(ib)) 3375 pngau = pngau + pawtab%wvl%pngau(ib) 3376 ! Conjugate gaussian 3377 pawtab%wvl%parg(1,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3378 & paw_setuploc%projector_fit(ib)%expos(1,1:pawtab%wvl%pngau(ib)) 3379 pawtab%wvl%parg(2,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3380 & -paw_setuploc%projector_fit(ib)%expos(2,1:pawtab%wvl%pngau(ib)) 3381 pawtab%wvl%pfac(1,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3382 & paw_setuploc%projector_fit(ib)%factors(1,1:pawtab%wvl%pngau(ib)) 3383 pawtab%wvl%pfac(2,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3384 & -paw_setuploc%projector_fit(ib)%factors(2,1:pawtab%wvl%pngau(ib)) 3385 pngau = pngau + pawtab%wvl%pngau(ib) 3386 pawtab%wvl%pngau(ib) = pawtab%wvl%pngau(ib) * 2 3387 end do 3388 pawtab%has_wvl=2 3389 else 3390 !Nullify wavelet objects for safety: 3391 pawtab%has_wvl=0 3392 call wvlpaw_free(pawtab%wvl) 3393 end if 3394 do ib=1,pawtab%basis_size 3395 if (ib==1) then 3396 do imsh=1,nmesh 3397 if(trim(paw_setuploc%projector_function(1)%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3398 iprojmesh=imsh 3399 exit 3400 end if 3401 end do 3402 call pawrad_copy(radmesh(iprojmesh),tproj_mesh) 3403 LIBPAW_ALLOCATE(tproj,(tproj_mesh%mesh_size,pawtab%basis_size)) 3404 else if (trim(paw_setuploc%projector_function(ib)%grid)/=trim(paw_setuploc%radial_grid(iprojmesh)%id)) then 3405 write(msg, '(a,a,a)' )& 3406 & 'All tprojectors must be given on the same radial mesh !',ch10,& 3407 & 'Action: check your pseudopotential file.' 3408 LIBPAW_ERROR(msg) 3409 end if 3410 shft=mesh_shift(iprojmesh) 3411 tproj(1+shft:tproj_mesh%mesh_size,ib)=paw_setuploc%projector_function(ib)%data(1:tproj_mesh%mesh_size-shft)& 3412 & *tproj_mesh%rad(1+shft:tproj_mesh%mesh_size) 3413 if (shft==1) tproj(1,ib)=zero 3414 end do 3415 write(msg,'(a,i1)') & 3416 & ' Radial grid used for projectors is grid ',iprojmesh 3417 call wrtout(ab_out,msg,'COLL') 3418 call wrtout(std_out, msg,'COLL') 3419 3420 !--------------------------------- 3421 !Read core density (coredens) 3422 3423 do imsh=1,nmesh 3424 if(trim(paw_setuploc%ae_core_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3425 icoremesh=imsh 3426 exit 3427 end if 3428 end do 3429 call pawrad_copy(radmesh(icoremesh),core_mesh) 3430 if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.& 3431 & (radmesh(icoremesh)%rstep /=pawrad%rstep) .or.& 3432 & (radmesh(icoremesh)%lstep /=pawrad%lstep)) then 3433 write(msg, '(a,a,a,a,a)' )& 3434 & 'Ncore must be given on a radial mesh with the same',ch10,& 3435 & 'type and step(s) than the main radial mesh (mesh for Phi) !',ch10,& 3436 & 'Action: check your pseudopotential file.' 3437 LIBPAW_ERROR(msg) 3438 end if 3439 LIBPAW_ALLOCATE(ncore,(core_mesh%mesh_size)) 3440 shft=mesh_shift(icoremesh) 3441 ncore(1+shft:core_mesh%mesh_size)=paw_setuploc%ae_core_density%data(1:core_mesh%mesh_size-shft)/sqrt(fourpi) 3442 if (shft==1) call pawrad_deducer0(ncore,core_mesh%mesh_size,core_mesh) 3443 3444 !Construct and save VH[z_NC] if requested 3445 if (pawtab%has_vhnzc==1) then 3446 LIBPAW_ALLOCATE(pawtab%VHnZC,(pawtab%mesh_size)) 3447 LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size)) 3448 call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl) 3449 pawtab%VHnZC(1:pawtab%mesh_size)=vhnzc(1:pawtab%mesh_size) 3450 pawtab%has_vhnzc=2 3451 LIBPAW_DEALLOCATE(vhnzc) 3452 end if 3453 3454 pawtab%core_mesh_size=pawtab%mesh_size 3455 if(save_core_msz) pawtab%core_mesh_size=core_mesh%mesh_size 3456 LIBPAW_ALLOCATE(pawtab%coredens,(pawtab%core_mesh_size)) 3457 pawtab%rcore=core_mesh%rad(pawtab%core_mesh_size) 3458 pawtab%coredens(1:pawtab%core_mesh_size)=ncore(1:pawtab%core_mesh_size) 3459 3460 !--------------------------------- 3461 !Read pseudo core density (tcoredens) 3462 3463 do imsh=1,nmesh 3464 if(trim(paw_setuploc%pseudo_core_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3465 iread1=imsh 3466 exit 3467 end if 3468 end do 3469 if (iread1/=icoremesh) then 3470 write(msg, '(a,a,a,a,a,a,a,a)' )& 3471 & 'Pseudized core density (tNcore) must be given',ch10,& 3472 & 'on the same radial mesh as core density (Ncore) !',ch10,& 3473 & 'Action: check your pseudopotential file.' 3474 LIBPAW_ERROR(msg) 3475 end if 3476 LIBPAW_ALLOCATE(tncore,(core_mesh%mesh_size)) 3477 shft=mesh_shift(icoremesh) 3478 tncore(1+shft:core_mesh%mesh_size)=paw_setuploc%pseudo_core_density%data(1:core_mesh%mesh_size-shft)/sqrt(fourpi) 3479 if (shft==1) call pawrad_deducer0(tncore,core_mesh%mesh_size,core_mesh) 3480 if(save_core_msz) then 3481 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,6)) 3482 else 3483 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,1)) 3484 end if 3485 if (maxval(abs(tncore(:)))<tol6) then 3486 pawtab%usetcore=0 3487 pawtab%tcoredens(1:pawtab%core_mesh_size,:)=zero 3488 else 3489 pawtab%usetcore=1 3490 pawtab%tcoredens(1:pawtab%core_mesh_size,1)=tncore(1:pawtab%core_mesh_size) 3491 end if 3492 write(msg,'(a,i1)') & 3493 & ' Radial grid used for (t)core density is grid ',icoremesh 3494 call wrtout(ab_out,msg,'COLL') 3495 call wrtout(std_out, msg,'COLL') 3496 3497 !--------------------------------- 3498 !Read core kinetic density (coretau) 3499 3500 if (paw_setuploc%ae_core_kinetic_energy_density%tread.and.pawtab%has_coretau>=1) then 3501 do imsh=1,nmesh 3502 if(trim(paw_setuploc%ae_core_kinetic_energy_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3503 icoretaumesh=imsh 3504 exit 3505 end if 3506 end do 3507 call pawrad_copy(radmesh(icoretaumesh),coretau_mesh) 3508 if (icoretaumesh/=icoremesh) then 3509 write(msg, '(5a)' )& 3510 & 'Core kinetic density (TAUcore) must be given',ch10,& 3511 & 'on the same radial mesh as core density (Ncore) !',ch10,& 3512 & 'Action: check your pseudopotential file.' 3513 LIBPAW_ERROR(msg) 3514 end if 3515 LIBPAW_ALLOCATE(coretau,(coretau_mesh%mesh_size)) 3516 shft=mesh_shift(icoretaumesh) 3517 coretau(1+shft:coretau_mesh%mesh_size)= & 3518 & paw_setuploc%ae_core_kinetic_energy_density%data(1:coretau_mesh%mesh_size-shft)/sqrt(fourpi) 3519 if (shft==1) call pawrad_deducer0(coretau,coretau_mesh%mesh_size,coretau_mesh) 3520 pawtab%coretau_mesh_size=pawtab%mesh_size 3521 if(save_core_msz) pawtab%coretau_mesh_size=coretau_mesh%mesh_size 3522 LIBPAW_ALLOCATE(pawtab%coretau,(pawtab%coretau_mesh_size)) 3523 pawtab%rcoretau=coretau_mesh%rad(pawtab%coretau_mesh_size) 3524 pawtab%coretau(1:pawtab%coretau_mesh_size)=coretau(1:pawtab%coretau_mesh_size) 3525 else if (pawtab%has_coretau>=1) then 3526 write(msg, '(5a)' )& 3527 & 'metaGGA exchange-correlation is requested but the core kinetic energy density',ch10,& 3528 & 'is not present in the pseudopotential file!',ch10,& 3529 & 'Action: check your pseudopotential file.' 3530 LIBPAW_ERROR(msg) 3531 end if 3532 3533 !--------------------------------- 3534 !Read pseudo core kinetic energy density (tcoretau) 3535 3536 if (paw_setuploc%pseudo_core_kinetic_energy_density%tread.and.pawtab%has_coretau>=1) then 3537 do imsh=1,nmesh 3538 if(trim(paw_setuploc%pseudo_core_kinetic_energy_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3539 iread1=imsh 3540 exit 3541 end if 3542 end do 3543 if (iread1/=icoretaumesh) then 3544 write(msg, '(5a)' )& 3545 & 'Pseudized core kinetic energy density (tTAUcore) must be given',ch10,& 3546 & 'on the same radial mesh as core kinetic density (TAUcore) !',ch10,& 3547 & 'Action: check your pseudopotential file.' 3548 LIBPAW_ERROR(msg) 3549 end if 3550 LIBPAW_ALLOCATE(tcoretau,(coretau_mesh%mesh_size)) 3551 shft=mesh_shift(icoretaumesh) 3552 tcoretau(1+shft:coretau_mesh%mesh_size)= & 3553 & paw_setuploc%pseudo_core_kinetic_energy_density%data(1:coretau_mesh%mesh_size-shft)/sqrt(fourpi) 3554 if (shft==1) call pawrad_deducer0(tcoretau,coretau_mesh%mesh_size,coretau_mesh) 3555 LIBPAW_ALLOCATE(pawtab%tcoretau,(pawtab%coretau_mesh_size)) 3556 pawtab%tcoretau(1:pawtab%coretau_mesh_size)=tcoretau(1:pawtab%coretau_mesh_size) 3557 pawtab%has_coretau=2 3558 write(msg,'(a,i1)') & 3559 & ' Radial grid used for (t)coretau kinetic density is grid ',icoretaumesh 3560 call wrtout(ab_out,msg,'COLL') 3561 call wrtout(std_out, msg,'COLL') 3562 else if (pawtab%has_coretau>=1) then 3563 write(msg, '(5a)' )& 3564 & 'metaGGA exchange-correlation is requested but the pseudo core kinetic energy density',ch10,& 3565 & 'is not present in the pseudopotential file!',ch10,& 3566 & 'Action: check your pseudopotential file.' 3567 LIBPAW_ERROR(msg) 3568 end if 3569 3570 !--------------------------------- 3571 !Read local pseudopotential=Vh(tn_zc) or Vbare 3572 3573 if ((paw_setuploc%blochl_local_ionic_potential%tread).and.& 3574 & (pawtab%usexcnhat==-1.or.pawtab%usexcnhat==0.or.(pawtab%usexcnhat==1.and.& 3575 & ((.not.paw_setuploc%zero_potential%tread).or.(.not.paw_setuploc%kresse_joubert_local_ionic_potential%tread))))) then 3576 usexcnhat=0;vlocopt=2 3577 do imsh=1,nmesh 3578 if(trim(paw_setuploc%blochl_local_ionic_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3579 iread1=imsh 3580 exit 3581 end if 3582 end do 3583 ivlocmesh=iread1 3584 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 3585 LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 3586 shft=mesh_shift(ivlocmesh) 3587 vlocr(1+shft:vloc_mesh%mesh_size)=paw_setuploc%blochl_local_ionic_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3588 if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh) 3589 else if((paw_setuploc%kresse_joubert_local_ionic_potential%tread).and.& 3590 & (pawtab%usexcnhat==-1.or.pawtab%usexcnhat==1.or.(pawtab%usexcnhat==0.and.& 3591 & (.not.paw_setuploc%zero_potential%tread)))) then 3592 usexcnhat=1;vlocopt=1 3593 do imsh=1,nmesh 3594 if(trim(paw_setuploc%kresse_joubert_local_ionic_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3595 iread1=imsh 3596 exit 3597 end if 3598 end do 3599 ivlocmesh=iread1 3600 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 3601 LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 3602 shft=mesh_shift(ivlocmesh) 3603 vlocr(1+shft:vloc_mesh%mesh_size)= & 3604 & paw_setuploc%kresse_joubert_local_ionic_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3605 if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh) 3606 else if(paw_setuploc%zero_potential%tread) then 3607 usexcnhat=0;vlocopt=0 3608 do imsh=1,nmesh 3609 if(trim(paw_setuploc%zero_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3610 iread1=imsh 3611 exit 3612 end if 3613 end do 3614 ivlocmesh=iread1 3615 ! vloc_mesh%mesh_type=radmesh(ivlocmesh)%mesh_type 3616 ! vloc_mesh%rstep=radmesh(ivlocmesh)%rstep 3617 ! vloc_mesh%lstep=radmesh(ivlocmesh)%lstep 3618 ! vloc_mesh%mesh_size=radmesh(ivlocmesh)%mesh_size 3619 ! vloc_mesh%mesh_size=pawrad_ifromr(radmesh(ivlocmesh),rmax_vloc) 3620 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 3621 LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 3622 vlocr=zero 3623 shft=mesh_shift(ivlocmesh) 3624 vlocr(1+shft:vloc_mesh%mesh_size)=paw_setuploc%zero_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3625 if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh) 3626 else 3627 write(msg, '(a,a,a,a,a)' )& 3628 & 'At least one local potential must be given',ch10,& 3629 & 'Action: check your pseudopotential file.' 3630 LIBPAW_ERROR(msg) 3631 end if 3632 3633 write(msg,'(a,i1)') & 3634 & ' Radial grid used for Vloc is grid ',ivlocmesh 3635 call wrtout(ab_out,msg,'COLL') 3636 call wrtout(std_out, msg,'COLL') 3637 3638 !------------------------------------------------- 3639 !Read LDA-1/2 potential 3640 3641 if (paw_setuploc%LDA_minus_half_potential%tread) then 3642 do imsh=1,nmesh 3643 if(trim(paw_setuploc%LDA_minus_half_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3644 iread1=imsh 3645 exit 3646 end if 3647 end do 3648 if(iread1/=ivlocmesh) then 3649 write(msg, '(a)' )& 3650 & 'The LDA-1/2 potential must be given on the same grid as the local potential.' 3651 LIBPAW_ERROR(msg) 3652 end if 3653 has_v_minushalf=1 3654 LIBPAW_ALLOCATE(pawtab%vminushalf,(vloc_mesh%mesh_size)) 3655 shft=mesh_shift(ivlocmesh) 3656 pawtab%vminus_mesh_size=vloc_mesh%mesh_size 3657 pawtab%vminushalf(1+shft:vloc_mesh%mesh_size)= & 3658 & paw_setuploc%LDA_minus_half_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3659 if (shft==1) call pawrad_deducer0(pawtab%vminushalf,vloc_mesh%mesh_size,vloc_mesh) 3660 write(msg,'(a,i1)') & 3661 & ' Radial grid used for LDA-1/2 potential is grid ',ivlocmesh 3662 call wrtout(ab_out,msg,'COLL') 3663 call wrtout(std_out, msg,'COLL') 3664 else 3665 has_v_minushalf=0 3666 end if 3667 if(has_v_minushalf==0.and.pawtab%has_vminushalf==1) then 3668 write(msg, '(a)' )& 3669 & 'The LDA-1/2 potential must be given in the XML PAW datafile.' 3670 LIBPAW_ERROR(msg) 3671 end if 3672 3673 !--------------------------------- 3674 !Eventually read "numeric" shapefunctions (if shape_type=-1) 3675 3676 if (pawtab%shape_type==-1) then 3677 LIBPAW_ALLOCATE(pawtab%shapefunc,(pawtab%mesh_size,pawtab%l_size)) 3678 do imsh=1,nmesh 3679 if(trim(paw_setuploc%shape_function%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3680 iread1=imsh 3681 exit 3682 end if 3683 end do 3684 call pawrad_copy(radmesh(iread1),shpf_mesh) 3685 ishpfmesh=iread1 3686 LIBPAW_ALLOCATE(shpf,(shpf_mesh%mesh_size,pawtab%l_size)) 3687 shft=mesh_shift(ishpfmesh) 3688 shpf(1,1)=one 3689 do ir=2,shpf_mesh%mesh_size 3690 shpf(ir,1)=paw_setuploc%shape_function%data(ir-shft,1) 3691 end do 3692 sz10=size(paw_setuploc%shape_function%data,2) 3693 if(sz10>=2) then 3694 do il=2,pawtab%l_size 3695 shpf(1,il)=zero 3696 do ir=2,shpf_mesh%mesh_size 3697 shpf(ir,il)=paw_setuploc%shape_function%data(ir-shft,il) 3698 end do 3699 end do 3700 else 3701 do il=2,pawtab%l_size 3702 shpf(1,il)=zero 3703 do ir=2,shpf_mesh%mesh_size 3704 shpf(ir,il)=paw_setuploc%shape_function%data(ir-shft,1)*shpf_mesh%rad(ir)**(il-1) 3705 end do 3706 end do 3707 end if 3708 write(msg,'(a,i1)') & 3709 & ' Radial grid used for shape functions is grid ',iread1 3710 call wrtout(ab_out,msg,'COLL') 3711 call wrtout(std_out, msg,'COLL') 3712 3713 ! Has to spline shape functions if mesh is not the "main" mesh 3714 if (ishpfmesh/=imainmesh) then 3715 msz=shpf_mesh%mesh_size 3716 LIBPAW_ALLOCATE(work1,(msz)) 3717 LIBPAW_ALLOCATE(work2,(msz)) 3718 LIBPAW_ALLOCATE(work3,(msz)) 3719 LIBPAW_ALLOCATE(work4,(pawrad%mesh_size)) 3720 work3(1:msz)=shpf_mesh%rad(1:msz) 3721 work4(1:pawrad%mesh_size)=pawrad%rad(1:pawrad%mesh_size) 3722 do il=1,pawtab%l_size 3723 call bound_deriv(shpf(1:msz,il),shpf_mesh,msz,yp1,ypn) 3724 call paw_spline(work3,shpf(:,il),msz,yp1,ypn,work1) 3725 call paw_splint(msz,work3,shpf(:,il),work1,pawrad%mesh_size,work4,pawtab%shapefunc(:,il)) 3726 end do 3727 LIBPAW_DEALLOCATE(work1) 3728 LIBPAW_DEALLOCATE(work2) 3729 LIBPAW_DEALLOCATE(work3) 3730 LIBPAW_DEALLOCATE(work4) 3731 else 3732 pawtab%shapefunc(:,:)=shpf(:,:) 3733 end if 3734 LIBPAW_DEALLOCATE(shpf) 3735 end if 3736 3737 !--------------------------------- 3738 !Read pseudo valence density 3739 3740 if (paw_setuploc%pseudo_valence_density%tread) then 3741 do imsh=1,nmesh 3742 if(trim(paw_setuploc%pseudo_valence_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3743 iread1=imsh 3744 exit 3745 end if 3746 end do 3747 ivalemesh=iread1 3748 call pawrad_copy(radmesh(iread1),vale_mesh) 3749 LIBPAW_ALLOCATE(tnvale,(vale_mesh%mesh_size)) 3750 shft=mesh_shift(ivalemesh) 3751 tnvale(1+shft:vale_mesh%mesh_size)=paw_setuploc%pseudo_valence_density%data(1:vale_mesh%mesh_size-shft)/sqrt(fourpi) 3752 if (shft==1) call pawrad_deducer0(tnvale,vale_mesh%mesh_size,vale_mesh) 3753 pawtab%has_tvale=1 3754 write(msg,'(a,i1)') & 3755 & ' Radial grid used for pseudo valence density is grid ',ivalemesh 3756 call wrtout(ab_out,msg,'COLL') 3757 call wrtout(std_out, msg,'COLL') 3758 else 3759 pawtab%has_tvale=0 3760 LIBPAW_ALLOCATE(tnvale,(0)) 3761 end if 3762 3763 !--------------------------------- 3764 !Read initial guess of rhoij (rhoij0) 3765 3766 LIBPAW_ALLOCATE(pawtab%rhoij0,(pawtab%lmn2_size)) 3767 pawtab%rhoij0=zero 3768 ilmn0=0 3769 do ib=1,pawtab%basis_size 3770 il=2*pawtab%orbitals(ib)+1 3771 occ=paw_setuploc%valence_states%state(ib)%ff 3772 if (occ<zero)occ=zero 3773 do ilmn=ilmn0+1,ilmn0+il 3774 pawtab%rhoij0(ilmn*(ilmn+1)/2)=occ/dble(il) 3775 end do 3776 ilmn0=ilmn0+il 3777 end do 3778 3779 !--------------------------------- 3780 !Read Kij terms (kij0) and deduce eventually Dij0 3781 3782 LIBPAW_ALLOCATE(kij,(pawtab%lmn2_size)) 3783 kij=zero 3784 nval=paw_setuploc%valence_states%nval 3785 do jlmn=1,pawtab%lmn_size 3786 j0lmn=jlmn*(jlmn-1)/2 3787 jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn) 3788 do ilmn=1,jlmn 3789 klmn=j0lmn+ilmn 3790 ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn) 3791 if (ilm==jlm) kij(klmn)=paw_setuploc%kinetic_energy_differences%data(jln+(iln-1)*nval) 3792 end do 3793 end do 3794 if (vlocopt>0) then 3795 LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size)) 3796 if (allocated(pawtab%vminushalf).and.pawtab%has_vminushalf==1) then 3797 vlocr(1:vloc_mesh%mesh_size)=vlocr(1:vloc_mesh%mesh_size)+pawtab%vminushalf(1:vloc_mesh%mesh_size) 3798 end if 3799 call atompaw_dij0(pawtab%indlmn,kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,& 3800 & vloc_mesh,vlocr,znucl) 3801 end if 3802 3803 !Keep eventualy Kij in memory 3804 if (pawtab%has_kij==1.or.vlocopt==0) then 3805 LIBPAW_ALLOCATE(pawtab%kij,(pawtab%lmn2_size)) 3806 pawtab%kij(:)=kij(:) 3807 if (vlocopt> 0) pawtab%has_kij=2 3808 ! This -1 means that pawtab%kij will be freed later 3809 if (vlocopt==0) pawtab%has_kij=-1 3810 end if 3811 3812 LIBPAW_DEALLOCATE(kij) 3813 3814 !--------------------------------- 3815 !Read exact-exchange Fock terms for core-valence interactions (ex_cvij) 3816 3817 if (paw_setuploc%exact_exchange_matrix%tread.eqv..true.) then 3818 pawtab%has_fock=2 3819 LIBPAW_ALLOCATE(pawtab%ex_cvij,(pawtab%lmn2_size)) 3820 pawtab%ex_cvij=zero 3821 nval=paw_setuploc%valence_states%nval 3822 do jlmn=1,pawtab%lmn_size 3823 j0lmn=jlmn*(jlmn-1)/2 3824 jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn) 3825 do ilmn=1,jlmn 3826 klmn=j0lmn+ilmn 3827 ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn) 3828 if (ilm==jlm) pawtab%ex_cvij(klmn)=paw_setuploc%exact_exchange_matrix%data(jln+(iln-1)*nval) 3829 end do 3830 end do 3831 pawtab%ex_cc=paw_setuploc%ex_cc 3832 end if 3833 3834 !---------------------------------------- 3835 ! store Lamb shielding 3836 pawtab%lamb_shielding=paw_setuploc%lamb_shielding 3837 3838 !========================================================== 3839 !Compute additional atomic data only depending on present DATASET 3840 3841 call pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,hyb_mixing,ixc,lnmax,& 3842 & mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,& 3843 & qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,& 3844 & vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl,& 3845 & tcoretau=tcoretau,coretau_mesh=coretau_mesh,xc_taupos=my_xc_taupos) 3846 3847 if(usewvl==1 .or. icoulomb > 0) then 3848 ! Calculate up to the 5th derivative of tcoredens 3849 call pawpsp_calc_d5(core_mesh,pawtab%core_mesh_size,pawtab%tcoredens) 3850 ! Other wvl related operations 3851 call pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr) 3852 else if (pawtab%has_wvl>0) then 3853 call wvlpaw_rholoc_nullify(pawtab%wvl%rholoc) 3854 end if 3855 3856 !========================================================== 3857 !Free temporary allocated space 3858 3859 call pawrad_free(radmesh) 3860 LIBPAW_DATATYPE_DEALLOCATE(radmesh) 3861 LIBPAW_DEALLOCATE(mesh_shift) 3862 3863 call pawrad_free(tproj_mesh) 3864 call pawrad_free(core_mesh) 3865 call pawrad_free(vloc_mesh) 3866 call pawrad_free(coretau_mesh) 3867 3868 if (allocated(vlocr)) then 3869 LIBPAW_DEALLOCATE(vlocr) 3870 end if 3871 if (allocated(ncore)) then 3872 LIBPAW_DEALLOCATE(ncore) 3873 end if 3874 if (allocated(tncore)) then 3875 LIBPAW_DEALLOCATE(tncore) 3876 end if 3877 if (allocated(coretau)) then 3878 LIBPAW_DEALLOCATE(coretau) 3879 end if 3880 if (allocated(tcoretau)) then 3881 LIBPAW_DEALLOCATE(tcoretau) 3882 end if 3883 if (allocated(tproj)) then 3884 LIBPAW_DEALLOCATE(tproj) 3885 end if 3886 3887 if(pawtab%shape_type==-1) then 3888 call pawrad_free(shpf_mesh) 3889 end if 3890 if (paw_setuploc%pseudo_valence_density%tread) then 3891 call pawrad_free(vale_mesh) 3892 end if 3893 if (paw_setuploc%ae_core_kinetic_energy_density%tread.and.pawtab%has_coretau>=1) then 3894 call pawrad_free(coretau_mesh) 3895 end if 3896 if (allocated(tnvale)) then 3897 LIBPAW_DEALLOCATE(tnvale) 3898 end if 3899 3900 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 xc_denpos= lowest allowed density (usually for the computation of the XC functionals) [xc_taupos]= lowest allowed kinetic energy density (for mGGA XC functionals) 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
xcccrc=XC core correction cutoff radius (bohr) from psp file
NOTES
Spin-orbit not yet implemented (to be done)
SOURCE
3940 subroutine pawpsp_7in(epsatm,ffspl,icoulomb,hyb_mixing,ixc,& 3941 & lmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,& 3942 & pawrad,pawtab,pawxcdev,qgrid_ff,qgrid_vl,& 3943 & usewvl,usexcnhat_in,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl,& 3944 & xc_taupos) ! Optional argument 3945 3946 !Arguments ------------------------------------ 3947 !scalars 3948 integer, intent(in):: icoulomb,ixc 3949 integer, intent(in):: lmax,lnmax,mmax 3950 integer, intent(in):: mqgrid_ff,mqgrid_vl,pawxcdev 3951 integer, intent(in):: usewvl,usexcnhat_in,xclevel 3952 real(dp), intent(in):: hyb_mixing,xc_denpos,zion,znucl 3953 real(dp), intent(in),optional:: xc_taupos 3954 real(dp), intent(out):: epsatm,xcccrc 3955 type(pawrad_type), intent(inout):: pawrad 3956 type(pawtab_type), intent(inout) :: pawtab 3957 !arrays 3958 real(dp),intent(in):: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl) 3959 real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax) 3960 real(dp),intent(out) :: vlspl(mqgrid_vl,2) 3961 3962 !Local variables ------------------------------ 3963 !scalars 3964 integer :: imainmesh,nmesh 3965 integer :: pspversion,usexcnhat,vlocopt 3966 logical :: save_core_msz 3967 real(dp) :: my_xc_taupos 3968 type(pawrad_type) :: core_mesh,tproj_mesh,vale_mesh,vloc_mesh 3969 !arrays 3970 real(dp),pointer :: ncore(:),tncore(:),tcoretau(:),tnvale(:),tproj(:,:),vlocr(:) 3971 type(pawrad_type),pointer :: radmesh(:) 3972 3973 !************************************************************************ 3974 3975 !Destroy everything in pawtab but optional flags 3976 call pawtab_free(pawtab) 3977 !Destroy everything in pawrad 3978 call pawrad_free(pawrad) 3979 3980 save_core_msz=(usewvl==1 .or. icoulomb .ne. 0) 3981 nullify(ncore);nullify(tncore);nullify(tcoretau);nullify(tnvale) 3982 nullify(tproj);nullify(vlocr) 3983 nullify(radmesh) 3984 3985 call pawpsp_read(core_mesh,tmp_unit,imainmesh,lmax,& 3986 & ncore,nmesh,pawrad,pawtab,pspversion,radmesh,save_core_msz,& 3987 & tcoretau,tncore,tnvale,tproj,tproj_mesh,usexcnhat_in,usexcnhat,& 3988 & vale_mesh,vlocopt,vlocr,vloc_mesh,znucl) 3989 3990 my_xc_taupos=xc_denpos;if(present(xc_taupos)) my_xc_taupos=xc_taupos 3991 call pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,hyb_mixing,ixc,lnmax,& 3992 & mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,& 3993 & qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,& 3994 & vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl,& 3995 & tcoretau=tcoretau,coretau_mesh=core_mesh,xc_taupos=my_xc_taupos) 3996 3997 if(usewvl==1 .or. icoulomb > 0) then 3998 ! Calculate up to the 5th derivative of tcoredens 3999 call pawpsp_calc_d5(core_mesh,pawtab%core_mesh_size,pawtab%tcoredens) 4000 ! Other wvl related operations 4001 call pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr) 4002 else if (pawtab%has_wvl>0) then 4003 call wvlpaw_rholoc_nullify(pawtab%wvl%rholoc) 4004 end if 4005 4006 !========================================================== 4007 !Free temporary allocated space 4008 call pawrad_free(radmesh) 4009 call pawrad_free(tproj_mesh) 4010 call pawrad_free(core_mesh) 4011 call pawrad_free(vloc_mesh) 4012 LIBPAW_DATATYPE_DEALLOCATE(radmesh) 4013 if (associated(vlocr)) then 4014 LIBPAW_POINTER_DEALLOCATE(vlocr) 4015 end if 4016 if (associated(ncore)) then 4017 LIBPAW_POINTER_DEALLOCATE(ncore) 4018 end if 4019 if (associated(tncore)) then 4020 LIBPAW_POINTER_DEALLOCATE(tncore) 4021 end if 4022 if (associated(tnvale)) then 4023 LIBPAW_POINTER_DEALLOCATE(tnvale) 4024 end if 4025 if (associated(tcoretau)) then 4026 LIBPAW_POINTER_DEALLOCATE(tcoretau) 4027 end if 4028 if (associated(tproj)) then 4029 LIBPAW_POINTER_DEALLOCATE(tproj) 4030 end if 4031 if (pspversion>=4) then 4032 call pawrad_free(vale_mesh) 4033 end if 4034 4035 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
SOURCE
4721 subroutine pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc) 4722 4723 !Arguments ------------------------------------ 4724 integer,intent(in) :: comm_mpi 4725 real(dp),intent(inout) :: epsatm,xcccrc 4726 real(dp),intent(inout) :: ffspl(:,:,:),vlspl(:,:) 4727 type(pawrad_type),intent(inout) :: pawrad 4728 type(pawtab_type),intent(inout) :: pawtab 4729 4730 !Local variables------------------------------- 4731 integer :: ierr,ii,me,nn_dpr 4732 integer :: siz_ffspl,siz1_ffspl,siz2_ffspl,siz3_ffspl,siz_vlspl,siz1_vlspl,siz2_vlspl 4733 integer,allocatable :: list_int(:) 4734 real(dp),allocatable :: list_dpr(:) 4735 4736 !************************************************************************* 4737 4738 me=xmpi_comm_rank(comm_mpi) 4739 4740 !Broadcast pawrad 4741 call pawrad_bcast(pawrad,comm_mpi) 4742 4743 !Broadcast pawtab (only data read from file) 4744 call pawtab_bcast(pawtab,comm_mpi,only_from_file=.true.) 4745 4746 !Broadcast the sizes of the arrays 4747 LIBPAW_ALLOCATE(list_int,(5)) 4748 if (me==0) then 4749 siz1_vlspl=size(vlspl,1); list_int(1)=siz1_vlspl 4750 siz2_vlspl=size(vlspl,2); list_int(2)=siz2_vlspl 4751 siz1_ffspl=size(ffspl,1); list_int(3)=siz1_ffspl 4752 siz2_ffspl=size(ffspl,2); list_int(4)=siz2_ffspl 4753 siz3_ffspl=size(ffspl,3); list_int(5)=siz3_ffspl 4754 end if 4755 call xmpi_bcast(list_int,0,comm_mpi,ierr) 4756 if (me/=0) then 4757 siz1_vlspl=list_int(1) 4758 siz2_vlspl=list_int(2) 4759 siz1_ffspl=list_int(3) 4760 siz2_ffspl=list_int(4) 4761 siz3_ffspl=list_int(5) 4762 end if 4763 siz_vlspl=siz1_vlspl*siz2_vlspl 4764 siz_ffspl=siz1_ffspl*siz2_ffspl*siz3_ffspl 4765 LIBPAW_DEALLOCATE(list_int) 4766 4767 !Broadcast the reals 4768 nn_dpr=2+siz_vlspl+siz_ffspl 4769 LIBPAW_ALLOCATE(list_dpr,(nn_dpr)) 4770 if (me==0) then 4771 ii=1 4772 list_dpr(ii)=epsatm ;ii=ii+1 4773 list_dpr(ii)=xcccrc ;ii=ii+1 4774 list_dpr(ii:ii+siz_vlspl-1)=reshape(vlspl,(/siz_vlspl/)) ;ii=ii+siz_vlspl 4775 list_dpr(ii:ii+siz_ffspl-1)=reshape(ffspl,(/siz_ffspl/)) ;ii=ii+siz_ffspl 4776 end if 4777 call xmpi_bcast(list_dpr,0,comm_mpi,ierr) 4778 if (me/=0) then 4779 ii=1 4780 epsatm=list_dpr(ii) ;ii=ii+1 4781 xcccrc=list_dpr(ii) ;ii=ii+1 4782 vlspl=reshape(list_dpr(ii:ii+siz_vlspl-1),(/siz1_vlspl,siz2_vlspl/)) 4783 ii=ii+siz_vlspl 4784 ffspl=reshape(list_dpr(ii:ii+siz_ffspl-1),(/siz1_ffspl,siz2_ffspl,siz3_ffspl/)) 4785 ii=ii+siz_ffspl 4786 end if 4787 LIBPAW_DEALLOCATE(list_dpr) 4788 4789 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 [coretau_mesh<type(pawrad_type)>]=radial mesh for the core kinetic energy 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 [tcoretau(coretau_mesh%mesh_size)]= pseudo core kinetic energy 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) [xc_taupos]= lowest allowed kinetic energy density (for mGGA 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
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
SOURCE
1770 subroutine pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,hyb_mixing,ixc,lnmax,& 1771 & mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,& 1772 & qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,& 1773 & vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl,& 1774 & tcoretau,coretau_mesh,xc_taupos) !optional 1775 1776 !Arguments ------------------------------------ 1777 !scalars 1778 integer,intent(in) :: imainmesh,ixc,lnmax,mqgrid_ff,mqgrid_vl 1779 integer,intent(in) :: nmesh,pawxcdev,pspversion,usexcnhat,vlocopt 1780 integer,intent(in) ::mmax 1781 integer,intent(in) :: xclevel 1782 real(dp),intent(in) :: hyb_mixing,xc_denpos,zion,znucl 1783 real(dp),intent(in),optional :: xc_taupos 1784 real(dp),intent(out) :: epsatm,xcccrc 1785 type(pawrad_type),intent(in) :: core_mesh,tproj_mesh,vale_mesh 1786 type(pawrad_type),intent(in),optional :: coretau_mesh 1787 type(pawrad_type),intent(inout) ::pawrad,vloc_mesh 1788 type(pawtab_type),intent(inout) :: pawtab 1789 !arrays 1790 real(dp),intent(in) :: ncore(core_mesh%mesh_size),tncore(core_mesh%mesh_size) 1791 real(dp),intent(in),optional :: tcoretau(:) 1792 real(dp),intent(in) :: qgrid_vl(mqgrid_vl),qgrid_ff(mqgrid_ff) 1793 real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax) 1794 real(dp),intent(out) :: vlspl(mqgrid_vl,2) 1795 real(dp),intent(inout) :: tnvale(vale_mesh%mesh_size*pawtab%has_tvale) 1796 real(dp),intent(inout) :: tproj(tproj_mesh%mesh_size,pawtab%basis_size) 1797 real(dp),intent(inout) :: vlocr(vloc_mesh%mesh_size) 1798 type(pawrad_type),intent(in) :: radmesh(nmesh) 1799 1800 !Local variables ------------------------------ 1801 !scalars 1802 integer,parameter :: reduced_mshsz=2501 1803 integer :: ib,il,ilm,ilmn,iln,ir,isnotzero,itest 1804 integer :: j0lmn,jlm,jlmn,jln,klmn,msz,msz1,msz_tmp,mst_tmp,nspden,usekden 1805 logical :: has_dij0,non_magnetic_xc,reduced_ncor,reduced_taucor,reduced_nval,reduced_vloc,testval 1806 real(dp),parameter :: reduced_rstep=0.00025_dp,rm_vloc=20.0_dp 1807 real(dp) :: d2nvdq0,intg,intvh,lstep_tmp,my_xc_taupos,qcore,qq,rstep_tmp,yp1,ypn 1808 character(len=500) :: msg 1809 type(pawang_type) :: pawang_tmp 1810 type(pawrad_type) :: rcore_mesh,rcoretau_mesh,rvale_mesh,rvloc_mesh,tproj_mesh_new 1811 !arrays 1812 real(dp) :: tmp_qgrid(1),tmp_q2vq(1) 1813 real(dp),allocatable :: ncorwk(:),nhat(:),nhatwk(:),nwk(:),r2k(:) 1814 real(dp),allocatable :: rtncor(:),rttaucor(:),rtnval(:),rvlocr(:) 1815 real(dp),allocatable :: vbare(:),vh(:),vhnzc(:),vxc1(:),vxc2(:) 1816 real(dp),allocatable,target :: work1(:),work2(:),work3(:) 1817 real(dp),pointer :: tmp1(:),tmp2(:) 1818 logical :: tmp_lmselect(1) 1819 1820 ! ************************************************************************* 1821 1822 !========================================================== 1823 !Perfom tests on meshes 1824 1825 !Are radial meshes for Phi and Vloc compatibles ? 1826 ! if (vloc_mesh%rmax<pawrad%rmax) then 1827 ! write(msg, '(a,a,a)' )& 1828 !& 'Rmax for Vloc < Rmax !',ch10,& 1829 !& 'Action : check your pseudopotential (increase Vloc meshSize).' 1830 ! LIBPAW_ERROR(msg) 1831 ! end if 1832 1833 !Check optional arguments 1834 usekden=merge(0,1,pawtab%has_coretau==0) 1835 my_xc_taupos=xc_denpos;if (present(xc_taupos)) my_xc_taupos=xc_taupos 1836 if (present(tcoretau)) then 1837 if (usekden>=1) then 1838 if (.not.(present(coretau_mesh))) then 1839 msg='tcoretau present but not coretau_mesh!' 1840 LIBPAW_BUG(msg) 1841 end if 1842 if (size(tcoretau)>coretau_mesh%mesh_size) then 1843 msg='wrong size for tcoretau!' 1844 LIBPAW_BUG(msg) 1845 end if 1846 if (coretau_mesh%mesh_size<pawtab%mesh_size) then 1847 write(msg, '(a,a,a,a,a)' )& 1848 & 'Mesh size for core kinetic energy density must be equal or larger',ch10,& 1849 & 'than mesh size for PAW augmentation regions !',ch10,& 1850 & 'Action : check your pseudopotential (increase TAUcore meshSize).' 1851 LIBPAW_ERROR(msg) 1852 end if 1853 end if 1854 end if 1855 1856 !Are mmax and mesh_size for partial waves compatibles ? 1857 if (mmax/=pawtab%partialwave_mesh_size) then 1858 write(msg, '(a,a,a)' )& 1859 & 'mmax /= phi_mesh_size in psp file !',ch10,& 1860 & 'Action: check your pseudopotential file.' 1861 LIBPAW_ERROR(msg) 1862 end if 1863 1864 !Are radial meshes for (t)Ncore / tTAU and Phi compatibles ? 1865 if (usekden>=1.and.core_mesh%mesh_size<pawtab%mesh_size) then 1866 write(msg, '(a,a,a,a,a)' )& 1867 & 'Mesh size for core density must be equal or larger',ch10,& 1868 & 'than mesh size for PAW augmentation regions !',ch10,& 1869 & 'Action : check your pseudopotential (increase Ncore meshSize).' 1870 LIBPAW_ERROR(msg) 1871 end if 1872 1873 !Are radial meshes for (t)Nvale and Phi compatibles ? 1874 if ((pawtab%has_tvale==1).and.(vale_mesh%rmax<pawrad%rmax)) then 1875 write(msg, '(a,a,a)' )& 1876 & 'Rmax for tNvale < Rmax for Phi !',ch10,& 1877 & 'Action : check your pseudopotential (increase tNvale meshSize).' 1878 LIBPAW_ERROR(msg) 1879 end if 1880 1881 !Is PAW radius included inside radial mesh ? 1882 if (pawtab%rpaw>pawrad%rmax+tol8) then 1883 write(msg, '(a,a,a)' )& 1884 & 'Radius of PAW sphere is outside the radial mesh !',ch10,& 1885 & 'Action: check your pseudopotential file.' 1886 LIBPAW_ERROR(msg) 1887 end if 1888 1889 !Max. radius of mesh for Vloc has to be "small" in order to avoid numeric noise ? 1890 if (vloc_mesh%rmax>rm_vloc) then 1891 msz_tmp=pawrad_ifromr(vloc_mesh,rm_vloc);mst_tmp=vloc_mesh%mesh_type 1892 rstep_tmp=vloc_mesh%rstep;lstep_tmp=vloc_mesh%lstep 1893 call pawrad_free(vloc_mesh) 1894 call pawrad_init(vloc_mesh,mesh_size=msz_tmp,mesh_type=mst_tmp,& 1895 & rstep=rstep_tmp,lstep=lstep_tmp,r_for_intg=rm_vloc) 1896 write(msg, '(a,i4,a)' ) ' Mesh size for Vloc has been set to ', & 1897 & vloc_mesh%mesh_size,' to avoid numerical noise.' 1898 call wrtout(std_out,msg,'COLL') 1899 call wrtout(ab_out,msg,'COLL') 1900 end if 1901 1902 !This test has been disable... MT 2006-25-10 1903 !For Simpson rule, it is better to have odd mesh sizes 1904 !itest=0 1905 !do imsh=1,nmesh 1906 !if (mod(radmesh(imsh)%mesh_size,2)==0.and.radmesh(imsh)%mesh_type==1) itest=1 1907 !end do 1908 !if (itest==1) then 1909 ! write(msg, '(5a)' ) & 1910 !& 'Regular radial meshes should have odd number of points ',ch10,& 1911 !& 'for better accuracy of integration sheme (Simpson rule).',ch10,& 1912 !& 'Althought it''s not compulsory, you should change mesh sizes in psp file.' 1913 ! LIBPAW_WARNING(msg) 1914 !end if 1915 1916 !Test the compatibilty between Rpaw and mesh for (t)Phi 1917 if (pspversion>=3) then 1918 itest=pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw) 1919 ! This test has been disable... MT 2015-02-12 1920 ! if (itest+2>radmesh(imainmesh)%mesh_size) then 1921 ! write(msg, '(9a)' ) & 1922 !& 'Atomic data could produce inaccurate results:',ch10,& 1923 !& 'Wavefunctions and pseudo-wavefunctions should',ch10,& 1924 !& 'be given on a radial mesh larger than the PAW',ch10,& 1925 !& 'spheres (at least 2 additional points) !',ch10,& 1926 !& 'Action: check your pseudopotential file.' 1927 ! LIBPAW_WARNING(msg) 1928 ! end if 1929 if (abs(pawtab%rpaw-radmesh(imainmesh)%rad(itest))<tol8) itest=itest-1 1930 ib=0;isnotzero=0 1931 do while ((isnotzero==0).and.(ib<pawtab%basis_size)) 1932 ib=ib+1;ir=itest 1933 do while ((isnotzero==0).and.(ir<pawtab%mesh_size)) 1934 ir=ir+1;if (abs(pawtab%phi(ir,ib)-pawtab%tphi(ir,ib))>tol8) isnotzero=1 1935 end do 1936 end do 1937 if (isnotzero>0) then 1938 write(msg, '(7a)' )& 1939 & 'Atomic data are inconsistent:',ch10,& 1940 & 'For r>=r_paw, pseudo wavefunctions are not',ch10,& 1941 & 'equal to wave functions (Phi(r)/=tPhi(r)) !',ch10,& 1942 & 'Action: check your pseudopotential file.' 1943 LIBPAW_ERROR(msg) 1944 end if 1945 else 1946 ! For compatibility reasons set PAW radius at the end of mesh (older versions) 1947 if (pawtab%rpaw/=pawrad%rmax) then 1948 msz_tmp=pawrad%mesh_size;mst_tmp=pawrad%mesh_type 1949 rstep_tmp=pawrad%rstep;lstep_tmp=pawrad%lstep 1950 call pawrad_free(pawrad) 1951 call pawrad_init(pawrad,mesh_size=msz_tmp,mesh_type=mst_tmp,rstep=rstep_tmp,lstep=lstep_tmp) 1952 pawtab%rpaw=pawrad%rmax 1953 end if 1954 end if 1955 !If Vloc is a "Vbare" potential, it has to be localized inside PAW spheres 1956 if (vlocopt==0.and.(vloc_mesh%rmax>pawtab%rpaw+tol10)) then 1957 if(vlocr(pawrad_ifromr(vloc_mesh,pawtab%rpaw))>tol10) then 1958 write(msg, '(7a)' )& 1959 & 'Atomic data are inconsistent:',ch10,& 1960 & 'Local potential is a "Vbare" potential',ch10,& 1961 & 'and is not localized inside PAW sphere !',ch10,& 1962 & 'Vbare is set to zero if r>rpaw.' 1963 LIBPAW_WARNING(msg) 1964 do ir=pawrad_ifromr(vloc_mesh,pawtab%rpaw),vloc_mesh%mesh_size 1965 vlocr(ir)=zero 1966 end do 1967 end if 1968 end if 1969 1970 !========================================================== 1971 !Initializations 1972 1973 has_dij0=(allocated(pawtab%dij0)) 1974 1975 !Allocate/initialize some dummy variables 1976 tmp_lmselect(1)=.true. 1977 non_magnetic_xc=.false. 1978 if (pawxcdev==0) then 1979 pawang_tmp%l_size_max=1;pawang_tmp%angl_size=1;pawang_tmp%ylm_size=1 1980 pawang_tmp%use_ls_ylm=0;pawang_tmp%gnt_option=0;pawang_tmp%ngnt=0;pawang_tmp%nsym=0 1981 LIBPAW_ALLOCATE(pawang_tmp%angwgth,(1)) 1982 pawang_tmp%angwgth(1)=one 1983 LIBPAW_ALLOCATE(pawang_tmp%anginit,(3,1)) 1984 pawang_tmp%anginit(1,1)=one 1985 pawang_tmp%anginit(2:3,1)=zero 1986 LIBPAW_ALLOCATE(pawang_tmp%ylmr,(1,1)) 1987 pawang_tmp%ylmr(1,1)=1._dp/sqrt(four_pi) 1988 LIBPAW_ALLOCATE(pawang_tmp%ylmrgr,(9,1,1)) 1989 pawang_tmp%ylmrgr(1:9,1,1)=zero 1990 end if 1991 1992 !========================================================== 1993 !Compute ffspl(q) (and derivatives) 1994 1995 ffspl=zero 1996 if (mqgrid_ff>0) then 1997 call pawpsp_nl(ffspl,pawtab%indlmn,pawtab%lmn_size,lnmax,mqgrid_ff,qgrid_ff,& 1998 & tproj_mesh,tproj) 1999 end if 2000 2001 !========================================================== 2002 !Compute eventually compensation charge radius (i.e. radius for shape functions) 2003 2004 if (pawtab%shape_type>0.and.pawtab%rshp<1.d-8) then 2005 pawtab%rshp=pawtab%rpaw 2006 else if (pawtab%shape_type==-1) then 2007 ir=pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw)+1;isnotzero=0 2008 do while ((isnotzero==0).and.(ir>1)) 2009 ir=ir-1;il=0 2010 do while ((isnotzero==0).and.(il<pawtab%l_size)) 2011 il=il+1;if (pawtab%shapefunc(ir,il)>tol16) isnotzero=1 2012 end do 2013 end do 2014 ir=min(ir+1,pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw)) 2015 pawtab%rshp=radmesh(imainmesh)%rad(ir) 2016 do il=1,pawtab%l_size 2017 if (pawtab%shapefunc(ir,il)>tol6) then 2018 write(msg, '(a,a,a)' )& 2019 & 'Shape function is not zero at PAW radius !',ch10,& 2020 & 'Action: check your pseudopotential file.' 2021 LIBPAW_ERROR(msg) 2022 end if 2023 end do 2024 end if 2025 2026 !========================================================== 2027 !Compute compensation charge density (nhat) 2028 !Add it to pseudo valence density 2029 2030 if (pawtab%has_tvale==1) then 2031 msz=vale_mesh%mesh_size 2032 LIBPAW_ALLOCATE(nhat,(msz)) 2033 ! A-Has to compute norm of nhat (Int[n-tild_n]) 2034 testval=(abs(tnvale(msz))<tol9) 2035 ! A1-If tnvale is not given with enough points, 2036 ! try to compute it from rhoij0 and tphi 2037 if (.not.testval) then 2038 msz1=pawtab%mesh_size 2039 ! Compute n and tild_n from phi and tphi 2040 LIBPAW_ALLOCATE(work1,(msz1)) 2041 LIBPAW_ALLOCATE(work2,(msz1)) 2042 work1=zero 2043 work2=zero 2044 do jlmn=1,pawtab%lmn_size 2045 j0lmn=jlmn*(jlmn-1)/2;jln=pawtab%indlmn(5,jlmn) 2046 do ilmn=1,jlmn 2047 klmn=j0lmn+ilmn;iln=pawtab%indlmn(5,ilmn) 2048 yp1=two;if (ilmn==jlmn) yp1=one 2049 work1(1:msz1)=work1(1:msz1)+yp1*pawtab%rhoij0(klmn) & 2050 & *pawtab% phi(1:msz1,iln)*pawtab% phi(1:msz1,jln) 2051 work2(1:msz1)=work2(1:msz1)+yp1*pawtab%rhoij0(klmn) & 2052 & *pawtab%tphi(1:msz1,iln)*pawtab%tphi(1:msz1,jln) 2053 end do 2054 end do 2055 ! Spline tnvale onto pawrad if needed 2056 LIBPAW_ALLOCATE(nwk,(msz1)) 2057 if ((vale_mesh%mesh_type/=pawrad%mesh_type).or.(vale_mesh%rstep/=pawrad%rstep).or.& 2058 & (vale_mesh%lstep/=pawrad%lstep)) then 2059 LIBPAW_ALLOCATE(work3,(vale_mesh%mesh_size)) 2060 call bound_deriv(tnvale(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn) 2061 call paw_spline(vale_mesh%rad,tnvale,vale_mesh%mesh_size,yp1,ypn,work3) 2062 call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,tnvale,work3,msz1,pawrad%rad(1:msz1),nwk(1:msz1)) 2063 LIBPAW_DEALLOCATE(work3) 2064 else 2065 nwk(1:msz1)=tnvale(1:msz1) 2066 end if 2067 ! Compare tild_n and tnvale (inside aug. region) 2068 if (maxval(abs((nwk(1:msz1)*four_pi*pawrad%rad(1:msz1)**2)-work2(1:msz1)))<tol6) then 2069 ! If equality then compute Int[n-tild_n] 2070 work1=work1-work2 2071 call simp_gen(qq,work1,pawrad) 2072 qq=qq/four_pi 2073 else 2074 ! If not equality, will use tnvale 2075 testval=.true. 2076 write(msg, '(3a)' ) & 2077 & 'Valence density is not given with enough points',ch10,& 2078 & 'in psp file. Some charge estimations will be coarse.' 2079 LIBPAW_WARNING(msg) 2080 end if 2081 LIBPAW_DEALLOCATE(nwk) 2082 LIBPAW_DEALLOCATE(work1) 2083 LIBPAW_DEALLOCATE(work2) 2084 end if 2085 ! A2-If tnvale is given with enough points, use it 2086 if (testval) then 2087 nhat(1:msz)=tnvale(1:msz)*vale_mesh%rad(1:msz)**2 2088 call simp_gen(qq,nhat,vale_mesh) 2089 qq=zion/four_pi-qq 2090 end if 2091 ! B-Compute nhat and add it to pseudo valence density 2092 call atompaw_shpfun(0,vale_mesh,intg,pawtab,nhat) 2093 nhat(1:msz)=qq*nhat(1:msz) 2094 tnvale(1:msz)=tnvale(1:msz)+nhat(1:msz) 2095 end if 2096 2097 !========================================================== 2098 !If Vloc potential is in "Vbare" format, translate it into VH(tnzc) format 2099 2100 if (vlocopt==0) then 2101 write(msg,'(a)') ' Local potential is in "Vbare" format... ' 2102 call wrtout(ab_out,msg,'COLL') 2103 call wrtout(std_out, msg,'COLL') 2104 msz=core_mesh%mesh_size 2105 LIBPAW_ALLOCATE(r2k,(msz)) 2106 call atompaw_shpfun(0,core_mesh,intg,pawtab,r2k) 2107 r2k(1:msz)=r2k(1:msz)*core_mesh%rad(1:msz)**2 2108 ! Compute VH[4pi.r2.n(r)=4pi.r2.tncore(r)+(Qcore-Z).r2.k(r)] 2109 LIBPAW_ALLOCATE(nwk,(core_mesh%mesh_size)) 2110 LIBPAW_ALLOCATE(vh,(core_mesh%mesh_size)) 2111 if (core_mesh%mesh_type==5) then 2112 nwk(1:msz)=tncore(1:msz)*four_pi*core_mesh%rad(1:msz)**2 2113 call simp_gen(qcore,nwk,core_mesh) 2114 qcore=znucl-zion-qcore 2115 else 2116 nwk(1:msz)=(ncore(1:msz)-tncore(1:msz))*four_pi*core_mesh%rad(1:msz)**2 2117 ib=1 2118 do ir=msz,2,-1 2119 if(abs(nwk(ir))<tol14)ib=ir 2120 end do 2121 call simp_gen(qcore,nwk,core_mesh,r_for_intg=core_mesh%rad(ib)) 2122 nwk(1:msz)=tncore(1:msz)*four_pi*core_mesh%rad(1:msz)**2 2123 end if 2124 nwk(1:msz)=nwk(1:msz)+r2k(1:msz)*(qcore-znucl) 2125 call poisson(nwk,0,core_mesh,vh) 2126 vh(2:msz)=vh(2:msz)/core_mesh%rad(2:msz) 2127 call pawrad_deducer0(vh,msz,core_mesh) 2128 2129 LIBPAW_DEALLOCATE(nwk) 2130 ! Eventually spline Vbare 2131 LIBPAW_ALLOCATE(vbare,(core_mesh%mesh_size)) 2132 if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2133 & (core_mesh%rstep /=vloc_mesh%rstep) .or.& 2134 & (core_mesh%lstep /=vloc_mesh%lstep)) then 2135 msz=core_mesh%mesh_size;if (vloc_mesh%rmax<core_mesh%rmax) msz=pawrad_ifromr(core_mesh,vloc_mesh%rmax) 2136 call bound_deriv(vlocr(1:vloc_mesh%mesh_size),vloc_mesh,vloc_mesh%mesh_size,yp1,ypn) 2137 LIBPAW_ALLOCATE(work1,(vloc_mesh%mesh_size)) 2138 LIBPAW_ALLOCATE(work2,(vloc_mesh%mesh_size)) 2139 call paw_spline(vloc_mesh%rad,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1) 2140 call paw_splint(vloc_mesh%mesh_size,vloc_mesh%rad,vlocr,work1,msz,core_mesh%rad(1:msz),vbare) 2141 LIBPAW_DEALLOCATE(work1) 2142 LIBPAW_DEALLOCATE(work2) 2143 else 2144 msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size) 2145 vbare(1:msz)=vlocr(1:msz) 2146 end if 2147 ! Build VH(tnzc) from Vbare 2148 vlocr(1:msz)=vbare(1:msz)+vh(1:msz) 2149 if(vloc_mesh%mesh_size>msz)then 2150 vlocr(msz+1:vloc_mesh%mesh_size)=vh(msz)*vloc_mesh%rad(msz)/vloc_mesh%rad(msz+1:vloc_mesh%mesh_size) 2151 end if 2152 LIBPAW_DEALLOCATE(vbare) 2153 LIBPAW_DEALLOCATE(vh) 2154 2155 ! Compute <tPhi_i|VH(tnzc)|tPhi_j> and int[VH(tnzc)*Qijhat(r)dr] parts of Dij0 2156 ! Note: it is possible as core_mesh and radmesh(imainmesh) have the same steps 2157 if (has_dij0) then 2158 msz=radmesh(imainmesh)%mesh_size 2159 LIBPAW_ALLOCATE(work1,(msz)) 2160 work1(1:msz)=vlocr(1:msz)*r2k(1:msz) 2161 call simp_gen(intvh,work1,radmesh(imainmesh)) 2162 do jlmn=1,pawtab%lmn_size 2163 j0lmn=jlmn*(jlmn-1)/2;jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn) 2164 do ilmn=1,jlmn 2165 klmn=j0lmn+ilmn;ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn) 2166 if (jlm==ilm) then 2167 work1(1:msz)=pawtab%tphi(1:msz,iln)*pawtab%tphi(1:msz,jln)*(vlocr(1:msz)-intvh) & 2168 & -pawtab%phi (1:msz,iln)*pawtab%phi (1:msz,jln)*intvh 2169 call simp_gen(intg,work1,radmesh(imainmesh)) 2170 pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg 2171 end if 2172 end do 2173 end do 2174 LIBPAW_DEALLOCATE(work1) 2175 end if 2176 LIBPAW_DEALLOCATE(r2k) 2177 end if 2178 2179 !========================================================== 2180 !If usexcnhat in psp file is different from usexcnhat chosen 2181 !by user, convert VH(tnzc) and Dij0 2182 2183 if (pawtab%usexcnhat==-1) then 2184 pawtab%usexcnhat=usexcnhat 2185 else if (usexcnhat/=pawtab%usexcnhat) then 2186 if (pawtab%has_tvale==0) then 2187 write(msg, '(5a)' ) & 2188 & 'It is only possible to modify the use of compensation charge density',ch10,& 2189 & 'for a file format containing the pseudo valence density (format>=paw4 or XML)!',ch10,& 2190 & 'Action: use usexcnhat=-1 in input file or change psp file format.' 2191 LIBPAW_ERROR(msg) 2192 else if (usekden>=1) then 2193 write(msg, '(5a)' ) & 2194 & 'It is not possible to modify the use of compensation charge density',ch10,& 2195 & 'within the metaGGA XC functional (need valence kinetic density)!',ch10,& 2196 & 'Action: use usexcnhat=-1 in input file or change psp file format.' 2197 LIBPAW_ERROR(msg) 2198 else 2199 msz=vloc_mesh%mesh_size 2200 ! Retrieve tvale and nhat onto vloc mesh 2201 LIBPAW_ALLOCATE(nwk,(msz)) 2202 LIBPAW_ALLOCATE(ncorwk,(msz)) 2203 LIBPAW_ALLOCATE(nhatwk,(msz)) 2204 nwk=zero;ncorwk=zero;nhatwk=zero 2205 if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2206 & (core_mesh%rstep /=vloc_mesh%rstep) .or.& 2207 & (core_mesh%lstep /=vloc_mesh%lstep)) then 2208 LIBPAW_ALLOCATE(work1,(core_mesh%mesh_size)) 2209 msz1=msz;if (core_mesh%rmax<vloc_mesh%rmax) msz1=pawrad_ifromr(vloc_mesh,core_mesh%rmax) 2210 call bound_deriv(tncore(1:core_mesh%mesh_size),core_mesh,core_mesh%mesh_size,yp1,ypn) 2211 call paw_spline(core_mesh%rad,tncore,core_mesh%mesh_size,yp1,ypn,work1) 2212 call paw_splint(core_mesh%mesh_size,core_mesh%rad,tncore,work1,msz1,vloc_mesh%rad(1:msz1),ncorwk(1:msz1)) 2213 LIBPAW_DEALLOCATE(work1) 2214 else 2215 msz1=min(core_mesh%mesh_size,msz) 2216 ncorwk(1:msz1)=tncore(1:msz1) 2217 end if 2218 if ((vale_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2219 & (vale_mesh%rstep /=vloc_mesh%rstep) .or.& 2220 & (vale_mesh%lstep /=vloc_mesh%lstep)) then 2221 LIBPAW_ALLOCATE(work1,(vale_mesh%mesh_size)) 2222 msz1=msz;if (vale_mesh%rmax<vloc_mesh%rmax) msz1=pawrad_ifromr(vloc_mesh,vale_mesh%rmax) 2223 call bound_deriv(tnvale(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn) 2224 call paw_spline(vale_mesh%rad,tnvale,vale_mesh%mesh_size,yp1,ypn,work1) 2225 call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,tnvale,work1,msz1,vloc_mesh%rad(1:msz1),nwk(1:msz1)) 2226 call bound_deriv(nhat(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn) 2227 call paw_spline(vale_mesh%rad,nhat,vale_mesh%mesh_size,yp1,ypn,work1) 2228 call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,nhat,work1,msz1,vloc_mesh%rad(1:msz1),nhatwk(1:msz1)) 2229 LIBPAW_DEALLOCATE(work1) 2230 else 2231 msz1=min(vale_mesh%mesh_size,msz) 2232 nwk (1:msz1)=tnvale(1:msz1) 2233 nhatwk(1:msz1)=nhat (1:msz1) 2234 end if 2235 2236 nwk=nwk-nhatwk 2237 nwk=sqrt(four_pi)*nwk;nhatwk=sqrt(four_pi)*nhatwk ! 0th-order moment of densities 2238 2239 ! Compute Vxc without nhat (vxc1) and with nhat (vxc2) 2240 nspden=1 2241 #if defined LIBPAW_HAVE_LIBXC 2242 if (ixc<0) nspden=libxc_functionals_nspin() 2243 #endif 2244 if (ixc<0) then 2245 LIBPAW_ALLOCATE(vxc1,(msz*nspden)) 2246 LIBPAW_ALLOCATE(vxc2,(msz*nspden)) 2247 LIBPAW_ALLOCATE(work1,(msz)) 2248 LIBPAW_ALLOCATE(work2,(msz*nspden)) 2249 LIBPAW_ALLOCATE(work3,(msz*nspden)) 2250 tmp1 => work1 2251 work2(1:msz)=nwk 2252 work3(1:msz)=nhatwk 2253 if (nspden==2) work2(msz+1:2*msz)=half*nwk 2254 if (nspden==2) work3(msz+1:2*msz)=half*nhatwk 2255 if (pawxcdev/=0) then 2256 call pawxcm(ncorwk,yp1,ypn,0,hyb_mixing,ixc,work1,1,tmp_lmselect,work3,0,non_magnetic_xc,msz,nspden,5,& 2257 & pawang_tmp,vloc_mesh,pawxcdev,work2,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2258 call pawxcm(ncorwk,yp1,ypn,0,hyb_mixing,ixc,work1,1,tmp_lmselect,work3,0,non_magnetic_xc,msz,nspden,5,& 2259 & pawang_tmp,vloc_mesh,pawxcdev,work2,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2260 vxc1=vxc1/sqrt(four_pi);vxc2=vxc2/sqrt(four_pi) ! Deduce Vxc from its first moment 2261 else 2262 call pawxc(ncorwk,yp1,ypn,hyb_mixing,ixc,work1,tmp1,1,tmp_lmselect,work3,0,0,non_magnetic_xc,msz,nspden,5,& 2263 & pawang_tmp,vloc_mesh,work2,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2264 call pawxc(ncorwk,yp1,ypn,hyb_mixing,ixc,work1,tmp1,1,tmp_lmselect,work3,0,0,non_magnetic_xc,msz,nspden,5,& 2265 & pawang_tmp,vloc_mesh,work2,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2266 end if 2267 LIBPAW_DEALLOCATE(nwk) 2268 LIBPAW_DEALLOCATE(ncorwk) 2269 LIBPAW_DEALLOCATE(nhatwk) 2270 LIBPAW_DEALLOCATE(work1) 2271 LIBPAW_DEALLOCATE(work2) 2272 LIBPAW_DEALLOCATE(work3) 2273 else 2274 LIBPAW_ALLOCATE(vxc1,(msz)) 2275 LIBPAW_ALLOCATE(vxc2,(msz)) 2276 LIBPAW_ALLOCATE(work1,(msz)) 2277 tmp1 => work1 2278 if (pawxcdev/=0) then 2279 call pawxcm(ncorwk,yp1,ypn,0,hyb_mixing,ixc,work1,1,tmp_lmselect,nhatwk,0,non_magnetic_xc,msz,1,5,& 2280 & pawang_tmp,vloc_mesh,pawxcdev,nwk,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2281 call pawxcm(ncorwk,yp1,ypn,0,hyb_mixing,ixc,work1,1,tmp_lmselect,nhatwk,0,non_magnetic_xc,msz,1,5,& 2282 & pawang_tmp,vloc_mesh,pawxcdev,nwk,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2283 vxc1=vxc1/sqrt(four_pi);vxc2=vxc2/sqrt(four_pi) ! Deduce Vxc from its first moment 2284 else 2285 call pawxc(ncorwk,yp1,ypn,hyb_mixing,ixc,work1,tmp1,1,tmp_lmselect,nhatwk,0,0,non_magnetic_xc,msz,1,5,& 2286 & pawang_tmp,vloc_mesh,nwk,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2287 call pawxc(ncorwk,yp1,ypn,hyb_mixing,ixc,work1,tmp1,1,tmp_lmselect,nhatwk,0,0,non_magnetic_xc,msz,1,5,& 2288 & pawang_tmp,vloc_mesh,nwk,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2289 end if 2290 LIBPAW_DEALLOCATE(nwk) 2291 LIBPAW_DEALLOCATE(ncorwk) 2292 LIBPAW_DEALLOCATE(nhatwk) 2293 LIBPAW_DEALLOCATE(work1) 2294 endif 2295 ! Compute difference of XC potentials 2296 if (usexcnhat==0.and.pawtab%usexcnhat/=0) vxc1(1:msz)=vxc2(1:msz)-vxc1(1:msz) 2297 if (usexcnhat/=0.and.pawtab%usexcnhat==0) vxc1(1:msz)=vxc1(1:msz)-vxc2(1:msz) 2298 ! Modify VH(tnzc) 2299 vlocr(1:msz)=vlocr(1:msz)-vxc1(1:msz) 2300 if (has_dij0) then 2301 ! Modify Dij0 2302 LIBPAW_ALLOCATE(work2,(pawtab%lmn2_size)) 2303 call atompaw_kij(pawtab%indlmn,work2,pawtab%lmn_size,ncore,0,0,pawtab,pawrad,& 2304 & core_mesh,vloc_mesh,vxc1(1:msz),znucl) 2305 pawtab%dij0=work2 2306 LIBPAW_DEALLOCATE(work2) 2307 end if 2308 LIBPAW_DEALLOCATE(vxc1) 2309 LIBPAW_DEALLOCATE(vxc2) 2310 end if ! has_tvale/=0 2311 end if 2312 if (pawtab%usexcnhat==0) then 2313 write(msg,'(a)') & 2314 & ' Compensation charge density is not taken into account in XC energy/potential' 2315 call wrtout(ab_out,msg,'COLL') 2316 call wrtout(std_out, msg,'COLL') 2317 end if 2318 if (pawtab%usexcnhat==1) then 2319 write(msg,'(a)') & 2320 & ' Compensation charge density is taken into account in XC energy/potential' 2321 call wrtout(ab_out,msg,'COLL') 2322 call wrtout(std_out, msg,'COLL') 2323 end if 2324 2325 !========================================================== 2326 ! Calculate the coefficient beta = \int { vH[nZc](r) - vloc(r) } 4pi r^2 dr 2327 ! 2328 LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size)) 2329 LIBPAW_ALLOCATE(nwk,(core_mesh%mesh_size)) 2330 ! get vH[nZc] 2331 call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl) 2332 2333 !Transpose vlocr mesh into core mesh 2334 nwk(:)=zero 2335 if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2336 & (core_mesh%rstep /=vloc_mesh%rstep) .or.& 2337 & (core_mesh%lstep /=vloc_mesh%lstep)) then 2338 msz=core_mesh%mesh_size;if (vloc_mesh%rmax<core_mesh%rmax) msz=pawrad_ifromr(core_mesh,vloc_mesh%rmax) 2339 call bound_deriv(vlocr(1:vloc_mesh%mesh_size),vloc_mesh,vloc_mesh%mesh_size,yp1,ypn) 2340 LIBPAW_ALLOCATE(work1,(vloc_mesh%mesh_size)) 2341 LIBPAW_ALLOCATE(work2,(vloc_mesh%mesh_size)) 2342 call paw_spline(vloc_mesh%rad,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1) 2343 call paw_splint(vloc_mesh%mesh_size,vloc_mesh%rad,vlocr,work1,msz,core_mesh%rad(1:msz),nwk) 2344 LIBPAW_DEALLOCATE(work1) 2345 LIBPAW_DEALLOCATE(work2) 2346 else 2347 msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size) 2348 nwk(1:msz)=vlocr(1:msz) 2349 end if 2350 2351 !Difference 2352 nwk(1:msz)=vhnzc(1:msz)-nwk(1:msz) 2353 if (msz<core_mesh%mesh_size) nwk(msz+1:core_mesh%mesh_size)=zero 2354 2355 !Perform the spherical integration 2356 nwk(1:msz)=nwk(1:msz)*four_pi*core_mesh%rad(1:msz)**2 2357 2358 call simp_gen(pawtab%beta,nwk,core_mesh) 2359 2360 LIBPAW_DEALLOCATE(vhnzc) 2361 LIBPAW_DEALLOCATE(nwk) 2362 2363 write(msg,'(a,e18.6)') & 2364 & ' beta integral value: ',pawtab%beta 2365 call wrtout(std_out,msg,'COLL') 2366 2367 2368 !========================================================== 2369 !Try to optimize CPU time: 2370 !If Vloc mesh size is big, spline Vloc into a smaller log. mesh 2371 2372 reduced_vloc=(vloc_mesh%mesh_size>int(reduced_mshsz)) 2373 if (reduced_vloc) then 2374 msz=vloc_mesh%mesh_size 2375 lstep_tmp=log(0.9999999_dp*vloc_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2376 call pawrad_init(rvloc_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2377 & rstep=reduced_rstep,lstep=lstep_tmp) 2378 LIBPAW_ALLOCATE(rvlocr,(reduced_mshsz)) 2379 call bound_deriv(vlocr(1:msz),vloc_mesh,msz,yp1,ypn) 2380 LIBPAW_ALLOCATE(work1,(msz)) 2381 LIBPAW_ALLOCATE(work2,(msz)) 2382 LIBPAW_ALLOCATE(work3,(msz)) 2383 work3(1:msz)=vloc_mesh%rad(1:msz) 2384 call paw_spline(work3,vlocr,msz,yp1,ypn,work1) 2385 call paw_splint(msz,work3,vlocr,work1,reduced_mshsz,rvloc_mesh%rad,rvlocr) 2386 LIBPAW_DEALLOCATE(work1) 2387 LIBPAW_DEALLOCATE(work2) 2388 LIBPAW_DEALLOCATE(work3) 2389 end if 2390 2391 !Keep VH(tnZc) eventually in memory 2392 if (pawtab%has_vhtnzc==1) then 2393 LIBPAW_ALLOCATE(pawtab%vhtnzc,(pawtab%mesh_size)) 2394 if ((reduced_vloc).and.(rvloc_mesh%mesh_type==pawrad%mesh_type)& 2395 & .and.(rvloc_mesh%rstep==pawrad%rstep).and.(rvloc_mesh%lstep==pawrad%lstep)) then 2396 pawtab%vhtnzc(1:pawtab%mesh_size)=rvlocr(1:pawtab%mesh_size) 2397 pawtab%has_vhtnzc=2 2398 else if ((vloc_mesh%mesh_type==pawrad%mesh_type)& 2399 & .and.(vloc_mesh%rstep==pawrad%rstep).and.(vloc_mesh%lstep==pawrad%lstep)) then 2400 pawtab%vhtnzc(1:pawtab%mesh_size)=vlocr(1:pawtab%mesh_size) 2401 pawtab%has_vhtnzc=2 2402 else 2403 msg = 'Vloc mesh is not right !' 2404 LIBPAW_ERROR(msg) 2405 end if 2406 end if 2407 2408 !========================================================== 2409 !Try to optimize CPU time: 2410 !If ncore mesh size is big, spline tncore into a smaller log. mesh 2411 2412 reduced_ncor=(core_mesh%mesh_size>int(reduced_mshsz)).and.(pawtab%usetcore/=0) 2413 if (reduced_ncor) then 2414 msz=core_mesh%mesh_size 2415 lstep_tmp=log(0.9999999_dp*core_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2416 call pawrad_init(rcore_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2417 & rstep=reduced_rstep,lstep=lstep_tmp) 2418 LIBPAW_ALLOCATE(rtncor,(reduced_mshsz)) 2419 call bound_deriv(tncore(1:msz),core_mesh,msz,yp1,ypn) 2420 LIBPAW_ALLOCATE(work1,(msz)) 2421 LIBPAW_ALLOCATE(work2,(msz)) 2422 LIBPAW_ALLOCATE(work3,(msz)) 2423 work3(1:msz)=core_mesh%rad(1:msz) 2424 call paw_spline(work3,tncore,msz,yp1,ypn,work1) 2425 call paw_splint(msz,work3,tncore,work1,reduced_mshsz,rcore_mesh%rad,rtncor) 2426 LIBPAW_DEALLOCATE(work1) 2427 LIBPAW_DEALLOCATE(work2) 2428 LIBPAW_DEALLOCATE(work3) 2429 end if 2430 2431 !========================================================== 2432 !Try to optimize CPU time: 2433 !If coretau mesh size is big, spline tcoretau into a smaller log. mesh 2434 2435 reduced_taucor=.false. 2436 if (usekden>=1.and.present(tcoretau)) then 2437 reduced_taucor=(coretau_mesh%mesh_size>int(reduced_mshsz)).and.(pawtab%usetcore/=0) 2438 if (reduced_taucor) then 2439 msz=coretau_mesh%mesh_size 2440 lstep_tmp=log(0.9999999_dp*coretau_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2441 call pawrad_init(rcoretau_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2442 & rstep=reduced_rstep,lstep=lstep_tmp) 2443 LIBPAW_ALLOCATE(rttaucor,(reduced_mshsz)) 2444 call bound_deriv(tcoretau(1:msz),coretau_mesh,msz,yp1,ypn) 2445 LIBPAW_ALLOCATE(work1,(msz)) 2446 LIBPAW_ALLOCATE(work2,(msz)) 2447 LIBPAW_ALLOCATE(work3,(msz)) 2448 work3(1:msz)=coretau_mesh%rad(1:msz) 2449 call paw_spline(work3,tcoretau,msz,yp1,ypn,work1) 2450 call paw_splint(msz,work3,tcoretau,work1,reduced_mshsz,rcoretau_mesh%rad,rttaucor) 2451 LIBPAW_DEALLOCATE(work1) 2452 LIBPAW_DEALLOCATE(work2) 2453 LIBPAW_DEALLOCATE(work3) 2454 end if 2455 end if 2456 2457 !========================================================== 2458 !Try to optimize CPU time: 2459 !If vale mesh size is big, spline tnvale into a smaller log. mesh 2460 2461 if (pawtab%has_tvale==1) then 2462 reduced_nval=(vale_mesh%mesh_size>int(reduced_mshsz)) 2463 if (reduced_nval) then 2464 msz=vale_mesh%mesh_size 2465 lstep_tmp=log(0.9999999_dp*vale_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2466 call pawrad_init(rvale_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2467 & rstep=reduced_rstep,lstep=lstep_tmp) 2468 LIBPAW_ALLOCATE(rtnval,(reduced_mshsz)) 2469 call bound_deriv(tnvale(1:msz),vale_mesh,msz,yp1,ypn) 2470 LIBPAW_ALLOCATE(work1,(msz)) 2471 LIBPAW_ALLOCATE(work2,(msz)) 2472 LIBPAW_ALLOCATE(work3,(msz)) 2473 work3(1:msz)=vale_mesh%rad(1:msz) 2474 call paw_spline(work3,tnvale,msz,yp1,ypn,work1) 2475 call paw_splint(msz,work3,tnvale,work1,reduced_mshsz,rvale_mesh%rad,rtnval) 2476 LIBPAW_DEALLOCATE(work1) 2477 LIBPAW_DEALLOCATE(work2) 2478 LIBPAW_DEALLOCATE(work3) 2479 end if 2480 else 2481 reduced_nval=.false. 2482 end if 2483 !========================================================== 2484 !Compute Vlspl(q) (and second derivative) from Vloc(r) 2485 2486 !Compute Vlspl(q)=q^2.Vloc(q) from vloc(r) 2487 if(mqgrid_vl>0) then 2488 if (reduced_vloc) then 2489 call pawpsp_lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),rvloc_mesh,rvlocr,yp1,ypn,zion) 2490 else 2491 call pawpsp_lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),vloc_mesh,vlocr,yp1,ypn,zion) 2492 end if 2493 ! Compute second derivative of Vlspl(q) 2494 call paw_spline(qgrid_vl,vlspl(:,1),mqgrid_vl,yp1,ypn,vlspl(:,2)) 2495 else 2496 ! Only to compute epsatm 2497 epsatm=zero 2498 if (reduced_vloc) then 2499 call pawpsp_lo(epsatm,1,tmp_qgrid,tmp_q2vq,rvloc_mesh,rvlocr,yp1,ypn,zion) 2500 else 2501 call pawpsp_lo(epsatm,1,tmp_qgrid,tmp_q2vq,vloc_mesh,vlocr,yp1,ypn,zion) 2502 end if 2503 end if 2504 !========================================================== 2505 !Compute tcorespl(q) (and second derivative) from tNcore(r) 2506 2507 pawtab%mqgrid=mqgrid_vl 2508 xcccrc=core_mesh%rmax 2509 LIBPAW_ALLOCATE(pawtab%tcorespl,(pawtab%mqgrid,2)) 2510 2511 if(mqgrid_vl>0.and.pawtab%usetcore/=0) then 2512 ! Compute tcorespl(q)=tNc(q) from tNcore(r) 2513 if (reduced_ncor) then 2514 call pawpsp_cg(pawtab%dncdq0,pawtab%d2ncdq0,mqgrid_vl,qgrid_vl,pawtab%tcorespl(:,1),rcore_mesh,rtncor,yp1,ypn) 2515 else 2516 call pawpsp_cg(pawtab%dncdq0,pawtab%d2ncdq0,mqgrid_vl,qgrid_vl,pawtab%tcorespl(:,1),core_mesh,tncore,yp1,ypn) 2517 end if 2518 ! Compute second derivative of tcorespl(q) 2519 call paw_spline(qgrid_vl,pawtab%tcorespl(:,1),mqgrid_vl,yp1,ypn,pawtab%tcorespl(:,2)) 2520 else 2521 pawtab%tcorespl=zero 2522 pawtab%dncdq0=zero 2523 pawtab%d2ncdq0=zero 2524 end if 2525 2526 !========================================================== 2527 !Compute tcoretauspl(q) 2528 2529 if (present(tcoretau)) then 2530 LIBPAW_ALLOCATE(pawtab%tcoretauspl,(pawtab%mqgrid,2*usekden)) 2531 if (usekden==1) then 2532 if (coretau_mesh%rmax/=xcccrc) then 2533 write(msg, '(a,a,a)' )& 2534 & 'Core density and core kinetic density should be given on the same grid!',ch10,& 2535 & 'Action : check your pseudopotential (increase tNvale meshSize).' 2536 LIBPAW_ERROR(msg) 2537 end if 2538 if(mqgrid_vl>0) then 2539 ! Compute tcorespl(q)=tNc(q) from tNcore(r) 2540 if (reduced_taucor) then 2541 call pawpsp_cg(pawtab%dtaucdq0,qq,mqgrid_vl,qgrid_vl,pawtab%tcoretauspl(:,1),rcoretau_mesh,rttaucor,yp1,ypn) 2542 else 2543 call pawpsp_cg(pawtab%dtaucdq0,qq,mqgrid_vl,qgrid_vl,pawtab%tcoretauspl(:,1),coretau_mesh,tcoretau,yp1,ypn) 2544 end if 2545 ! Compute second derivative of tcorespl(q) 2546 call paw_spline(qgrid_vl,pawtab%tcoretauspl(:,1),mqgrid_vl,yp1,ypn,pawtab%tcoretauspl(:,2)) 2547 else 2548 pawtab%tcoretauspl=zero 2549 pawtab%dtaucdq0=zero 2550 end if 2551 end if 2552 end if 2553 2554 !========================================================== 2555 !Compute tvalespl(q) (and second derivative) from tNvale(r) 2556 2557 if (pawtab%has_tvale/=0.and.mqgrid_vl>0) then 2558 LIBPAW_ALLOCATE(pawtab%tvalespl,(pawtab%mqgrid,2)) 2559 if (reduced_nval) then 2560 call pawpsp_cg(pawtab%dnvdq0,d2nvdq0,mqgrid_vl,qgrid_vl,pawtab%tvalespl(:,1),rvale_mesh,rtnval,yp1,ypn) 2561 pawtab%tnvale_mesh_size=rvale_mesh%mesh_size 2562 else 2563 call pawpsp_cg(pawtab%dnvdq0,d2nvdq0,mqgrid_vl,qgrid_vl,pawtab%tvalespl(:,1),vale_mesh,tnvale,yp1,ypn) 2564 pawtab%tnvale_mesh_size=vale_mesh%mesh_size 2565 end if 2566 ! Compute second derivative of tvalespl(q) 2567 call paw_spline(qgrid_vl,pawtab%tvalespl(:,1),mqgrid_vl,yp1,ypn,pawtab%tvalespl(:,2)) 2568 else 2569 pawtab%dnvdq0=zero 2570 pawtab%tnvale_mesh_size=0 2571 end if 2572 2573 !================================================== 2574 !Compute Ex-correlation energy for the core density 2575 2576 nspden=1 2577 #if defined LIBPAW_HAVE_LIBXC 2578 if (ixc<0) nspden=libxc_functionals_nspin() 2579 #endif 2580 2581 LIBPAW_ALLOCATE(work1,(core_mesh%mesh_size*nspden)) 2582 LIBPAW_ALLOCATE(work2,(core_mesh%mesh_size)) 2583 LIBPAW_ALLOCATE(work3,(1)) 2584 work1(:)=zero;work2(:)=zero;work3(:)=zero 2585 tmp1 => work1 ; tmp2 => work1 2586 2587 if (pawxcdev/=0) then 2588 call pawxcm(ncore,pawtab%exccore,yp1,0,hyb_mixing,ixc,work2,1,tmp_lmselect,work3,0,non_magnetic_xc,core_mesh%mesh_size,& 2589 & nspden,4,pawang_tmp,core_mesh,pawxcdev,work1,1,0,tmp1,xclevel,xc_denpos) 2590 else 2591 if (present(tcoretau)) then 2592 call pawxc(ncore,pawtab%exccore,yp1,hyb_mixing,ixc,work2,work1,1,tmp_lmselect,work3,0,0,non_magnetic_xc,core_mesh%mesh_size,& 2593 & nspden,4,pawang_tmp,core_mesh,tmp1,1,0,tmp2,xclevel,xc_denpos,coretau=tcoretau,xc_taupos=my_xc_taupos) 2594 else 2595 call pawxc(ncore,pawtab%exccore,yp1,hyb_mixing,ixc,work2,work1,1,tmp_lmselect,work3,0,0,non_magnetic_xc,core_mesh%mesh_size,& 2596 & nspden,4,pawang_tmp,core_mesh,tmp1,1,0,tmp2,xclevel,xc_denpos) 2597 end if 2598 end if 2599 2600 LIBPAW_DEALLOCATE(work1) 2601 LIBPAW_DEALLOCATE(work2) 2602 LIBPAW_DEALLOCATE(work3) 2603 2604 !================================================== 2605 !Compute atomic contribution to Dij (Dij0) 2606 !if not already in memory 2607 if ((.not.has_dij0).and.(pawtab%has_kij==2.or.pawtab%has_kij==-1)) then 2608 LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size)) 2609 if (reduced_vloc) then 2610 call atompaw_dij0(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,& 2611 & rvloc_mesh,rvlocr,znucl) 2612 else 2613 call atompaw_dij0(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,& 2614 & vloc_mesh,vlocr,znucl) 2615 end if 2616 has_dij0=.true. 2617 end if 2618 !================================================== 2619 !Compute kinetic operator contribution to Dij 2620 2621 if (pawtab%has_kij==1.and.has_dij0) then 2622 LIBPAW_ALLOCATE(pawtab%kij,(pawtab%lmn2_size)) 2623 call atompaw_kij(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,1,pawtab,pawrad,core_mesh,& 2624 & vloc_mesh,vlocr,znucl) 2625 pawtab%has_kij=2 2626 end if 2627 2628 !pawtab%has_kij=-1 means that kij does not have to be kept in memory 2629 if (pawtab%has_kij==-1) then 2630 LIBPAW_DEALLOCATE(pawtab%kij) 2631 pawtab%has_kij=0 2632 end if 2633 2634 !========================================================== 2635 !If projectors have to be kept in memory, we need 2636 !them on the main radial mesh (so, spline them if necessary) 2637 2638 if (pawtab%has_tproj>0) then 2639 if ((tproj_mesh%mesh_type/=pawrad%mesh_type).or.& 2640 & (tproj_mesh%rstep /=pawrad%rstep).or.& 2641 & (tproj_mesh%lstep /=pawrad%lstep)) then 2642 ir=pawrad_ifromr(pawrad,tproj_mesh%rmax) 2643 call pawrad_init(tproj_mesh_new,mesh_size=ir,mesh_type=pawrad%mesh_type,& 2644 & rstep=pawrad%rstep,lstep=pawrad%lstep) 2645 LIBPAW_ALLOCATE(pawtab%tproj,(tproj_mesh_new%mesh_size,pawtab%basis_size)) 2646 LIBPAW_ALLOCATE(work1,(tproj_mesh%mesh_size)) 2647 do ib=1,pawtab%basis_size 2648 call bound_deriv(tproj(:,ib),tproj_mesh,tproj_mesh%mesh_size,yp1,ypn) 2649 call paw_spline(tproj_mesh%rad,tproj(:,ib),tproj_mesh%mesh_size,yp1,ypn,work1) 2650 call paw_splint(tproj_mesh%mesh_size,tproj_mesh%rad,tproj(:,ib),work1,& 2651 & tproj_mesh_new%mesh_size,tproj_mesh_new%rad,pawtab%tproj(:,ib)) 2652 end do 2653 LIBPAW_DEALLOCATE(work1) 2654 call pawrad_free(tproj_mesh_new) 2655 else 2656 LIBPAW_ALLOCATE(pawtab%tproj,(tproj_mesh%mesh_size,pawtab%basis_size)) 2657 pawtab%tproj(:,:)=tproj(:,:) 2658 end if 2659 pawtab%has_tproj=2 2660 end if 2661 2662 !========================================================== 2663 !Free temporary allocated space 2664 2665 if (pawtab%has_tvale==1) then 2666 LIBPAW_DEALLOCATE(nhat) 2667 end if 2668 if (reduced_vloc) then 2669 call pawrad_free(rvloc_mesh) 2670 LIBPAW_DEALLOCATE(rvlocr) 2671 end if 2672 if (reduced_ncor) then 2673 call pawrad_free(rcore_mesh) 2674 LIBPAW_DEALLOCATE(rtncor) 2675 end if 2676 if (reduced_taucor) then 2677 call pawrad_free(rcoretau_mesh) 2678 LIBPAW_DEALLOCATE(rttaucor) 2679 end if 2680 if (reduced_nval) then 2681 call pawrad_free(rvale_mesh) 2682 LIBPAW_DEALLOCATE(rtnval) 2683 end if 2684 if (pawxcdev==0) then 2685 LIBPAW_DEALLOCATE(pawang_tmp%angwgth) 2686 LIBPAW_DEALLOCATE(pawang_tmp%anginit) 2687 LIBPAW_DEALLOCATE(pawang_tmp%ylmr) 2688 LIBPAW_DEALLOCATE(pawang_tmp%ylmrgr) 2689 end if 2690 2691 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
SOURCE
2714 subroutine pawpsp_calc_d5(mesh,mesh_size,tcoredens) 2715 2716 !Arguments ------------------------------------ 2717 integer,intent(in) :: mesh_size 2718 type(pawrad_type),intent(in) :: mesh 2719 real(dp),intent(inout) :: tcoredens(mesh_size,6) 2720 2721 !Local variables------------------------------- 2722 integer,parameter :: it=1 !number of steps for smoothing function 2723 logical,parameter :: use_smooth=.true. 2724 2725 ! ************************************************************************* 2726 2727 !calculate first derivative from density, 2728 !and store it 2729 call nderiv_gen(tcoredens(:,2),tcoredens(:,1),mesh) 2730 2731 !get second derivative from density, and store it 2732 call paw_spline(mesh%rad,tcoredens(:,1),mesh_size,& 2733 & zero,zero,tcoredens(:,3)) 2734 2735 !smooth functions, to avoid numerical instabilities 2736 if(use_smooth) then 2737 call paw_smooth(tcoredens(:,2),mesh_size,it) 2738 call paw_smooth(tcoredens(:,3),mesh_size,it) 2739 end if 2740 2741 !get third derivative from first derivative: 2742 call paw_spline(mesh%rad,tcoredens(:,2),mesh_size,& 2743 & zero,zero,tcoredens(:,4)) 2744 2745 !get fourth derivative from second derivative: 2746 call paw_spline(mesh%rad,tcoredens(:,3),mesh_size,& 2747 & zero,zero,tcoredens(:,5)) 2748 2749 !smooth 3rd and 4th order derivatives 2750 if(use_smooth) then 2751 call paw_smooth(tcoredens(:,4),mesh_size,it) 2752 call paw_smooth(tcoredens(:,5),mesh_size,it) 2753 end if 2754 2755 !get fifth derivative from third derivative: 2756 call paw_spline(mesh%rad,tcoredens(:,4),mesh_size,& 2757 & zero,zero,tcoredens(:,6)) 2758 2759 !smooth 5th order derivative 2760 if(use_smooth) then 2761 call paw_smooth(tcoredens(:,6),mesh_size,it) 2762 end if 2763 2764 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).
SOURCE
458 subroutine pawpsp_cg(dnqdq0,d2nqdq0,mqgrid,qgrid,nq,radmesh,nr,yp1,ypn) 459 460 !Arguments---------------------------------------------------------- 461 !scalars 462 integer,intent(in) :: mqgrid 463 real(dp),intent(out) :: dnqdq0,d2nqdq0,yp1,ypn 464 type(pawrad_type),intent(in) :: radmesh 465 !arrays 466 real(dp),intent(in) :: nr(:) 467 real(dp),intent(in) :: qgrid(mqgrid) 468 real(dp),intent(out) :: nq(mqgrid) 469 470 !Local variables------------------------------- 471 !scalars 472 integer :: iq,ir,mesh_size 473 real(dp) :: aexp,arg,bexp,dn,r0tor1,r1torm,rm,rmtoin 474 logical :: begin_r0 475 !character(len=500) :: msg 476 !arrays 477 real(dp),allocatable :: ff(:),rnr(:) 478 479 ! ************************************************************************* 480 481 mesh_size=min(size(nr),radmesh%mesh_size) 482 LIBPAW_ALLOCATE(ff,(mesh_size)) 483 LIBPAW_ALLOCATE(rnr,(mesh_size)) 484 ff=zero;rnr=zero 485 486 do ir=1,mesh_size 487 rnr(ir)=radmesh%rad(ir)*nr(ir) 488 end do 489 490 !Is mesh beginning with r=0 ? 491 begin_r0=(radmesh%rad(1)<1.d-20) 492 493 !Adjustment of an exponentional at r_max (n_exp(r)=aexp*Exp[-bexp*r]) 494 rm=radmesh%rad(mesh_size) 495 dn=one/(12._dp*radmesh%stepint*radmesh%radfact(mesh_size)) & 496 & *( 3._dp*nr(mesh_size-4) & 497 & -16._dp*nr(mesh_size-3) & 498 & +36._dp*nr(mesh_size-2) & 499 & -48._dp*nr(mesh_size-1) & 500 & +25._dp*nr(mesh_size)) 501 if (dn<0._dp.and. & 502 & abs(radmesh%rad(mesh_size)*nr(mesh_size))>1.d-20) then 503 bexp=-dn/nr(mesh_size) 504 if (bexp * rm > 50._dp) then 505 ! This solves the problem with the weird core charge used in v4[62] in which bexp x rm ~= 10^3 506 !write(msg,"(a,es16.8)")"Tooooo large bexp * rm: ", bexp*rm, ", setting aexp to 0" 507 !LIBPAW_WARNING(msg) 508 bexp=0.001_dp;aexp=zero 509 else 510 aexp=nr(mesh_size)*exp(bexp*rm) 511 if (abs(aexp)<1.d-20) then 512 bexp=0.001_dp;aexp=zero 513 end if 514 end if 515 else 516 bexp=0.001_dp;aexp=zero 517 end if 518 519 !=========================================== 520 !=== Compute n(q) for q=0 separately 521 !=========================================== 522 523 !Integral from 0 to r1 (only if r1<>0) 524 r0tor1=zero 525 if (.not.begin_r0) r0tor1=(rnr(1)*radmesh%rad(1)**2)/3.d0 526 527 !Integral from r1 to rmax 528 do ir=1,mesh_size 529 if (abs(rnr(ir))>1.d-20) ff(ir)=rnr(ir)*radmesh%rad(ir) 530 end do 531 call simp_gen(r1torm,ff,radmesh) 532 533 !Integral from rmax to infinity 534 !This part is approximated using an exponential density aexp*Exp[-bexp*r] 535 !(formulae obtained with mathematica) 536 rmtoin=aexp*exp(-bexp*rm)/bexp**3*(two+two*bexp*rm+bexp*bexp*rm*rm) 537 538 !Some of the three parts 539 nq(1)=four_pi*(r0tor1+r1torm+rmtoin) 540 541 !=========================================== 542 !=== Compute n(q) for other q''s 543 !=========================================== 544 545 !Loop over q values 546 do iq=2,mqgrid 547 arg=two_pi*qgrid(iq) 548 549 ! Integral from 0 to r1 (only if r1<>0) 550 r0tor1=zero;if (.not.begin_r0) & 551 & r0tor1=nr(1)*(sin(arg*radmesh%rad(1))/arg/arg& 552 & -radmesh%rad(1)*cos(arg*radmesh%rad(1))/arg) 553 554 ! Integral from r1 to rmax 555 do ir=1,mesh_size 556 if (abs(rnr(ir))>1.d-20) ff(ir)=sin(arg*radmesh%rad(ir))*rnr(ir) 557 end do 558 call simp_gen(r1torm,ff,radmesh) 559 560 ! Integral from rmax to infinity 561 ! This part is approximated using an exponential density aexp*Exp[-bexp*r] 562 ! (formulae obtained with mathematica) 563 rmtoin=aexp*exp(-bexp*rm)/(arg**2+bexp**2)**2 & 564 & *(arg*(two*bexp+arg**2*rm+bexp**2*rm)*cos(arg*rm) & 565 & +(arg**2*(bexp*rm-one)+bexp**2*(bexp*rm+one))*sin(arg*rm)) 566 567 ! Store q^2 v(q) 568 nq(iq)=two/qgrid(iq)*(r0tor1+r1torm+rmtoin) 569 end do 570 571 !=========================================== 572 !=== Compute derivatives of n(q) 573 !=== at ends of interval 574 !=========================================== 575 576 !yp(0)=zero 577 yp1=zero 578 579 !yp(qmax)=$ 2\int_0^\infty[(-\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r) r n(r) dr]$ 580 arg=two_pi*qgrid(mqgrid) 581 582 !Integral from 0 to r1 (only if r1<>0) 583 r0tor1=zero;if (.not.begin_r0) & 584 & r0tor1=two_pi*nr(1)*(3.d0*radmesh%rad(1)/arg /arg*cos(arg*radmesh%rad(1))+ & 585 & (radmesh%rad(1)**2/arg-3.0d0/arg**3)*sin(arg*radmesh%rad(1))) 586 587 !Integral from r1 to rmax 588 do ir=1,mesh_size 589 if (abs(rnr(ir))>1.d-20) ff(ir)=(two_pi*radmesh%rad(ir)*cos(arg*radmesh%rad(ir)) & 590 & - sin(arg*radmesh%rad(ir))/qgrid(mqgrid)) *rnr(ir) 591 end do 592 call simp_gen(r1torm,ff,radmesh) 593 594 !Integral from rmax to infinity 595 !This part is approximated using an exponential density aexp*Exp[-bexp*r] 596 !(formulae obtained with mathematica) 597 rmtoin=-one/(qgrid(mqgrid)*(arg**2+bexp**2)**3) & 598 & *aexp*exp(-bexp*rm) & 599 & *((arg**5*rm-two_pi*arg**4*qgrid(mqgrid)*rm*(bexp*rm-two) & 600 & +two*arg**3*bexp*(bexp*rm+one)+arg*bexp**3*(bexp*rm+two) & 601 & -four_pi*arg**2*bexp*qgrid(mqgrid)*(bexp**2*rm**2-three) & 602 & -two_pi*bexp**3*qgrid(mqgrid)*(bexp**2*rm**2+two*bexp*rm+two))*cos(arg*rm) & 603 & +(two*arg**2*bexp**3*rm+two_pi*arg**5*qgrid(mqgrid)*rm**2 & 604 & +arg**4*(bexp*rm-one)+bexp**4*(bexp*rm+one) & 605 & +four_pi*arg**3*qgrid(mqgrid)*(bexp**2*rm**2+two*bexp*rm-one) & 606 & +two_pi*arg*bexp**2*qgrid(mqgrid)*(bexp**2*rm**2+four*bexp*rm+6._dp))*sin(arg*rm)) 607 608 !Some of the three parts 609 ypn=two/qgrid(mqgrid)*(r0tor1+r1torm+rmtoin) 610 611 !=========================================== 612 !=== Compute 1/q dn(q)/dq at q=0 613 !=========================================== 614 615 !Integral from 0 to r1 (only if r1<>0) 616 r0tor1=zero 617 if (.not.begin_r0) r0tor1=(rnr(1)*radmesh%rad(1)**4)/5.d0 618 619 !Integral from r1 to rmax 620 do ir=1,mesh_size 621 if (abs(rnr(ir))>1.d-20) ff(ir)=rnr(ir)*radmesh%rad(ir)**3 622 end do 623 call simp_gen(r1torm,ff,radmesh) 624 625 !Integral from rmax to infinity 626 !This part is approximated using an exponential density aexp*Exp[-bexp*r] 627 !(formulae obtained with mathematica) 628 rmtoin=aexp*exp(-bexp*rm)/bexp**5 & 629 & *(24._dp+24._dp*bexp*rm+12._dp*bexp**2*rm**2+four*bexp**3*rm**3+bexp**4*rm**4) 630 631 !Some of the three parts 632 dnqdq0=-(2.d0/3.d0)*two_pi**3*(r0tor1+r1torm+rmtoin) 633 634 LIBPAW_DEALLOCATE(ff) 635 LIBPAW_DEALLOCATE(rnr) 636 637 d2nqdq0 = 1_dp 638 639 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
88 type, public :: pawpsp_header_type 89 90 !Integer scalars 91 integer :: basis_size ! Number of elements of the wf basis ((l,n) quantum numbers) 92 integer :: l_size ! Maximum value of l+1 leading to a non zero Gaunt coefficient 93 integer :: lmn_size ! Number of elements of the paw basis 94 integer :: mesh_size ! Dimension of (main) radial mesh 95 integer :: pawver ! Version number of paw psp format 96 integer :: shape_type ! Type of shape function 97 real(dp) :: rpaw ! Radius for paw spheres 98 real(dp) :: rshp ! Cut-off radius of shape function 99 100 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).
SOURCE
295 subroutine pawpsp_lo(epsatm,mqgrid,qgrid,q2vq,radmesh,vloc,yp1,ypn,zion) 296 297 !Arguments---------------------------------------------------------- 298 !scalars 299 integer,intent(in) :: mqgrid 300 real(dp),intent(in) :: zion 301 real(dp),intent(out) :: epsatm,yp1,ypn 302 type(pawrad_type),intent(in) :: radmesh 303 !arrays 304 real(dp),intent(in) :: qgrid(mqgrid) 305 real(dp),intent(in) :: vloc(:) 306 real(dp),intent(out) :: q2vq(mqgrid) 307 308 !Local variables ------------------------------ 309 !scalars 310 integer :: iq,ir,irmax,mesh_size 311 real(dp) :: arg,r0tor1,r1torm,rmtoin 312 logical :: begin_r0 313 !arrays 314 real(dp),allocatable :: ff(:),rvpz(:) 315 316 !************************************************************************ 317 318 mesh_size=size(vloc) 319 irmax=pawrad_ifromr(radmesh,min(20._dp,radmesh%rmax)) 320 irmax=min(irmax,mesh_size) 321 322 !Particular case of a zero potential 323 if (maxval(abs(vloc(1:irmax)))<=1.e-20_dp) then 324 q2vq=zero;yp1=zero;ypn=zero;epsatm=zero 325 return 326 end if 327 328 LIBPAW_ALLOCATE(ff,(mesh_size)) 329 LIBPAW_ALLOCATE(rvpz,(mesh_size)) 330 ff=zero;rvpz=zero 331 332 !Is mesh beginning with r=0 ? 333 begin_r0=(radmesh%rad(1)<1.e-20_dp) 334 335 !Store r.V+Z 336 do ir=1,irmax 337 rvpz(ir)=radmesh%rad(ir)*vloc(ir)+zion 338 end do 339 340 !=========================================== 341 !=== Compute q^2 v(q) for q=0 separately 342 !=========================================== 343 344 !Integral from 0 to r1 (only if r1<>0) 345 r0tor1=zero;if (.not.begin_r0) & 346 & r0tor1=(zion*0.5_dp+radmesh%rad(1)*vloc(1)/3._dp)*radmesh%rad(1)**2 347 348 !Integral from r1 to rmax 349 do ir=1,irmax 350 if (abs(rvpz(ir))>1.e-20_dp) then 351 ff(ir)=radmesh%rad(ir)*rvpz(ir) 352 end if 353 end do 354 355 call simp_gen(r1torm,ff,radmesh) 356 357 !Integral from rmax to infinity 358 !This part is neglected... might be improved. 359 rmtoin=zero 360 361 !Some of the three parts 362 epsatm=four_pi*(r0tor1+r1torm+rmtoin) 363 364 q2vq(1)=-zion/pi 365 366 !=========================================== 367 !=== Compute q^2 v(q) for other q''s 368 !=========================================== 369 370 !Loop over q values 371 do iq=2,mqgrid 372 arg=two_pi*qgrid(iq) 373 374 ! Integral from 0 to r1 (only if r1<>0) 375 r0tor1=zero;if (.not.begin_r0) & 376 & r0tor1=( vloc(1)/arg*sin(arg*radmesh%rad(1)) & 377 & -rvpz(1) *cos(arg*radmesh%rad(1)) +zion )/pi 378 379 ! Integral from r1 to rmax 380 do ir=1,irmax 381 if (abs(rvpz(ir))>1.e-20_dp) ff(ir)=sin(arg*radmesh%rad(ir))*rvpz(ir) 382 end do 383 call simp_gen(r1torm,ff,radmesh) 384 385 ! Integral from rmax to infinity 386 ! This part is neglected... might be improved. 387 rmtoin=zero 388 389 ! Store q^2 v(q) 390 q2vq(iq)=-zion/pi + two*qgrid(iq)*(r0tor1+r1torm+rmtoin) 391 end do 392 393 !=========================================== 394 !=== Compute derivatives of q^2 v(q) 395 !=== at ends of interval 396 !=========================================== 397 398 !yp(0)=zero 399 yp1=zero 400 401 !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]$ 402 arg=two_pi*qgrid(mqgrid) 403 404 !Integral from 0 to r1 (only if r1<>0) 405 r0tor1=zero;if (.not.begin_r0) & 406 & r0tor1=zion*radmesh%rad(1) *sin(arg*radmesh%rad(1)) & 407 & +three*radmesh%rad(1)*vloc(1)/arg *cos(arg*radmesh%rad(1)) & 408 & +(radmesh%rad(1)**2-one/arg**2)*vloc(1)*sin(arg*radmesh%rad(1)) 409 410 !Integral from r1 to rmax 411 do ir=1,irmax 412 if (abs(rvpz(ir))>1.e-20_dp) ff(ir)=( arg*radmesh%rad(ir)*cos(arg*radmesh%rad(ir)) & 413 & + sin(arg*radmesh%rad(ir))) *rvpz(ir) 414 end do 415 call simp_gen(r1torm,ff,radmesh) 416 417 !Integral from rmax to infinity 418 !This part is neglected... might be improved. 419 rmtoin=zero 420 421 !Some of the three parts 422 ypn=two*(r0tor1+r1torm+rmtoin) 423 424 LIBPAW_DEALLOCATE(ff) 425 LIBPAW_DEALLOCATE(rvpz) 426 427 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 for the calculation of XC potential
(if density<xc_denpos, density=zero)
[xc_taupos]=tolerance on kinetic energy density for the calculation of XC potential (mGGA)
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
SOURCE
4840 subroutine pawpsp_main( & 4841 & pawrad,pawtab,& 4842 & filpsp,usewvl,icoulomb,hyb_mixing,ixc,xclevel,pawxcdev,usexcnhat,& 4843 & qgrid_ff,qgrid_vl,ffspl,vlspl,epsatm,xcccrc,zionpsp,znuclpsp,& 4844 & wvl_ngauss,psxml,comm_mpi,xc_denpos,xc_taupos) ! Optional arguments 4845 4846 !Arguments ------------------------------------ 4847 !scalars 4848 integer,intent(in) :: icoulomb,ixc 4849 integer,intent(in) :: pawxcdev,usewvl,usexcnhat,xclevel 4850 integer,optional,intent(in) :: comm_mpi 4851 real(dp),intent(in):: hyb_mixing,zionpsp,znuclpsp 4852 real(dp),optional,intent(in) :: xc_denpos,xc_taupos 4853 real(dp),intent(out) :: epsatm,xcccrc 4854 character(len=fnlen),intent(in):: filpsp ! name of the psp file 4855 type(pawrad_type),intent(inout) :: pawrad 4856 type(pawtab_type),intent(inout) :: pawtab 4857 type(paw_setup_t),optional,intent(in) :: psxml 4858 !arrays 4859 integer,optional,intent(in) :: wvl_ngauss(2) 4860 real(dp),intent(in) :: qgrid_ff(:),qgrid_vl(:) 4861 real(dp),intent(inout) :: ffspl(:,:,:) 4862 real(dp),intent(out) :: vlspl(:,:) 4863 4864 !Local variables------------------------------- 4865 integer :: has_coretau,has_tproj,has_wvl,ipsp,lmax,lloc,lnmax,mmax,me,mqgrid_ff,mqgrid_vl 4866 integer :: pspcod,pspxc,usexml 4867 real(dp),parameter :: xc_denpos_default=tol14 4868 real(dp) :: my_xc_denpos,my_xc_taupos,r2well,zion,znucl 4869 character(len=500) :: msg 4870 type(pawpsp_header_type) :: pawpsp_header 4871 !arrays 4872 4873 ! ************************************************************************* 4874 4875 !Check consistency of parameters 4876 if (icoulomb/= 0.or.usewvl==1) then 4877 if (.not.present(wvl_ngauss)) then 4878 msg='usewvl==1 or icoulomb/=0: a mandatory argument is missing!' 4879 LIBPAW_BUG(msg) 4880 end if 4881 end if 4882 4883 mqgrid_ff=size(qgrid_ff) 4884 mqgrid_vl=size(qgrid_vl) 4885 lnmax=size(ffspl,3) 4886 if (size(ffspl,1)/=mqgrid_ff.or.size(ffspl,2)/=2) then 4887 msg='invalid sizes for ffspl!' 4888 LIBPAW_BUG(msg) 4889 end if 4890 if (size(vlspl,1)/=mqgrid_vl.or.size(vlspl,2)/=2) then 4891 msg='invalid sizes for vlspl!' 4892 LIBPAW_BUG(msg) 4893 end if 4894 4895 my_xc_denpos=xc_denpos_default;if (present(xc_denpos)) my_xc_denpos=xc_denpos 4896 my_xc_taupos=my_xc_denpos;if (present(xc_taupos)) my_xc_taupos=xc_taupos 4897 pawtab%usexcnhat=usexcnhat 4898 me=0;if (present(comm_mpi))me=xmpi_comm_rank(comm_mpi) 4899 4900 has_wvl=0; if (usewvl==1.or.icoulomb/=0) has_wvl=1 4901 has_tproj=0; if (usewvl==1) has_tproj=1 4902 has_coretau=0 ; if (pawxc_get_usekden(ixc)>=1) has_coretau=1 4903 call pawtab_set_flags(pawtab,has_coretau=has_coretau,has_tvale=1,has_wvl=has_wvl,has_tproj=has_tproj) 4904 4905 if(me==0) then 4906 write(msg, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(filpsp) 4907 call wrtout(ab_out, msg,'COLL') 4908 call wrtout(std_out, msg,'COLL') 4909 4910 ! This checks if file is xml or UPF 4911 ! It sets usexml as well 4912 call pawpsp_check_xml_upf(filpsp) 4913 4914 ! ---------------------------------------------------------------------------- 4915 if (usexml /= 1) then 4916 ! Open the atomic data file, and read the three first lines 4917 open (unit=tmp_unit,file=filpsp,form='formatted',status='old') 4918 rewind (unit=tmp_unit) 4919 ! Read first 3 lines of psp file: 4920 call pawpsp_read_header(tmp_unit,lloc,lmax,mmax,pspcod,& 4921 & pspxc,r2well,zion,znucl) 4922 4923 else if (usexml == 1 .and. present(psxml)) then 4924 write(msg,'(a,a)') & 4925 & '- pawpsp : Reading pseudopotential header in XML form from ', trim(filpsp) 4926 call wrtout(ab_out,msg,'COLL') 4927 call wrtout(std_out, msg,'COLL') 4928 4929 ! Return header information 4930 call pawpsp_read_header_xml(lloc,lmax,pspcod,& 4931 & pspxc,psxml,r2well,zion,znucl) 4932 ! Fill in pawpsp_header object: 4933 call pawpsp_read_pawheader(pawpsp_header%basis_size,& 4934 & lmax,pawpsp_header%lmn_size,& 4935 & pawpsp_header%l_size,pawpsp_header%mesh_size,& 4936 & pawpsp_header%pawver,psxml,& 4937 & pawpsp_header%rpaw,pawpsp_header%rshp,pawpsp_header%shape_type) 4938 end if 4939 4940 ! Check data for consistency against main routine input 4941 call pawpsp_consistency() 4942 4943 ! Read rest of the PSP file 4944 if (pspcod==7) then 4945 ! ABINIT proprietary format 4946 call pawpsp_7in(epsatm,ffspl,icoulomb,hyb_mixing,ixc,& 4947 & lmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,& 4948 & pawrad,pawtab,pawxcdev,qgrid_ff,qgrid_vl,& 4949 & usewvl,usexcnhat,vlspl,xcccrc,xclevel,my_xc_denpos,zion,znucl,xc_taupos=my_xc_taupos) 4950 4951 else if (pspcod==17)then 4952 ! XML format 4953 ipsp=1 4954 call pawpsp_17in(epsatm,ffspl,icoulomb,ipsp,hyb_mixing,ixc,lmax,& 4955 & lnmax,mmax,mqgrid_ff,mqgrid_vl,pawpsp_header,pawrad,pawtab,& 4956 & pawxcdev,qgrid_ff,qgrid_vl,usewvl,usexcnhat,vlspl,xcccrc,& 4957 & xclevel,my_xc_denpos,zion,znucl,xc_taupos=my_xc_taupos) 4958 4959 end if 4960 end if!me==0 4961 4962 close(unit=tmp_unit) 4963 4964 write(msg,'(3a)') ' pawpsp: atomic psp has been read ',& 4965 & ' and splines computed',ch10 4966 call wrtout(ab_out,msg,'COLL') 4967 call wrtout(std_out, msg,'COLL') 4968 4969 !Communicate PAW objects 4970 if(present(comm_mpi)) then 4971 if(xmpi_comm_size(comm_mpi)>1) then 4972 call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc) 4973 end if 4974 end if 4975 4976 !WVL+PAW: 4977 if(icoulomb/=0.or.usewvl==1) then 4978 if(present(comm_mpi))then 4979 call pawpsp_wvl(filpsp,pawrad,pawtab,usewvl,wvl_ngauss,comm_mpi) 4980 else 4981 call pawpsp_wvl(filpsp,pawrad,pawtab,usewvl,wvl_ngauss) 4982 end if 4983 end if 4984 4985 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]$
SOURCE
135 subroutine pawpsp_nl(ffspl,indlmn,lmnmax,lnmax,mqgrid,qgrid,radmesh,wfll) 136 137 !Arguments ------------------------------------ 138 !scalars 139 integer,intent(in) :: lmnmax,lnmax,mqgrid 140 type(pawrad_type),intent(in) :: radmesh 141 !arrays 142 integer,intent(in) :: indlmn(6,lmnmax) 143 real(dp),intent(in) :: qgrid(mqgrid) 144 real(dp),intent(in) :: wfll(:,:) 145 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) 146 147 !Local variables------------------------------- 148 !scalars 149 integer :: ilmn,iln,iln0,iq,ir,ll,meshsz,mmax 150 real(dp),parameter :: eps=tol14**4,TOLJ=0.001_dp 151 real(dp) :: arg,argn,bes 152 real(dp) :: besp,qr 153 real(dp) :: yp1,ypn 154 character(len=100) :: msg 155 type(pawrad_type) :: tmpmesh 156 !arrays 157 real(dp),allocatable :: ff(:),gg(:),rr(:),rr2(:),rr2wf(:),rrwf(:),work(:) 158 159 !************************************************************************* 160 161 !Is mesh beginning with r=0 ? 162 if (radmesh%rad(1)>tol10) then 163 msg='Radial mesh cannot begin with r<>0!' 164 LIBPAW_BUG(msg) 165 end if 166 167 meshsz=size(wfll,1) 168 if (meshsz>radmesh%mesh_size) then 169 msg='wrong size for wfll!' 170 LIBPAW_BUG(msg) 171 end if 172 173 !Init. temporary arrays and variables 174 LIBPAW_ALLOCATE(ff,(meshsz)) 175 LIBPAW_ALLOCATE(gg,(meshsz)) 176 LIBPAW_ALLOCATE(rr,(meshsz)) 177 LIBPAW_ALLOCATE(rr2,(meshsz)) 178 LIBPAW_ALLOCATE(rrwf,(meshsz)) 179 LIBPAW_ALLOCATE(rr2wf,(meshsz)) 180 LIBPAW_ALLOCATE(work,(mqgrid)) 181 rr(1:meshsz) =radmesh%rad(1:meshsz) 182 rr2(1:meshsz)=two_pi*rr(1:meshsz)*rr(1:meshsz) 183 argn=two_pi*qgrid(mqgrid) 184 mmax=meshsz 185 186 !Loop on (l,n) projectors 187 iln0=0 188 do ilmn=1,lmnmax 189 iln=indlmn(5,ilmn) 190 if(iln>iln0) then 191 iln0=iln;ll=indlmn(1,ilmn) 192 193 ir=meshsz 194 do while (abs(wfll(ir,iln))<eps) 195 ir=ir-1 196 end do 197 mmax=min(ir+1,meshsz) 198 if (mmax/=radmesh%int_meshsz) then 199 call pawrad_init(tmpmesh,mesh_size=meshsz,mesh_type=radmesh%mesh_type, & 200 & rstep=radmesh%rstep,lstep=radmesh%lstep,r_for_intg=rr(mmax)) 201 else 202 call pawrad_copy(radmesh,tmpmesh) 203 end if 204 205 rrwf(:) =rr (:)*wfll(:,iln) 206 rr2wf(:)=rr2(:)*wfll(:,iln) 207 208 ! 1-Compute f_l(0<q<qmax) 209 if (mqgrid>2) then 210 do iq=2,mqgrid-1 211 arg=two_pi*qgrid(iq) 212 do ir=1,mmax 213 qr=arg*rr(ir) 214 call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ) 215 ff(ir)=bes*rrwf(ir) 216 end do 217 call simp_gen(ffspl(iq,1,iln),ff,tmpmesh) 218 end do 219 end if 220 221 ! 2-Compute f_l(q=0) and first derivative 222 ffspl(1,1,iln)=zero;yp1=zero 223 if (ll==0) then 224 call simp_gen(ffspl(1,1,iln),rrwf,tmpmesh) 225 end if 226 if (ll==1) then 227 call simp_gen(yp1,rr2wf,tmpmesh) 228 yp1=yp1*third 229 end if 230 231 ! 3-Compute f_l(q=qmax) and first derivative 232 if (mqgrid>1) then 233 ! if (ll==0.or.ll==1) then 234 do ir=1,mmax 235 qr=argn*rr(ir) 236 call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ) 237 ff(ir)=bes*rrwf(ir) 238 gg(ir)=besp*rr2wf(ir) 239 end do 240 call simp_gen(ffspl(mqgrid,1,iln),ff,tmpmesh) 241 call simp_gen(ypn,gg,tmpmesh) 242 else 243 ypn=yp1 244 end if 245 246 ! 4-Compute second derivative of f_l(q) 247 call paw_spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 248 249 call pawrad_free(tmpmesh) 250 251 ! End loop on (l,n) projectors 252 end if 253 end do 254 255 LIBPAW_DEALLOCATE(ff) 256 LIBPAW_DEALLOCATE(gg) 257 LIBPAW_DEALLOCATE(rr) 258 LIBPAW_DEALLOCATE(rr2) 259 LIBPAW_DEALLOCATE(rrwf) 260 LIBPAW_DEALLOCATE(rr2wf) 261 LIBPAW_DEALLOCATE(work) 262 263 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
SOURCE
724 subroutine pawpsp_read(core_mesh,funit,imainmesh,lmax,& 725 & ncore,nmesh,pawrad,pawtab,pspversion,radmesh,save_core_msz,& 726 & tcoretau,tncore,tnvale,tproj,tproj_mesh,usexcnhat_in,usexcnhat_out,vale_mesh,& 727 & vlocopt,vlocr,vloc_mesh,znucl) 728 729 !Arguments ------------------------------------ 730 integer,intent(in):: funit,lmax,usexcnhat_in 731 integer,intent(out) :: imainmesh,pspversion,usexcnhat_out,vlocopt 732 logical,intent(in) :: save_core_msz 733 real(dp),intent(in):: znucl 734 !arrays 735 real(dp),pointer :: ncore(:),tcoretau(:),tncore(:),tnvale(:),tproj(:,:),vlocr(:) 736 type(pawrad_type),intent(inout) :: pawrad 737 type(pawrad_type),intent(out) :: core_mesh,tproj_mesh,vale_mesh,vloc_mesh 738 type(pawrad_type),pointer :: radmesh(:) 739 type(pawtab_type),intent(inout) :: pawtab 740 integer,intent(out)::nmesh 741 742 !Local variables------------------------------- 743 integer :: creatorid,imsh 744 integer :: icoremesh,ishpfmesh,ivalemesh,ivlocmesh 745 integer :: ib,il,ilm,ilmn,iln,iprojmesh 746 integer :: ii,ir,iread1,iread2,jj 747 integer :: msz,pngau_,ptotgau_ 748 real(dp):: rc,rread1,rread2 749 real(dp) :: yp1,ypn 750 !arrays 751 integer,allocatable :: nprj(:) 752 real(dp),allocatable :: shpf(:,:),val(:),vhnzc(:) 753 real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:) 754 character :: blank=' ',numb=' ' 755 character(len=80) :: pspline 756 character(len=500) :: msg,submsg 757 logical :: read_gauss=.false. 758 type(pawrad_type)::shpf_mesh 759 760 ! ************************************************************************* 761 762 !========================================================== 763 !Read lines 4 to 11 of the header 764 765 !This is important for BigDFT in standalone mode 766 call pawpsp_read_header_2(funit,pspversion,pawtab%basis_size,pawtab%lmn_size) 767 768 !Check pspversion for wvl-paw 769 if(pspversion<4 .and. pawtab%has_wvl>0)then 770 write(msg, '(a,i2,a,a)' )& 771 & 'In reading atomic psp file, finds pspversion=',pspversion,ch10,& 772 & 'For WVL-PAW, pspversion >= 4 is required.' 773 LIBPAW_BUG(msg) 774 end if 775 776 777 !Have to maintain compatibility with Abinit v4.2.x 778 if (pspversion==1) then 779 LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size)) 780 read(funit,*) (pawtab%orbitals(ib), ib=1,pawtab%basis_size) 781 pawtab%l_size=2*maxval(pawtab%orbitals)+1 782 nmesh=3 783 LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh)) 784 read(funit,'(a80)') pspline 785 radmesh(1)%lstep=zero 786 read(unit=pspline,fmt=*,err=10,end=10) radmesh(1)%mesh_type,& 787 & radmesh(1)%rstep,radmesh(1)%lstep 788 10 read(funit,*) pawtab%rpaw 789 read(funit,*) radmesh(1)%mesh_size,radmesh(2)%mesh_size,& 790 & radmesh(3)%mesh_size 791 read(funit,'(a80)') pspline 792 pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99 793 read(unit=pspline,fmt=*,err=11,end=11) pawtab%shape_type,& 794 & pawtab%shape_lambda,pawtab%shape_sigma 795 11 read(funit,*) creatorid 796 if (pawtab%shape_type==3) pawtab%shape_type=-1 797 radmesh(2)%mesh_type=radmesh(1)%mesh_type 798 radmesh(3)%mesh_type=radmesh(1)%mesh_type 799 radmesh(2)%rstep=radmesh(1)%rstep 800 radmesh(3)%rstep=radmesh(1)%rstep 801 radmesh(2)%lstep=radmesh(1)%lstep 802 radmesh(3)%lstep=radmesh(1)%lstep 803 else 804 805 ! Here psp file for Abinit 4.3+ 806 LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size)) 807 read(funit,*) (pawtab%orbitals(ib), ib=1,pawtab%basis_size) 808 pawtab%l_size=2*maxval(pawtab%orbitals)+1 809 read(funit,*) nmesh 810 LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh)) 811 do imsh=1,nmesh 812 rread2=zero 813 read(funit,'(a80)') pspline 814 read(unit=pspline,fmt=*,err=20,end=20) ii,iread1,iread2,rread1,rread2 815 20 continue 816 if (ii<=nmesh) then 817 radmesh(ii)%mesh_type=iread1 818 radmesh(ii)%mesh_size=iread2 819 radmesh(ii)%rstep=rread1 820 radmesh(ii)%lstep=rread2 821 else 822 write(msg, '(3a)' )& 823 & 'Index of mesh out of range !',ch10,& 824 & 'Action : check your pseudopotential file.' 825 LIBPAW_ERROR(msg) 826 end if 827 end do 828 read(funit,*) pawtab%rpaw 829 read(funit,'(a80)') pspline 830 read(unit=pspline,fmt=*) pawtab%shape_type 831 pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99 832 end if 833 834 !Initialize radial meshes 835 do imsh=1,nmesh 836 call pawrad_init(radmesh(imsh)) 837 end do 838 839 !========================================================== 840 !Initialize various dims and indexes 841 842 pawtab%l_size=2*maxval(pawtab%orbitals)+1 843 pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2 844 pawtab%ij_size=pawtab%basis_size*(pawtab%basis_size+1)/2 845 pawtab%usexcnhat=usexcnhat_in 846 847 !indlmn calculation (indices for (l,m,n) basis) 848 if (allocated(pawtab%indlmn)) then 849 LIBPAW_DEALLOCATE(pawtab%indlmn) 850 end if 851 LIBPAW_ALLOCATE(pawtab%indlmn,(6,pawtab%lmn_size)) 852 LIBPAW_BOUND1_ALLOCATE(nprj,BOUNDS(0,maxval(pawtab%orbitals))) 853 pawtab%indlmn(:,:)=0 854 ilmn=0;iln=0;nprj=0 855 do ib=1,pawtab%basis_size 856 il=pawtab%orbitals(ib) 857 nprj(il)=nprj(il)+1 858 iln=iln+1 859 do ilm=1,2*il+1 860 pawtab%indlmn(1,ilmn+ilm)=il 861 pawtab%indlmn(2,ilmn+ilm)=ilm-(il+1) 862 pawtab%indlmn(3,ilmn+ilm)=nprj(il) 863 pawtab%indlmn(4,ilmn+ilm)=il*il+ilm 864 pawtab%indlmn(5,ilmn+ilm)=iln 865 pawtab%indlmn(6,ilmn+ilm)=1 866 end do 867 ilmn=ilmn+2*il+1 868 end do 869 LIBPAW_DEALLOCATE(nprj) 870 !Are ilmn (found here) and pawtab%lmn_size compatibles ? 871 if (ilmn/=pawtab%lmn_size) then 872 write(msg, '(a,a,a,a,a)' )& 873 & 'Calculated lmn size differs from',ch10,& 874 & 'lmn_size read from pseudo !',ch10,& 875 & 'Action: check your pseudopotential file.' 876 LIBPAW_ERROR(msg) 877 end if 878 879 !========================================================== 880 !Here reading shapefunction parameters 881 882 !Shapefunction parameters for Abinit 4.3...4.5 883 if (pspversion==2) then 884 if (pawtab%shape_type==1) read(unit=pspline,fmt=*) ii,pawtab%shape_lambda,pawtab%shape_sigma 885 if (pawtab%shape_type==3) pawtab%shape_type=-1 886 pawtab%rshp=zero 887 888 !Shapefunction parameters for Abinit 4.6+ 889 else if (pspversion>=3) then 890 pawtab%rshp=zero 891 if (pawtab%shape_type==-1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp 892 if (pawtab%shape_type== 1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp, & 893 & pawtab%shape_lambda,pawtab%shape_sigma 894 if (pawtab%shape_type== 2) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp 895 if (pawtab%shape_type== 3) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp 896 end if 897 21 continue 898 !If shapefunction type is gaussian, check exponent 899 if (pawtab%shape_type==1) then 900 if (pawtab%shape_lambda<2) then 901 write(msg, '(3a)' )& 902 & 'For a gaussian shape function, exponent lambda must be >1 !',ch10,& 903 & 'Action: check your psp file.' 904 LIBPAW_ERROR(msg) 905 end if 906 end if 907 !If shapefunction type is Bessel, deduce here its parameters from rc 908 if (pawtab%shape_type==3) then 909 LIBPAW_ALLOCATE(pawtab%shape_alpha,(2,pawtab%l_size)) 910 LIBPAW_ALLOCATE(pawtab%shape_q,(2,pawtab%l_size)) 911 rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw 912 do il=1,pawtab%l_size 913 call atompaw_shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc) 914 end do 915 end if 916 917 !========================================================== 918 !Mirror pseudopotential parameters to the output and log files 919 920 write(msg,'(a,i1)')' Pseudopotential format is: paw',pspversion 921 call wrtout(ab_out,msg,'COLL') 922 call wrtout(std_out, msg,'COLL') 923 write(msg,'(2(a,i3),a,64i4)') & 924 & ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',& 925 & pawtab%lmn_size,'), orbitals=',pawtab%orbitals(1:pawtab%basis_size) 926 call wrtout(ab_out,msg,'COLL') 927 call wrtout(std_out, msg,'COLL') 928 write(msg,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw 929 call wrtout(ab_out,msg,'COLL') 930 call wrtout(std_out, msg,'COLL') 931 write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:' 932 call wrtout(ab_out,msg,'COLL') 933 call wrtout(std_out, msg,'COLL') 934 do imsh=1,nmesh 935 if (radmesh(imsh)%mesh_type==1) & 936 & write(msg,'(a,i1,a,i4,a,g12.5)') & 937 & ' - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,& 938 & ' , step=',radmesh(imsh)%rstep 939 if (radmesh(imsh)%mesh_type==2) & 940 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 941 & ' - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,& 942 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 943 if (radmesh(imsh)%mesh_type==3) & 944 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 945 & ' - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,& 946 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 947 if (radmesh(imsh)%mesh_type==4) & 948 & write(msg,'(a,i1,a,i4,a,g12.5)') & 949 & ' - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,& 950 & ' , AA=',radmesh(imsh)%rstep 951 call wrtout(ab_out,msg,'COLL') 952 call wrtout(std_out, msg,'COLL') 953 end do 954 if (pawtab%shape_type==-1) then 955 write(msg,'(a)')& 956 ' Shapefunction is NUMERIC type: directly read from atomic data file' 957 call wrtout(ab_out,msg,'COLL') 958 call wrtout(std_out, msg,'COLL') 959 end if 960 if (pawtab%shape_type==1) then 961 write(msg,'(2a,a,f6.3,a,i3)')& 962 & ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,& 963 & ' with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda 964 call wrtout(ab_out,msg,'COLL') 965 call wrtout(std_out, msg,'COLL') 966 end if 967 if (pawtab%shape_type==2) then 968 write(msg,'(a)')& 969 ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2' 970 call wrtout(ab_out,msg,'COLL') 971 call wrtout(std_out, msg,'COLL') 972 end if 973 if (pawtab%shape_type==3) then 974 write(msg,'(a)')& 975 & ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)' 976 call wrtout(ab_out,msg,'COLL') 977 call wrtout(std_out, msg,'COLL') 978 end if 979 if (pawtab%rshp<1.d-8) then 980 write(msg,'(a)') ' Radius for shape functions = sphere core radius' 981 else 982 write(msg,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp 983 end if 984 call wrtout(ab_out,msg,'COLL') 985 call wrtout(std_out, msg,'COLL') 986 987 !========================================================== 988 !Perfom tests 989 990 !Are lmax and orbitals compatibles ? 991 if (lmax/=maxval(pawtab%orbitals)) then 992 write(msg, '(a,a,a)' )& 993 & 'lmax /= MAX(orbitals) !',ch10,& 994 & 'Action: check your pseudopotential file.' 995 LIBPAW_ERROR(msg) 996 end if 997 998 !Only mesh_type=1,2, 3 or 4 allowed 999 do imsh=1,nmesh 1000 if (radmesh(imsh)%mesh_type>4) then 1001 write(msg, '(a,a,a)' )& 1002 & 'Only mesh types 1,2, 3 or 4 allowed !',ch10,& 1003 & 'Action : check your pseudopotential or input file.' 1004 LIBPAW_ERROR(msg) 1005 end if 1006 end do 1007 1008 !========================================================== 1009 !Read tabulated atomic data 1010 1011 !--------------------------------- 1012 !Read wave-functions (phi) 1013 do ib=1,pawtab%basis_size 1014 read (funit,*) 1015 if (pspversion==1) iread1=1 1016 if (pspversion>1) read (funit,*) iread1 1017 if (ib==1) then 1018 call pawrad_free(pawrad) 1019 call pawrad_init(pawrad,mesh_size=radmesh(iread1)%mesh_size,mesh_type=radmesh(iread1)%mesh_type,& 1020 & rstep=radmesh(iread1)%rstep,lstep=radmesh(iread1)%lstep,r_for_intg=pawtab%rpaw) 1021 pawtab%partialwave_mesh_size=pawrad%mesh_size 1022 pawtab%mesh_size=pawrad_ifromr(pawrad,pawtab%rpaw)+5 1023 pawtab%mesh_size=min(pawtab%mesh_size,pawrad%mesh_size) 1024 if (pawtab%mesh_size>pawrad%mesh_size-2) pawtab%mesh_size=pawrad%mesh_size 1025 imainmesh=iread1 1026 LIBPAW_ALLOCATE(pawtab%phi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 1027 else if (iread1/=imainmesh) then 1028 write(msg, '(a,a,a)' )& 1029 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 1030 & 'Action: check your pseudopotential file.' 1031 LIBPAW_ERROR(msg) 1032 end if 1033 read (funit,*) (pawtab%phi(ir,ib),ir=1,pawtab%partialwave_mesh_size) 1034 end do 1035 1036 !--------------------------------- 1037 !Read pseudo wave-functions (tphi) 1038 LIBPAW_ALLOCATE(pawtab%tphi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 1039 do ib=1,pawtab%basis_size 1040 read (funit,*) 1041 if (pspversion==1) iread1=1 1042 if (pspversion>1) read (funit,*) iread1 1043 if (iread1/=imainmesh) then 1044 write(msg, '(a,a,a)' )& 1045 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 1046 & 'Action: check your pseudopotential file.' 1047 LIBPAW_ERROR(msg) 1048 end if 1049 read (funit,*) (pawtab%tphi(ir,ib),ir=1,pawtab%partialwave_mesh_size) 1050 end do 1051 write(msg,'(a,i1)') & 1052 & ' Radial grid used for partial waves is grid ',imainmesh 1053 call wrtout(ab_out,msg,'COLL') 1054 call wrtout(std_out, msg,'COLL') 1055 1056 !--------------------------------- 1057 !Read projectors (tproj) 1058 do ib=1,pawtab%basis_size 1059 read (funit,*) 1060 if (pspversion==1) iread1=2 1061 if (pspversion>1) read (funit,*) iread1 1062 if (ib==1) then 1063 iprojmesh=iread1 1064 call pawrad_copy(radmesh(iprojmesh),tproj_mesh) 1065 LIBPAW_POINTER_ALLOCATE(tproj,(tproj_mesh%mesh_size,pawtab%basis_size)) 1066 else if (iread1/=iprojmesh) then 1067 write(msg, '(a,a,a)' )& 1068 & 'All tprojectors must be given on the same radial mesh !',ch10,& 1069 & 'Action: check your pseudopotential file.' 1070 LIBPAW_ERROR(msg) 1071 end if 1072 ! read projectors from a mesh 1073 read (funit,*) (tproj(ir,ib),ir=1,tproj_mesh%mesh_size) 1074 end do 1075 write(msg,'(a,i2)') & 1076 & ' Radial grid used for projectors is grid ',iprojmesh 1077 call wrtout(ab_out,msg,'COLL') 1078 call wrtout(std_out, msg,'COLL') 1079 1080 !--------------------------------- 1081 !Read gaussian projectors for wavelets 1082 ! -- only if pawtab%has_wvl flag is on 1083 ! -- if not, we skip the lines 1084 read(funit,'(a80)') pspline 1085 if(index(trim(pspline),'GAUSSIAN')/=0) read_gauss=.true. 1086 if (read_gauss) then 1087 if (pawtab%has_wvl>0) then 1088 call wvlpaw_allocate(pawtab%wvl) 1089 jj=0 1090 do ib=1,pawtab%basis_size 1091 if(ib/=1) read(funit,*) pspline 1092 ! read Gaussian coefficients 1093 read(funit,*) pngau_, ptotgau_ !total number of gaussians 1094 if(ib==1) then 1095 pawtab%wvl%ptotgau=ptotgau_ 1096 LIBPAW_ALLOCATE(pawtab%wvl%pngau,(pawtab%basis_size)) 1097 LIBPAW_ALLOCATE(pawtab%wvl%parg,(2,pawtab%wvl%ptotgau)) 1098 LIBPAW_ALLOCATE(pawtab%wvl%pfac,(2,pawtab%wvl%ptotgau)) 1099 else 1100 if(pawtab%wvl%ptotgau/=ptotgau_) then 1101 write(msg,'(3a)')& 1102 & 'Total number of gaussians, should be the same for all projectors !',ch10,& 1103 & 'Action: check your pseudopotential file.' 1104 LIBPAW_ERROR(msg) 1105 end if 1106 end if !ib==1 1107 read(funit,*)(pawtab%wvl%parg(:,ii),ii=jj+1,jj+pngau_) 1108 read(funit,*)(pawtab%wvl%pfac(:,ii),ii=jj+1,jj+pngau_) 1109 pawtab%wvl%pngau(ib)=pngau_ 1110 jj=jj+pngau_ 1111 end do 1112 pawtab%has_wvl=2 1113 else 1114 ! If pawtab%has_wvl=0, we skip the lines 1115 do ib=1,pawtab%basis_size 1116 if(ib/=1) read(funit,*) 1117 read(funit,*) pngau_, ptotgau_ 1118 LIBPAW_ALLOCATE(val, (pngau_ *2)) 1119 read(funit,*) val 1120 read(funit,*) val 1121 LIBPAW_DEALLOCATE(val) 1122 end do 1123 end if 1124 end if 1125 1126 !--------------------------------- 1127 !Read core density (coredens) 1128 if(read_gauss) read (funit,*) !if not read_gauss, this line was already read 1129 if (pspversion==1) iread1=1 1130 if (pspversion>1) read (funit,*) iread1 1131 icoremesh=iread1 1132 call pawrad_copy(radmesh(icoremesh),core_mesh) 1133 if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.& 1134 & (radmesh(icoremesh)%rstep /=pawrad%rstep) .or.& 1135 & (radmesh(icoremesh)%lstep /=pawrad%lstep)) then 1136 write(msg, '(a,a,a,a,a)' )& 1137 & 'Ncore must be given on a radial mesh with the same',ch10,& 1138 & 'type and step(s) than the main radial mesh (mesh for Phi) !',ch10,& 1139 & 'Action: check your pseudopotential file.' 1140 LIBPAW_ERROR(msg) 1141 end if 1142 LIBPAW_POINTER_ALLOCATE(ncore,(core_mesh%mesh_size)) 1143 read (funit,*) (ncore(ir),ir=1,core_mesh%mesh_size) 1144 1145 !Construct and save VH[z_NC] if requested 1146 if (pawtab%has_vhnzc==1) then 1147 LIBPAW_ALLOCATE(pawtab%VHnZC,(pawtab%mesh_size)) 1148 LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size)) 1149 call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl) 1150 pawtab%VHnZC(1:pawtab%mesh_size)=vhnzc(1:pawtab%mesh_size) 1151 pawtab%has_vhnzc=2 1152 LIBPAW_DEALLOCATE(vhnzc) 1153 end if 1154 1155 pawtab%core_mesh_size=pawrad%mesh_size 1156 if(save_core_msz) pawtab%core_mesh_size=core_mesh%mesh_size 1157 LIBPAW_ALLOCATE(pawtab%coredens,(pawtab%core_mesh_size)) 1158 pawtab%rcore=core_mesh%rad(pawtab%core_mesh_size) 1159 pawtab%coredens(1:pawtab%core_mesh_size)=ncore(1:pawtab%core_mesh_size) 1160 1161 !--------------------------------- 1162 !Read pseudo core density (tcoredens) 1163 if(save_core_msz) then 1164 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,6)) 1165 else 1166 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,1)) 1167 end if 1168 pawtab%tcoredens=zero 1169 read (funit,*) 1170 if (pspversion==1) iread1=1 1171 if (pspversion>1) read (funit,*) iread1 1172 if (iread1/=icoremesh) then 1173 write(msg, '(a,a,a,a,a,a,a,a)' )& 1174 & 'Pseudized core density (tNcore) must be given',ch10,& 1175 & 'on the same radial mesh as core density (Ncore) !',ch10,& 1176 & 'Action: check your pseudopotential file.' 1177 LIBPAW_ERROR(msg) 1178 end if 1179 LIBPAW_POINTER_ALLOCATE(tncore,(core_mesh%mesh_size)) 1180 read (funit,*) (tncore(ir),ir=1,core_mesh%mesh_size) 1181 if (maxval(abs(tncore(:)))<tol6) then 1182 pawtab%usetcore=0 1183 else 1184 pawtab%usetcore=1 1185 pawtab%tcoredens(1:pawtab%core_mesh_size,1)=tncore(1:pawtab%core_mesh_size) 1186 end if 1187 write(msg,'(a,i1)') & 1188 & ' Radial grid used for (t)core density is grid ',icoremesh 1189 call wrtout(ab_out,msg,'COLL') 1190 call wrtout(std_out, msg,'COLL') 1191 1192 !--------------------------------- 1193 !Read frozen part of Dij terms (dij0) 1194 LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size)) 1195 read (funit,*) 1196 read (funit,*) (pawtab%dij0(ib),ib=1,pawtab%lmn2_size) 1197 1198 !--------------------------------- 1199 !Read initial guess of rhoij (rhoij0) 1200 LIBPAW_ALLOCATE(pawtab%rhoij0,(pawtab%lmn2_size)) 1201 read (funit,*) 1202 read (funit,*) (pawtab%rhoij0(ib),ib=1,pawtab%lmn2_size) 1203 1204 !--------------------------------- 1205 !Read local pseudopotential=Vh(tn_zc) or Vbare 1206 read (funit,*) 1207 if (pspversion==1) ivlocmesh=3 1208 vlocopt=1 1209 if (pspversion==2) then 1210 read (funit,*) ivlocmesh 1211 else if (pspversion>2) then 1212 ! read (funit,fmt=*,err=30,end=30) ivlocmesh,vlocopt 1213 msg=blank 1214 read (funit,fmt='(a)') msg 1215 read (msg,fmt=*) ivlocmesh 1216 write(numb,'(i1)')ivlocmesh 1217 ii=index(msg,numb) 1218 if(len_trim(trim(msg(ii+1:)))/=0)then 1219 submsg=trim(msg(ii+1:)) 1220 if(len_trim(submsg)/=0)then 1221 do ii=1,len_trim(submsg) 1222 numb=submsg(ii:ii) 1223 if(numb==blank)cycle 1224 jj=index('0123456789',numb) 1225 if(jj<1 .or. jj>10)exit 1226 vlocopt=jj-1 1227 end do 1228 end if 1229 end if 1230 end if 1231 usexcnhat_out=0;if (vlocopt==1) usexcnhat_out=1 1232 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 1233 LIBPAW_POINTER_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 1234 read (funit,*) (vlocr(ir),ir=1,vloc_mesh%mesh_size) 1235 write(msg,'(a,i1)') & 1236 & ' Radial grid used for Vloc is grid ',ivlocmesh 1237 call wrtout(ab_out,msg,'COLL') 1238 call wrtout(std_out, msg,'COLL') 1239 1240 !--------------------------------- 1241 !Eventually read "numeric" shapefunctions (if shape_type=-1) 1242 if (pawtab%shape_type==-1) then 1243 LIBPAW_ALLOCATE(pawtab%shapefunc,(pawtab%mesh_size,pawtab%l_size)) 1244 do il=1,pawtab%l_size 1245 read (funit,*) 1246 if (pspversion==1) iread1=1 1247 if (pspversion>1) read (funit,*) iread1 1248 if (il==1) then 1249 call pawrad_copy(radmesh(iread1),shpf_mesh) 1250 ishpfmesh=iread1 1251 LIBPAW_ALLOCATE(shpf,(shpf_mesh%mesh_size,pawtab%l_size)) 1252 else if (iread1/=ishpfmesh) then 1253 write(msg, '(a,a,a)' )& 1254 & 'All shape functions must be given on the same radial mesh !',ch10,& 1255 & 'Action: check your pseudopotential file.' 1256 LIBPAW_ERROR(msg) 1257 end if 1258 read (funit,*) (shpf(ir,il),ir=1,shpf_mesh%mesh_size) 1259 end do 1260 write(msg,'(a,i1)') & 1261 & ' Radial grid used for shape functions is grid ',iread1 1262 call wrtout(ab_out,msg,'COLL') 1263 call wrtout(std_out, msg,'COLL') 1264 1265 ! Has to spline shape functions if mesh is not the "main" mesh 1266 if (ishpfmesh/=imainmesh) then 1267 msz=shpf_mesh%mesh_size 1268 LIBPAW_ALLOCATE(work1,(msz)) 1269 LIBPAW_ALLOCATE(work2,(msz)) 1270 LIBPAW_ALLOCATE(work3,(msz)) 1271 LIBPAW_ALLOCATE(work4,(pawtab%mesh_size)) 1272 work3(1:pawtab%mesh_size)=shpf_mesh%rad(1:pawtab%mesh_size) 1273 work4(1:pawtab%mesh_size)=pawrad%rad(1:pawtab%mesh_size) 1274 do il=1,pawtab%l_size 1275 call bound_deriv(shpf(1:msz,il),shpf_mesh,msz,yp1,ypn) 1276 call paw_spline(work3,shpf(:,il),msz,yp1,ypn,work1) 1277 call paw_splint(msz,work3,shpf(:,il),work1,pawtab%mesh_size,work4,pawtab%shapefunc(:,il)) 1278 end do 1279 LIBPAW_DEALLOCATE(work1) 1280 LIBPAW_DEALLOCATE(work2) 1281 LIBPAW_DEALLOCATE(work3) 1282 LIBPAW_DEALLOCATE(work4) 1283 else 1284 pawtab%shapefunc(:,:)=shpf(:,:) 1285 end if 1286 LIBPAW_DEALLOCATE(shpf) 1287 call pawrad_free(shpf_mesh) 1288 end if 1289 1290 !--------------------------------- 1291 !Read pseudo valence density (if psp version >=4) 1292 if (pspversion>=4) then 1293 read (funit,*) 1294 read (funit,*) iread1 1295 ivalemesh=iread1 1296 call pawrad_copy(radmesh(iread1),vale_mesh) 1297 LIBPAW_POINTER_ALLOCATE(tnvale,(vale_mesh%mesh_size)) 1298 read (funit,*) (tnvale(ir),ir=1,vale_mesh%mesh_size) 1299 pawtab%has_tvale=1 1300 write(msg,'(a,i1)') & 1301 & ' Radial grid used for pseudo valence density is grid ',ivalemesh 1302 call wrtout(ab_out,msg,'COLL') 1303 call wrtout(std_out, msg,'COLL') 1304 else 1305 pawtab%has_tvale=0 1306 LIBPAW_POINTER_ALLOCATE(tnvale,(0)) 1307 end if 1308 1309 !--------------------------------- 1310 !Initialize (to zero) kinetic energy densities (for testing purpose) 1311 if (pawtab%has_coretau>0) then 1312 write(msg,'(5a)' )& 1313 & 'Kinetic energy density is requested but the core kinetic energy density',ch10,& 1314 & 'is not present in the pseudopotential file!',ch10,& 1315 & 'We assume that it is zero (for testing purpose).' 1316 LIBPAW_WARNING(msg) 1317 pawtab%coretau_mesh_size=pawtab%mesh_size 1318 if(save_core_msz) pawtab%coretau_mesh_size=core_mesh%mesh_size 1319 LIBPAW_ALLOCATE(pawtab%coretau,(pawtab%coretau_mesh_size)) 1320 LIBPAW_ALLOCATE(pawtab%tcoretau,(pawtab%coretau_mesh_size)) 1321 LIBPAW_POINTER_ALLOCATE(tcoretau,(core_mesh%mesh_size)) 1322 pawtab%rcoretau=core_mesh%rad(pawtab%coretau_mesh_size) 1323 pawtab%coretau=zero ; pawtab%tcoretau=zero ; tcoretau=zero 1324 endif 1325 1326 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
SOURCE
1347 subroutine pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,radmesh,phi_cor,& 1348 & filename,kappacor) ! optional arguments 1349 1350 !Arguments ------------------------------------ 1351 integer,intent(out) :: lmncmax,nphicor 1352 character(len=*),optional :: filename 1353 !arrays 1354 integer,allocatable,intent(inout) :: indlmn_core(:,:),lcor(:),ncor(:) 1355 integer,allocatable,intent(inout),optional :: kappacor(:) 1356 real(dp),allocatable,intent(inout) :: phi_cor(:,:),energy_cor(:) 1357 type(pawrad_type),intent(in) :: radmesh 1358 1359 !Local variables------------------------------- 1360 integer :: ib,i1,i2,il,im,ilm,ilmn,iln,ios,jln 1361 integer :: nmesh,npts,unt,flagrel,tmp1,tmp2,tmp3,kappa,spinor,i2j,i2mj 1362 real(dp) :: noccor,r1,r2 1363 logical :: ex,oldformat,usexml,diracrel 1364 character(len=8) :: dum,dum1,dum2,dum3,dum4 1365 character(len=80) :: fline 1366 character(len=500) :: msg 1367 character(len=fnlen) :: filename_ 1368 1369 !arrays 1370 integer,allocatable :: meshtp(:),meshsz(:) 1371 real(dp),allocatable :: rad(:),radstp(:),work(:) 1372 real(dp),allocatable :: logstp(:),phitmp(:) 1373 type(pawrad_type) :: tmpmesh 1374 1375 ! ************************************************************************ 1376 1377 !Check for core WF file existence and XML format 1378 usexml=.false.;oldformat=.false. 1379 if (present(filename)) then 1380 ! Core WF file given as optional argument 1381 filename_=filename;ex=.false. 1382 inquire(file=trim(filename_),iostat=ios,exist=ex) 1383 if (ios/=0) then 1384 write(msg,'(2a)') 'INQUIRE returns an error for file ',trim(filename_) 1385 LIBPAW_ERROR(msg) 1386 end if 1387 if (.not.ex) then 1388 write(msg,'(3a)') 'This file does not exist: ',trim(filename_),'!' 1389 LIBPAW_ERROR(msg) 1390 end if 1391 unt = libpaw_get_free_unit() 1392 open(unit=unt,file=trim(filename_),form='formatted',status='old', action="read") 1393 read(unt,*) fline 1394 close(unt) 1395 usexml=(fline(1:5)=='<?xml') 1396 else 1397 ! Core WF file: new format 1398 filename_='corewf.abinit';ex=.false. 1399 inquire(file=trim(filename_),iostat=ios,exist=ex) 1400 if (ios/=0) then 1401 write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!' 1402 LIBPAW_ERROR(msg) 1403 end if 1404 if (.not.ex) then 1405 ! Core WF file: new format XML 1406 filename_='corewf.xml';ex=.false. 1407 inquire(file=trim(filename_),iostat=ios,exist=ex) 1408 if (ios/=0) then 1409 write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!' 1410 LIBPAW_ERROR(msg) 1411 end if 1412 usexml=ex 1413 if (.not.ex) then 1414 ! Core WF file: old format 1415 filename_='corewf.dat';ex=.false. 1416 inquire(file=trim(filename_),iostat=ios,exist=ex) 1417 if (ios/=0) then 1418 write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!' 1419 LIBPAW_ERROR(msg) 1420 end if 1421 oldformat=ex 1422 if (.not.ex) then 1423 ! No core WF file found 1424 write(msg, '(3a)' )& 1425 & 'Checks for existence of files corewf.abinit[.xml] or corewf.dat',ch10,& 1426 & 'but INQUIRE finds file does not exist!' 1427 LIBPAW_ERROR(msg) 1428 end if 1429 end if 1430 end if 1431 end if 1432 1433 !Core WF file is in new XML format 1434 if ((.not.oldformat).and.(usexml)) then 1435 if(present(kappacor)) then 1436 call rdpawpsxml_core(energy_cor,trim(filename_),lcor,ncor,nphicor,radmesh,phi_cor, & 1437 & kappacor=kappacor) 1438 else 1439 call rdpawpsxml_core(energy_cor,trim(filename_),lcor,ncor,nphicor,radmesh,phi_cor) 1440 endif 1441 endif 1442 1443 !Core WF file is in (proprietary) format 1444 if ((.not.oldformat).and.(.not.usexml)) then 1445 unt = libpaw_get_free_unit() 1446 open(unt,file=trim(filename_),form='formatted',action="read") 1447 read(unt,*) ! skip title 1448 1449 diracrel=.false. 1450 read(unit=unt,fmt=*,err=23,end=23) flagrel,tmp1,tmp2,tmp3 1451 if (flagrel==2) diracrel=.true. 1452 23 continue 1453 1454 if(present(kappacor).and.(.not.diracrel)) then 1455 write(msg,'(3a)') 'Error in pawpsp_read_core:',ch10, & 1456 & ' Diracrel corewf file has to be provided!' 1457 LIBPAW_ERROR(msg) 1458 endif 1459 if((.not.present(kappacor)).and.diracrel) then 1460 write(msg,'(3a)') 'Error in pawpsp_read_core:',ch10, & 1461 & ' Cannot use diracrelativistic corewf file!' 1462 ABI_ERROR(msg) 1463 endif 1464 1465 read(unt,*) ! skip zatom,zcore,pspdat 1466 read(unt,*) ! skip pspcod,pspxc,lmax 1467 read(unt,*) ! skip pspfmt,creatorID 1468 read(unt,*) nphicor 1469 read(unt,*) ! skip orbitals 1470 read(unt,*) nmesh 1471 LIBPAW_ALLOCATE(meshsz,(nmesh)) 1472 LIBPAW_ALLOCATE(meshtp,(nmesh)) 1473 LIBPAW_ALLOCATE(radstp,(nmesh)) 1474 LIBPAW_ALLOCATE(logstp,(nmesh)) 1475 do iln=1,nmesh 1476 r2=zero;read(unt,'(a80)') fline 1477 read(unit=fline,fmt=*,err=20,end=20) ib,i1,i2,r1,r2 1478 20 continue 1479 if (ib<=nmesh) then 1480 meshtp(ib)=i1;meshsz(ib)=i2 1481 radstp(ib)=r1;logstp(ib)=r2 1482 end if 1483 end do 1484 read(unt,*) ! skip rmax(core) 1485 LIBPAW_ALLOCATE(ncor,(nphicor)) 1486 LIBPAW_ALLOCATE(lcor,(nphicor)) 1487 LIBPAW_ALLOCATE(energy_cor,(nphicor)) 1488 LIBPAW_ALLOCATE(phi_cor,(radmesh%mesh_size,nphicor)) 1489 if(present(kappacor).and.diracrel) then 1490 LIBPAW_ALLOCATE(kappacor,(nphicor)) 1491 endif 1492 do iln=1,nphicor 1493 read(unt,*) ! skip comment 1494 read(unt,*) i1 1495 if(present(kappacor).and.diracrel) then 1496 read(unt,*) ncor(iln),lcor(iln),kappacor(iln) 1497 else 1498 read(unt,*) ncor(iln),lcor(iln) 1499 endif 1500 read(unt,*) energy_cor(iln) 1501 energy_cor(iln)=energy_cor(iln)*half ! For consistency reasons (in the legacy coreWF format, energies are in Ry) 1502 LIBPAW_ALLOCATE(phitmp,(meshsz(i1))) 1503 read(unt,*) phitmp 1504 if ((radmesh%mesh_type/=meshtp(i1)) & 1505 & .or.(radmesh%rstep/=radstp(i1)) & 1506 & .or.(radmesh%lstep/=logstp(i1))) then 1507 call pawrad_init(tmpmesh,mesh_size=meshsz(i1),mesh_type=meshtp(i1),rstep=radstp(i1),lstep=logstp(i1)) 1508 npts=radmesh%mesh_size 1509 if (tmpmesh%rmax<radmesh%rmax+tol8) npts=pawrad_ifromr(radmesh,tmpmesh%rmax)-1 1510 LIBPAW_ALLOCATE(work,(meshsz(i1))) 1511 call bound_deriv(phitmp,tmpmesh,meshsz(i1),r1,r2) 1512 call paw_spline(tmpmesh%rad,phitmp,meshsz(i1),r1,r2,work) 1513 call paw_splint(meshsz(i1),tmpmesh%rad,phitmp,work,npts,radmesh%rad(1:npts),phi_cor(1:npts,iln)) 1514 if (npts<radmesh%mesh_size) phi_cor(npts+1:radmesh%mesh_size,iln)=zero 1515 LIBPAW_DEALLOCATE(work) 1516 call pawrad_free(tmpmesh) 1517 else 1518 npts=min(meshsz(i1),radmesh%mesh_size) 1519 phi_cor(1:npts,iln)=phitmp(1:npts) 1520 if (npts<radmesh%mesh_size) phi_cor(npts+1:radmesh%mesh_size,iln)=zero 1521 end if 1522 LIBPAW_DEALLOCATE(phitmp) 1523 end do 1524 LIBPAW_DEALLOCATE(meshsz) 1525 LIBPAW_DEALLOCATE(meshtp) 1526 LIBPAW_DEALLOCATE(radstp) 1527 LIBPAW_DEALLOCATE(logstp) 1528 end if 1529 1530 close(unt) 1531 1532 !Core WF file is in old (proprietary) format 1533 if ((oldformat).and.(.not.usexml)) then 1534 unt = libpaw_get_free_unit() 1535 open(unt,file=trim(filename_),form='formatted',action="read") 1536 do while (dum/='atompaw ') 1537 read(unt,'(a8)') dum 1538 end do 1539 read(unt,'(2i4)') npts,nphicor 1540 LIBPAW_ALLOCATE(ncor,(nphicor)) 1541 LIBPAW_ALLOCATE(lcor,(nphicor)) 1542 LIBPAW_ALLOCATE(energy_cor,(nphicor)) 1543 LIBPAW_ALLOCATE(phi_cor,(npts,nphicor)) 1544 LIBPAW_ALLOCATE(rad,(npts)) 1545 do iln=1,nphicor 1546 read(unt,'(a4,i4,a3,i4,a6,f15.7,a8,f15.7)') & 1547 & dum1,ncor(iln),dum2,lcor(iln),dum3,noccor,dum4,energy_cor(iln) 1548 energy_cor(iln)=energy_cor(iln)*half ! For consistency reasons (in the legacy coreWF format, energies are in Ry) 1549 1550 do jln=1,npts 1551 read(unt,*) rad(jln),phi_cor(jln,iln) 1552 end do 1553 read(unt,*) 1554 end do 1555 LIBPAW_DEALLOCATE(rad) 1556 close(unt) 1557 end if 1558 1559 !Set an array 'a la' indlmn 1560 1561 1562 !===== DIRAC-RELATIVISTIC CASE ===== 1563 !Warning due to the nature of the dirac-relativistic solution used: 1564 ! These corewf have complex spherical harmonics, which need to be converted later! 1565 if(present(kappacor)) then 1566 1567 lmncmax=0 1568 do ib=1,nphicor 1569 il=lcor(ib) 1570 kappa=sign(1,kappacor(ib)) 1571 i2j=2*il-kappa!j=l-sgn(kappa)/2 1572 lmncmax=lmncmax+i2j+1 1573 end do 1574 lmncmax=lmncmax*2 1575 LIBPAW_ALLOCATE(indlmn_core,(8,lmncmax)) 1576 indlmn_core=0;ilmn=0;iln=0 1577 do ib=1,2*nphicor 1578 iln=iln+modulo(ib,2) 1579 il=lcor(iln) 1580 kappa=sign(1,kappacor(iln)) ! sgn(kappa)=+1 or -1 1581 spinor=2-modulo(ib,2) ! spinor= 1 or 2 1582 i2j=2*il-kappa ! j=l-sgn(kappa)/2 = l-1/2 or l+1/2 1583 do ilm=1,i2j+1 1584 !mj= -j,...,j 1585 i2mj=-i2j+2*(ilm-1) ! 2m_j= -jc ... +jc 1586 im=(i2mj-3+2*spinor)/2 ! m=m_j-1/2 (spinor=1) or m_j+1/2 (spinor=2) 1587 if(abs(im)<=il) then 1588 !Valid value for sph. harm., i.e. abs(m)<=l 1589 indlmn_core(1,ilmn+ilm)=il !l 1590 indlmn_core(2,ilmn+ilm)=im !m 1591 indlmn_core(3,ilmn+ilm)=kappa !sign of kappa 1592 indlmn_core(4,ilmn+ilm)=il*il+im+il+1 !lm 1593 indlmn_core(5,ilmn+ilm)=iln !ln also includes the two kappa values here 1594 indlmn_core(6,ilmn+ilm)=spinor !spinor index (1 up, 2 down) 1595 indlmn_core(7,ilmn+ilm)=i2j !2*j (times 2 to make it an integer) 1596 indlmn_core(8,ilmn+ilm)=i2mj !2*m_j (times 2 to make it an integer) 1597 else 1598 !Invalid value for sph. harm. ; will be multiplied by zero later 1599 indlmn_core(1,ilmn+ilm)=-1 !Invalid value that should be checked later 1600 indlmn_core(2:8,ilmn+ilm)=-1 ; indlmn_core(3,ilmn+ilm)=0 1601 endif 1602 end do 1603 ilmn=ilmn+i2j+1 1604 end do 1605 1606 !===== NON OR SCALAR-RELATIVISTIC CASE ===== 1607 else 1608 lmncmax=0 1609 do ib=1,nphicor 1610 il=lcor(ib) 1611 lmncmax=lmncmax+2*il+1 1612 end do 1613 LIBPAW_ALLOCATE(indlmn_core,(6,lmncmax)) 1614 indlmn_core=0;ilmn=0;iln=0 1615 do ib=1,nphicor 1616 il=lcor(ib) 1617 iln=iln+1 1618 do ilm=1,2*il+1 1619 indlmn_core(1,ilmn+ilm)=il 1620 indlmn_core(2,ilmn+ilm)=ilm-(il+1) 1621 indlmn_core(3,ilmn+ilm)=1 1622 indlmn_core(4,ilmn+ilm)=il*il+ilm 1623 indlmn_core(5,ilmn+ilm)=iln 1624 indlmn_core(6,ilmn+ilm)=1 1625 end do 1626 ilmn=ilmn+2*il+1 1627 end do 1628 1629 endif ! Relativistic? 1630 1631 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
SOURCE
4210 subroutine pawpsp_read_header(funit,lloc,lmax,mmax,pspcod,pspxc,r2well,zion,znucl) 4211 4212 !Arguments ------------------------------------ 4213 !scalars 4214 integer,intent(in):: funit 4215 integer,intent(out):: lloc,lmax,mmax,pspcod,pspxc 4216 real(dp),intent(out):: r2well,zion,znucl 4217 !Local variables------------------------------- 4218 integer:: pspdat 4219 character(len=fnlen):: title 4220 character(len=500) :: msg 4221 4222 ! ************************************************************************* 4223 4224 !Read and write some description of file from first line (character data) 4225 read (funit,'(a)') title 4226 write(msg, '(a,a)' ) '- ',trim(title) 4227 call wrtout(ab_out,msg,'COLL') 4228 call wrtout(std_out, msg,'COLL') 4229 4230 !Read and write more data describing psp parameters 4231 read (funit,*) znucl,zion,pspdat 4232 write(msg, '(a,f9.5,f10.5,2x,i8,t47,a)' ) & 4233 & '-',znucl,zion,pspdat,'znucl, zion, pspdat' 4234 call wrtout(ab_out,msg,'COLL') 4235 call wrtout(std_out, msg,'COLL') 4236 4237 read (funit,*) pspcod,pspxc,lmax,lloc,mmax,r2well 4238 if(pspxc<0) then 4239 write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) & 4240 & pspcod,pspxc,lmax,lloc,mmax,r2well,& 4241 & 'pspcod,pspxc,lmax,lloc,mmax,r2well' 4242 else 4243 write(msg, '(4i5,i10,f10.5,t47,a)' ) & 4244 & pspcod,pspxc,lmax,lloc,mmax,r2well,& 4245 & 'pspcod,pspxc,lmax,lloc,mmax,r2well' 4246 end if 4247 call wrtout(ab_out,msg,'COLL') 4248 call wrtout(std_out, msg,'COLL') 4249 4250 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
SOURCE
4274 subroutine pawpsp_read_header_2(funit,pspversion,basis_size,lmn_size) 4275 4276 !Arguments ------------------------------------ 4277 !scalars 4278 integer,intent(in):: funit 4279 integer,intent(out) :: pspversion,basis_size,lmn_size 4280 4281 !Local variables------------------------------- 4282 integer :: creatorid 4283 character(len=80) :: pspline 4284 character(len=500) :: msg 4285 4286 ! ************************************************************************* 4287 4288 !Read psp version in line 4 of the header 4289 pspversion=1 4290 read (funit,'(a80)') pspline;pspline=adjustl(pspline) 4291 if (pspline(1:3)=="paw".or.pspline(1:3)=="PAW") & 4292 & read(unit=pspline(4:80),fmt=*) pspversion 4293 if (pspversion<1.or.pspversion>5) then 4294 write(msg, '(a,i2,a,a,a)' )& 4295 & 'This version of PAW psp file (',pspversion,') is not compatible with',ch10,& 4296 & 'current version of Abinit.' 4297 LIBPAW_ERROR(msg) 4298 end if 4299 4300 if (pspversion==1) then 4301 read (unit=pspline,fmt=*) basis_size,lmn_size 4302 else 4303 ! Here psp file for Abinit 4.3+ 4304 read (unit=pspline(5:80),fmt=*) creatorid 4305 read (funit,*) basis_size,lmn_size 4306 end if 4307 4308 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
SOURCE
4442 subroutine pawpsp_read_header_xml(lloc,lmax,pspcod,pspxc,& 4443 & psxml,r2well,zion,znucl) 4444 4445 !Arguments ------------------------------------ 4446 !scalars 4447 type(paw_setup_t),intent(in) :: psxml 4448 integer,intent(out):: lloc,lmax,pspcod,pspxc 4449 real(dp),intent(out):: r2well,zion,znucl 4450 !Local variables------------------------------- 4451 integer :: il 4452 #if defined LIBPAW_HAVE_LIBXC 4453 integer :: ii,id 4454 #endif 4455 character(len=100) :: xclibxc 4456 character(len=500) :: msg 4457 !arrays 4458 4459 ! ************************************************************************* 4460 4461 lloc = 0 4462 r2well = 0 4463 pspcod=17 4464 znucl=psxml%atom%znucl 4465 zion =psxml%atom%zval 4466 4467 !lmax: 4468 lmax = 0 4469 do il=1,psxml%valence_states%nval 4470 if(psxml%valence_states%state(il)%ll>lmax) lmax=psxml%valence_states%state(il)%ll 4471 end do 4472 !pspxc 4473 select case(trim(psxml%xc_functional%name)) 4474 case('PZ') 4475 pspxc = 2 4476 #if defined LIBPAW_HAVE_LIBXC 4477 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4478 & +libxc_functionals_getid('XC_LDA_C_PZ')) 4479 #endif 4480 case('W') 4481 pspxc = 4 4482 #if defined LIBPAW_HAVE_LIBXC 4483 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4484 & +libxc_functionals_getid('XC_LDA_C_WIGNER')) 4485 #endif 4486 case('HL') 4487 pspxc = 5 4488 #if defined LIBPAW_HAVE_LIBXC 4489 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4490 & +libxc_functionals_getid('XC_LDA_C_HL')) 4491 #endif 4492 case('GL') 4493 #if defined LIBPAW_HAVE_LIBXC 4494 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4495 & +libxc_functionals_getid('XC_LDA_C_GL')) 4496 #else 4497 write(msg, '(7a)' )& 4498 & 'The exchange and correlation functional by Gunnarson-Lundqvist', ch10,& 4499 & 'is not implemented in Abinit.',ch10,& 4500 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4501 & ' generation or compile ABINIT with the libXC library.' 4502 LIBPAW_ERROR(msg) 4503 #endif 4504 case('VWN') 4505 #if defined LIBPAW_HAVE_LIBXC 4506 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4507 & +libxc_functionals_getid('XC_LDA_C_VWN')) 4508 #else 4509 write(msg, '(7a)' )& 4510 & 'The exchange and correlation functional by Vosko,Wilk and Nusair', ch10,& 4511 & 'is not implemented in Abinit.',ch10,& 4512 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4513 & ' generation or compile ABINIT with the libXC library.' 4514 LIBPAW_ERROR(msg) 4515 #endif 4516 case('PW') 4517 pspxc = 7 4518 #if defined LIBPAW_HAVE_LIBXC 4519 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4520 & +libxc_functionals_getid('XC_LDA_C_PW')) 4521 #endif 4522 case('PBE') 4523 pspxc = 11 4524 #if defined LIBPAW_HAVE_LIBXC 4525 pspxc = -(libxc_functionals_getid('XC_GGA_X_PBE')*1000 & 4526 & +libxc_functionals_getid('XC_GGA_C_PBE')) 4527 #endif 4528 case('revPBE') 4529 pspxc = 14 4530 #if defined LIBPAW_HAVE_LIBXC 4531 pspxc = -(libxc_functionals_getid('XC_GGA_X_PBE_R')*1000 & 4532 & +libxc_functionals_getid('XC_GGA_C_PBE')) 4533 #endif 4534 case('RPBE') 4535 pspxc = 15 4536 #if defined LIBPAW_HAVE_LIBXC 4537 pspxc = -(libxc_functionals_getid('XC_GGA_X_RPBE')*1000 & 4538 & +libxc_functionals_getid('XC_GGA_C_PBE')) 4539 #endif 4540 case('PW91') 4541 #if defined LIBPAW_HAVE_LIBXC 4542 pspxc = -(libxc_functionals_getid('XC_GGA_X_PW91')*1000 & 4543 & +libxc_functionals_getid('XC_GGA_C_PW91')) 4544 #else 4545 write(msg, '(7a)' )& 4546 & 'The exchange and correlation functional by Perdew and Wang 91', ch10,& 4547 & 'is not implemented in Abinit.',ch10,& 4548 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4549 & ' generation or compile ABINIT with the libXC library.' 4550 LIBPAW_ERROR(msg) 4551 #endif 4552 case('BLYP') 4553 #if defined LIBPAW_HAVE_LIBXC 4554 pspxc = -(libxc_functionals_getid('XC_GGA_X_B88')*1000 & 4555 & +libxc_functionals_getid('XC_GGA_C_LYP')) 4556 #else 4557 write(msg, '(7a)' )& 4558 & 'The exchange and correlation functional BLYP', ch10,& 4559 & 'is not implemented in Abinit.',ch10,& 4560 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4561 & ' generation or compile ABINIT with the libXC library.' 4562 LIBPAW_ERROR(msg) 4563 #endif 4564 case DEFAULT 4565 xclibxc=trim(psxml%xc_functional%name) 4566 if (xclibxc(1:3)=='XC_' .or.xclibxc(1:3)=='xc_' .or. & 4567 & xclibxc(1:5)=='LDA_X'.or.xclibxc(1:5)=='LDA_C'.or. & 4568 & xclibxc(1:5)=='lda_x'.or.xclibxc(1:5)=='lda_c'.or. & 4569 & xclibxc(1:5)=='GGA_X'.or.xclibxc(1:5)=='GGA_C'.or. & 4570 & xclibxc(1:5)=='gga_x'.or.xclibxc(1:5)=='gga_c'.or. & 4571 & xclibxc(1:6)=='MGGA_X'.or.xclibxc(1:6)=='MGGA_C'.or. & 4572 & xclibxc(1:6)=='mgga_x'.or.xclibxc(1:6)=='mgga_c') then 4573 #if defined LIBPAW_HAVE_LIBXC 4574 pspxc=0 4575 ii=index(xclibxc,'+') ; if (ii<=0) ii=0 4576 if (ii>0) then 4577 id=libxc_functionals_getid(xclibxc(1:ii-1)) 4578 if (id<=0) then 4579 write(msg, '(3a)' ) 'The ',xclibxc(1:ii-1), & 4580 & ' functional (read from PAW-XML file) was not found in the libXC library!' 4581 LIBPAW_ERROR(msg) 4582 end if 4583 pspxc=pspxc-id*1000 4584 end if 4585 id=libxc_functionals_getid(xclibxc(ii+1:)) 4586 if (id<=0) then 4587 write(msg, '(3a)' ) 'The ',xclibxc(ii+1:), & 4588 & ' functional (read from PAW-XML file) was not found in the libXC library!' 4589 LIBPAW_ERROR(msg) 4590 end if 4591 pspxc=pspxc-id 4592 #else 4593 msg='Cannot use LibXC functional because ABINIT is not compiled with LibXC !' 4594 LIBPAW_ERROR(msg) 4595 #endif 4596 ! To be eliminated later (temporary) 4597 else if(trim(psxml%xc_functional%functionaltype)=='LIBXC')then 4598 #if defined LIBPAW_HAVE_LIBXC 4599 xclibxc=trim(psxml%xc_functional%name) 4600 read(unit=xclibxc,fmt=*) pspxc 4601 pspxc=-pspxc 4602 #else 4603 msg='Cannot use LibXC functional because ABINIT is not compiled with LibXC !' 4604 LIBPAW_ERROR(msg) 4605 #endif 4606 else 4607 write(msg, '(3a)') 'Unknown XC functional in psp file: ',trim(xclibxc),' !' 4608 LIBPAW_ERROR(msg) 4609 end if 4610 end select 4611 4612 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
SOURCE
4633 subroutine pawpsp_read_pawheader(basis_size,lmax,lmn_size,& 4634 & l_size,mesh_size,pspversion,psxml,rpaw,rshp,shape_type) 4635 4636 !Arguments ------------------------------------ 4637 !scalars 4638 integer,intent(in):: lmax 4639 integer,intent(out):: basis_size,mesh_size,lmn_size,l_size 4640 integer,intent(out):: pspversion,shape_type 4641 real(dp),intent(out)::rpaw,rshp 4642 type(paw_setup_t),intent(in) :: psxml 4643 !Local variables------------------------------- 4644 integer::il 4645 4646 ! ************************************************************************* 4647 4648 !All of this was moved from pawpsxml2ab, 4649 !basis_size 4650 basis_size=psxml%valence_states%nval 4651 !mesh_size 4652 do il=1,psxml%ngrid 4653 if(psxml%radial_grid(il)%id==psxml%idgrid) & 4654 & mesh_size=psxml%radial_grid(il)%iend-psxml%radial_grid(il)%istart+1 4655 end do 4656 !lmn_size: 4657 lmn_size=0 4658 do il=1,psxml%valence_states%nval 4659 lmn_size=lmn_size+2*psxml%valence_states%state(il)%ll+1 4660 end do 4661 !lsize 4662 l_size=2*lmax+1 4663 !pspversion 4664 pspversion=10 4665 !rpaw: 4666 rpaw=0.d0 4667 if (psxml%rpaw<0.d0) then 4668 do il=1,psxml%valence_states%nval 4669 if(psxml%valence_states%state(il)%rc>rpaw) rpaw=psxml%valence_states%state(il)%rc 4670 end do 4671 else 4672 rpaw=psxml%rpaw 4673 end if 4674 !shape_type, rshp: 4675 select case(trim(psxml%shape_function%gtype)) 4676 case('gauss') 4677 shape_type=1 4678 rshp=rpaw 4679 case('bessel') 4680 shape_type=3 4681 rshp=psxml%shape_function%rc 4682 case('sinc') 4683 shape_type=2 4684 rshp=psxml%shape_function%rc 4685 case('exp') 4686 shape_type=1 4687 rshp=rpaw 4688 case('num') 4689 shape_type=-1 4690 rshp=rpaw 4691 end select 4692 4693 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
SOURCE
1652 subroutine pawpsp_rw_atompaw(basis_size,filpsp,wvl) 1653 1654 !Arguments ------------------------------------ 1655 integer,intent(in):: basis_size 1656 type(wvlpaw_type),intent(in)::wvl 1657 character(len=fnlen),intent(in)::filpsp 1658 !arrays 1659 character(strlen) :: pspline 1660 character(len=fnlen)::fname 1661 1662 !Local variables------------------------------- 1663 integer :: ib,ii,ios,jj,step,iunt,ount 1664 !arrays 1665 1666 ! ************************************************************************* 1667 iunt = libpaw_get_free_unit() 1668 ount = libpaw_get_free_unit() 1669 1670 step=0 1671 ! Open psp file for reading 1672 open(unit=iunt,file=trim(filpsp),form='formatted',status='old',action="read") 1673 ! Open the file for writing 1674 write(fname,'(2a)') libpaw_basename(trim(filpsp)),".wvl" 1675 open(unit=ount,file=fname,form='formatted',status='unknown',action="write") 1676 1677 read_loop: do 1678 if(step==0) then 1679 read(iunt,'(a)',IOSTAT=ios) pspline 1680 if ( ios /= 0 ) exit read_loop 1681 if(index(trim(pspline),'CORE_DENSITY')/=0 .and. & 1682 & index(trim(pspline),'PSEUDO_CORE_DENSITY')==0 .and. & 1683 & index(trim(pspline),'TCORE_DENSITY')==0 ) then 1684 step=1 1685 else 1686 write(ount,'(a)') trim(pspline) 1687 end if 1688 elseif(step==1) then 1689 !Write Gaussian projectors: 1690 jj=0 1691 do ib=1,basis_size 1692 write(ount,'(a,i1,a)') "===== GAUSSIAN_TPROJECTOR ",ib,& 1693 & " ===== " 1694 write(ount,'(i5,1x,i5,1x,a)')wvl%pngau(ib),wvl%ptotgau, ":ngauss, total ngauss" 1695 write(ount,'(3(1x,es23.16))')(wvl%parg(:,ii),& 1696 & ii=jj+1,jj+wvl%pngau(ib)) 1697 write(ount,'(3(1x,es23.16))')(wvl%pfac(:,ii),& 1698 & ii=jj+1,jj+wvl%pngau(ib)) 1699 jj=jj+wvl%pngau(ib) 1700 end do 1701 write(ount,'(a)') trim(pspline) 1702 step=0 1703 end if 1704 end do read_loop 1705 1706 close(iunt) 1707 close(ount) 1708 1709 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
SOURCE
2790 subroutine pawpsp_vhar2rho(radmesh,rho,vv) 2791 2792 !Arguments ------------------------------------ 2793 type(pawrad_type),intent(in) :: radmesh 2794 real(dp), intent(in) :: vv(:) 2795 real(dp), intent(out):: rho(:) 2796 2797 !Local variables------------------------------- 2798 integer :: nr 2799 real(dp) :: dfdr(radmesh%mesh_size),d2fdr(radmesh%mesh_size) 2800 2801 ! ************************************************************************* 2802 2803 nr=size(vv) 2804 if (nr/=size(rho)) then 2805 LIBPAW_BUG('wrong sizes!') 2806 end if 2807 2808 !Laplacian = 2809 !\frac{\partial^2}{\partial r^2} + 2/r \frac{\partial}{\partial r} 2810 2811 !Calculate derivatives 2812 call nderiv_gen(dfdr(1:nr),vv,radmesh,der2=d2fdr(1:nr)) 2813 2814 rho(2:nr)=d2fdr(2:nr) + 2._dp*dfdr(2:nr)/radmesh%rad(2:nr) 2815 call pawrad_deducer0(rho,nr,radmesh) 2816 2817 rho(1:nr)=-rho(1:nr)/(4._dp*pi) 2818 2819 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
SOURCE
4331 subroutine pawpsp_wvl(filpsp,pawrad, pawtab,usewvl, wvl_ngauss, comm_mpi) 4332 4333 !Arguments------------------------------------ 4334 !scalars 4335 integer, optional,intent(in):: comm_mpi 4336 integer, intent(in):: usewvl, wvl_ngauss(2) 4337 character(len=fnlen),intent(in)::filpsp 4338 type(pawrad_type),intent(in) :: pawrad 4339 type(pawtab_type),intent(inout):: pawtab 4340 !arrays 4341 4342 !Local variables------------------------------- 4343 !scalars 4344 integer:: ii, me, mparam, nterm_bounds(2) 4345 type(pawrad_type)::tproj_mesh 4346 character(len=500) :: msg 4347 !arrays 4348 integer,allocatable:: ngauss_param(:) 4349 real(dp),allocatable:: gauss_param(:,:) 4350 4351 ! ************************************************************************* 4352 4353 me=0; if (present(comm_mpi))me=xmpi_comm_rank(comm_mpi) 4354 4355 !If usewvl flag is on, we must have the pawtab%wvl pointer allocated 4356 if (usewvl==1.and.pawtab%has_wvl==0) then 4357 call wvlpaw_allocate(pawtab%wvl) 4358 pawtab%has_wvl=1 4359 end if 4360 4361 !Fit projectors to a sum of Gaussians: 4362 if (usewvl ==1 .and. pawtab%wvl%ptotgau==0 ) then 4363 4364 if (pawtab%has_tproj==0) then 4365 msg='pawtab%tproj must be allocated' 4366 LIBPAW_BUG(msg) 4367 end if 4368 4369 ! 1) fit projectors to gaussians 4370 write(msg,'(a,a)')ch10,'Fitting tproj to Gaussians' 4371 call wrtout(std_out,msg,'COLL') 4372 4373 ! See fit_gen (option==4): 4374 do ii=1,2 4375 nterm_bounds(ii)=ceiling(wvl_ngauss(ii)/2.0) 4376 end do 4377 mparam=nterm_bounds(2)*4 4378 LIBPAW_ALLOCATE(gauss_param,(mparam,pawtab%basis_size)) 4379 LIBPAW_ALLOCATE(ngauss_param,(pawtab%basis_size)) 4380 ! compute tproj_mesh 4381 call pawrad_init(tproj_mesh,mesh_size=size(pawtab%tproj,1),& 4382 & mesh_type=pawrad%mesh_type,rstep=pawrad%rstep, lstep=pawrad%lstep) 4383 4384 if(present(comm_mpi)) then 4385 call gaussfit_projector(pawtab%basis_size,mparam,& 4386 & ngauss_param,nterm_bounds,pawtab%orbitals,& 4387 & gauss_param,tproj_mesh,& 4388 & pawtab%rpaw,pawtab%tproj,comm_mpi) 4389 else 4390 call gaussfit_projector(pawtab%basis_size,mparam,& 4391 & ngauss_param,nterm_bounds,pawtab%orbitals,& 4392 & gauss_param,tproj_mesh,& 4393 & pawtab%rpaw,pawtab%tproj) 4394 end if 4395 ! tproj is now as a sum of sin+cos functions, 4396 ! convert it to a sum of complex gaussians and fill %wvl object: 4397 call pawpsp_wvl_sin2gauss(pawtab%basis_size,mparam,& 4398 & ngauss_param,gauss_param,pawtab%wvl) 4399 LIBPAW_DEALLOCATE(gauss_param) 4400 LIBPAW_DEALLOCATE(ngauss_param) 4401 4402 if(me==0) then 4403 call pawpsp_rw_atompaw(pawtab%basis_size,filpsp,pawtab%wvl) 4404 end if 4405 4406 pawtab%has_wvl=2 4407 4408 end if 4409 4410 !Projectors in real space are no more needed 4411 call pawrad_free(tproj_mesh) 4412 if(allocated(pawtab%tproj)) then 4413 LIBPAW_DEALLOCATE(pawtab%tproj) 4414 pawtab%has_tproj=0 4415 end if 4416 4417 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
SOURCE
2848 subroutine pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr) 2849 2850 !Arguments ------------------------------------ 2851 !scalars 2852 integer,intent(in)::usewvl 2853 type(pawrad_type),intent(in) :: vale_mesh 2854 type(pawtab_type),intent(inout) :: pawtab 2855 type(pawrad_type),intent(in) ::vloc_mesh 2856 2857 !arrays 2858 real(dp),intent(in) :: tnvale(vale_mesh%mesh_size*pawtab%has_tvale) 2859 real(dp),intent(in) :: vlocr(vloc_mesh%mesh_size) 2860 2861 2862 !Local variables ------------------------------ 2863 !scalars 2864 integer :: msz 2865 character(len=500) :: msg 2866 !arrays 2867 2868 ! ************************************************************************* 2869 2870 !If usewvl flag is on, we must have the pawtab%wvl pointer allocated 2871 if (pawtab%has_wvl==0) then 2872 msg='pawtab%has_wvl flag should be on o entry' 2873 LIBPAW_BUG(msg) 2874 end if 2875 call wvlpaw_allocate(pawtab%wvl) 2876 2877 !========================================================== 2878 !Change mesh_size of tvalespl 2879 !Compute second derivative from tNvale(r) 2880 2881 if (pawtab%has_tvale/=0) then 2882 if(usewvl==1) then 2883 if(allocated(pawtab%tvalespl)) then 2884 LIBPAW_DEALLOCATE(pawtab%tvalespl) 2885 end if 2886 LIBPAW_ALLOCATE(pawtab%tvalespl,(vale_mesh%mesh_size,2)) 2887 pawtab%tnvale_mesh_size=vale_mesh%mesh_size 2888 pawtab%tvalespl(:,1)=tnvale 2889 ! Compute second derivative of tvalespl(r) 2890 call paw_spline(vale_mesh%rad,pawtab%tvalespl(:,1),vale_mesh%mesh_size,zero,zero,pawtab%tvalespl(:,2)) 2891 end if 2892 else 2893 pawtab%dnvdq0=zero 2894 pawtab%tnvale_mesh_size=0 2895 end if 2896 2897 !========================================================== 2898 !Save rholoc: 2899 !Get local density from local potential 2900 !use the poisson eq. 2901 msz=vloc_mesh%mesh_size 2902 call wvlpaw_rholoc_free(pawtab%wvl%rholoc) 2903 LIBPAW_ALLOCATE(pawtab%wvl%rholoc%d,(msz,4)) 2904 LIBPAW_ALLOCATE(pawtab%wvl%rholoc%rad,(msz)) 2905 pawtab%wvl%rholoc%msz=msz 2906 pawtab%wvl%rholoc%rad(1:msz)=vloc_mesh%rad(1:msz) 2907 2908 !get rho from v: 2909 call pawpsp_vhar2rho(vloc_mesh,pawtab%wvl%rholoc%d(:,1),vlocr) 2910 ! 2911 !get second derivative, and store it 2912 call paw_spline(pawtab%wvl%rholoc%rad,pawtab%wvl%rholoc%d(:,1),pawtab%wvl%rholoc%msz,& 2913 & zero,zero,pawtab%wvl%rholoc%d(:,2)) 2914 2915 !save also vlocr: 2916 pawtab%wvl%rholoc%d(:,3)=vlocr 2917 2918 !get second derivative, and store it 2919 call paw_spline(pawtab%wvl%rholoc%rad,vlocr,pawtab%wvl%rholoc%msz,& 2920 & zero,zero,pawtab%wvl%rholoc%d(:,4)) 2921 2922 !Test 2923 !do ii=1,pawtab%wvl%rholoc%msz 2924 !write(503,'(3(f16.10,x))')pawtab%wvl%rholoc%rad(ii),pawtab%wvl%rholoc%d(ii,1),pawtab%wvl%rholoc%d(ii,3) 2925 !end do 2926 ! 2927 !Do splint 2928 ! 2929 !nmesh=4000 2930 !rread1= (9.9979999d0/real(nmesh-1,dp)) ! 0.0025001d0 !step 2931 !allocate(raux1(nmesh),raux2(nmesh)) 2932 !do ii=1,nmesh 2933 !raux1(ii)=rread1*real(ii-1,dp) !mesh 2934 !end do 2935 !call splint(pawtab%wvl%rholoc%msz,pawtab%wvl%rholoc%rad,pawtab%wvl%rholoc%d(:,1),pawtab%wvl%rholoc%d(:,2),& 2936 !& nmesh,raux1,raux2,ierr) 2937 !do ii=1,nmesh 2938 !write(401,'(10(f20.7,x))')raux1(ii),raux2(ii),raux2(ii)*raux1(ii)**2 2939 !end do 2940 !deallocate(raux1,raux2) 2941 2942 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
SOURCE
4063 subroutine pawpsp_wvl_sin2gauss(basis_size,mparam,nparam,& 4064 & param,wvl) 4065 4066 !Arguments ------------------------------------ 4067 integer,intent(in) :: mparam,basis_size 4068 integer,intent(in) :: nparam(basis_size) 4069 real(dp),intent(in) :: param(mparam,basis_size) 4070 type(wvlpaw_type),intent(inout):: wvl 4071 4072 !Local variables ------------------------------ 4073 integer :: i,ii,ib,ngauss,nterm 4074 real(dp) :: sep 4075 real(dp) :: a1(mparam),a2(mparam),a3(mparam),a4(mparam),a5(mparam) 4076 real(dp) :: b1r(mparam),b2r(mparam),b1i(mparam),b2i(mparam) 4077 character(len=500) :: message 4078 ! 4079 !extra variables, use to debug 4080 ! 4081 !integer::igau,nr,unitp 4082 !real(dp)::step,rmax 4083 !real(dp),allocatable::r(:), y(:) 4084 !complex::fac,arg 4085 !complex(dp),allocatable::f(:) 4086 !************************************************************************ 4087 4088 ! Convert from \sum(sin+cos) expressions to sums of complex gaussians 4089 ! (only works for option=4, see fit_gen) 4090 4091 ! get number of coefficients: 4092 ii=0 4093 do ib=1,basis_size 4094 nterm=nparam(ib)/4 !option=4, there are 4 parameters for each term 4095 ii=ii+nterm*2 !two gaussians for each term 4096 end do 4097 ! 4098 ! Allocate objects 4099 ! 4100 ngauss=ii 4101 wvl%ptotgau=ngauss !total number of complex gaussians 4102 LIBPAW_ALLOCATE(wvl%pfac,(2,ngauss)) 4103 LIBPAW_ALLOCATE(wvl%parg,(2,ngauss)) 4104 LIBPAW_ALLOCATE(wvl%pngau,(basis_size)) 4105 wvl%pngau(1:basis_size)=nparam(1:basis_size)/2 !option=4 4106 ! 4107 ii=0 4108 do ib=1,basis_size 4109 ! 4110 ! Get parameters in sin+cos expansion: 4111 ! Option4: \sum a1 exp(-a2 x^2) ( a3 sin(k x^2) + a4 cos(k x^2)) 4112 ! 4113 nterm=nparam(ib)/4 !option=4 4114 ! 4115 a1(1:nterm)=param(1:nterm,ib) 4116 a2(1:nterm)=param(nterm+1:nterm*2,ib) 4117 a3(1:nterm)=param(nterm*2+1:nterm*3,ib) 4118 a4(1:nterm)=param(nterm*3+1:nterm*4,ib) 4119 sep=1.1d0 4120 do i=1,nterm 4121 a5(i)=sep**(i) 4122 end do 4123 4124 ! First check that "a2" is a positive number (it is multiplied by -1, so 4125 ! that gaussians decay to zero: 4126 if( any(a2(1:nterm) < tol12) ) then 4127 message = 'Real part of Gaussians should be a negative number (they should go to zero at infty)' 4128 LIBPAW_ERROR(message) 4129 end if 4130 4131 ! 4132 ! Now translate them to a sum of complex gaussians: 4133 ! pngau(ib)=nterm*2 4134 ! Two gaussians by term: 4135 ! 4136 ! First gaussian 4137 b1r(1:nterm)= a1(1:nterm)*a4(1:nterm)/2.d0 !coefficient, real 4138 b1i(1:nterm)=-a1(1:nterm)*a3(1:nterm)/2.d0 !coefficient, imag 4139 b2r(1:nterm)=-a2(1:nterm) !exponential, real 4140 b2i(1:nterm)= a5(1:nterm) !exponential, imag 4141 ! 4142 wvl%pfac(1,ii+1:ii+nterm)=b1r(1:nterm) 4143 wvl%pfac(2,ii+1:ii+nterm)=b1i(1:nterm) 4144 wvl%parg(1,ii+1:ii+nterm)=b2r(1:nterm) 4145 wvl%parg(2,ii+1:ii+nterm)=b2i(1:nterm) 4146 ! Second gaussian 4147 wvl%pfac(1,ii+nterm+1:ii+nterm*2)= b1r(1:nterm) 4148 wvl%pfac(2,ii+nterm+1:ii+nterm*2)=-b1i(1:nterm) 4149 wvl%parg(1,ii+nterm+1:ii+nterm*2)= b2r(1:nterm) 4150 wvl%parg(2,ii+nterm+1:ii+nterm*2)=-b2i(1:nterm) 4151 ! 4152 ii=ii+nterm*2 4153 end do 4154 4155 ! begin debug 4156 ! write(*,*)'pawpsp_wvl_sin2gauss, comment me' 4157 ! nr=3000 4158 ! rmax=10.d0 4159 ! LIBPAW_ALLOCATE(r,(nr)) 4160 ! LIBPAW_ALLOCATE(f,(nr)) 4161 ! LIBPAW_ALLOCATE(y,(nr)) 4162 ! step=rmax/real(nr-1,dp) 4163 ! do ir=1,nr 4164 ! r(ir)=real(ir-1,dp)*step 4165 ! end do 4166 ! ! 4167 ! ii=0 4168 ! do ib=1,basis_size 4169 ! unitp=500+ib 4170 ! f(:)=czero 4171 ! ! 4172 ! do igau=1,wvl%pngau(ib) 4173 ! ii=ii+1 4174 ! arg=cmplx(wvl%parg(1,ii),wvl%parg(2,ii)) 4175 ! fac=cmplx(wvl%pfac(1,ii),wvl%pfac(2,ii)) 4176 ! f(:)=f(:)+fac*exp(arg*r(:)**2) 4177 ! end do 4178 ! do ir=1,nr 4179 ! write(unitp,'(3f16.7)')r(ir),real(f(ir))!,y(ir) 4180 ! end do 4181 ! end do 4182 ! LIBPAW_DEALLOCATE(r) 4183 ! LIBPAW_DEALLOCATE(f) 4184 ! LIBPAW_DEALLOCATE(y) 4185 ! end debug 4186 4187 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
SOURCE
5007 subroutine pawpsp_check_xml_upf(filpsp) 5008 5009 !Arguments ------------------------------------ 5010 !scalars 5011 character(len=fnlen),intent(in):: filpsp ! name of the psp file 5012 5013 !Local variables------------------------------- 5014 integer :: unt 5015 character(len=70):: testxml 5016 5017 ! ************************************************************************* 5018 5019 ! Check if the file pseudopotential file is written in XML 5020 usexml = 0 5021 unt = libpaw_get_free_unit() 5022 open (unit=unt,file=filpsp,form='formatted',status='old',action="read") 5023 rewind (unit=unt) 5024 read(unt,*) testxml 5025 if(testxml(1:5)=='<?xml')then 5026 usexml = 1 5027 read(unt,*) testxml 5028 if(testxml(1:4)/='<paw')then 5029 msg='Reading a NC pseudopotential for a PAW calculation?' 5030 LIBPAW_BUG(msg) 5031 end if 5032 else 5033 usexml = 0 5034 end if 5035 close (unit=unt) 5036 5037 ! Check if pseudopotential file is a Q-espresso UPF file 5038 unt = libpaw_get_free_unit() 5039 open (unit=unt,file=filpsp,form='formatted',status='old',action="read") 5040 rewind (unit=unt) 5041 read(unt,*) testxml ! just a string, no relation to xml. 5042 if(testxml(1:9)=='<PP_INFO>')then 5043 msg='UPF format not allowed with PAW (USPP part not read yet)!' 5044 LIBPAW_ERROR(msg) 5045 end if 5046 close (unit=unt) 5047 5048 end subroutine pawpsp_check_xml_upf
pawpsp_main/pawpsp_consistency [ Functions ]
[ Top ] [ pawpsp_main ] [ Functions ]
NAME
pawpsp_consistency
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
5070 subroutine pawpsp_consistency() 5071 5072 ! ************************************************************************* 5073 5074 !Check pspcod=7 or 17 5075 if(pspcod/=7 .and. pspcod/=17)then 5076 write(msg, '(a,i2,a,a)' )& 5077 & 'In reading atomic psp file, finds pspcod=',pspcod,ch10,& 5078 & 'This is not an allowed value within PAW.' 5079 LIBPAW_BUG(msg) 5080 end if 5081 5082 !Does nuclear charge znuclpsp agree with psp input znucl 5083 if (abs(znuclpsp-znucl)>tol8) then 5084 write(msg, '(a,f10.5,2a,f10.5,5a)' )& 5085 & 'Pseudopotential file znucl=',znucl,ch10,& 5086 & 'does not equal input znuclpsp=',znuclpsp,' better than 1e-08 .',ch10,& 5087 & 'znucl is read from the psp file in pspatm_abinit, while',ch10,& 5088 & 'znuclpsp is read in iofn2.' 5089 LIBPAW_BUG(msg) 5090 end if 5091 5092 !Does nuclear charge zionpsp agree with psp input zion 5093 if (abs(zionpsp-zion)>tol8) then 5094 write(msg, '(a,f10.5,2a,f10.5,5a)' )& 5095 & 'Pseudopotential file zion=',zion,ch10,& 5096 & 'does not equal input zionpsp=',zionpsp,' better than 1e-08 .',ch10,& 5097 & 'zion is read from the psp file in pawpsp_main, while',ch10,& 5098 & 'zionpsp is read in iofn2.' 5099 LIBPAW_BUG(msg) 5100 end if 5101 5102 !Check several choices for ixc against pspxc 5103 !ixc is from ABINIT code; pspxc is from atomic psp file 5104 if (ixc==0) then 5105 msg='Note that input ixc=0 => no xc is being used.' 5106 LIBPAW_WARNING(msg) 5107 else if(ixc/=pspxc) then 5108 write(msg, '(a,i8,a,a,a,i8,a,a,a,a,a,a,a,a,a,a)' ) & 5109 & 'Pseudopotential file pspxc=',pspxc,',',ch10,& 5110 & 'not equal to input ixc=',ixc,'.',ch10,& 5111 & 'These parameters must agree to get the same xc ',ch10,& 5112 & 'in ABINIT code as in psp construction.',ch10,& 5113 & 'Action : check psp design or input file.',ch10,& 5114 & 'Assume experienced user. Execution will continue.',ch10 5115 LIBPAW_WARNING(msg) 5116 end if 5117 5118 if (lloc>lmax ) then 5119 write(msg, '(a,2i12,a,a,a,a)' )& 5120 & 'lloc,lmax=',lloc,lmax,ch10,& 5121 & 'chosen l of local psp exceeds range from input data.',ch10,& 5122 & 'Action : check pseudopotential input file.' 5123 LIBPAW_ERROR(msg) 5124 end if 5125 5126 end subroutine pawpsp_consistency