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-2018 ABINIT group (MT, FJ,TR, GJ, FB, FrD, AF, GMR, DRH) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
FOR DEVELOPERS: in order to preserve the portability of libPAW library, please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
SOURCE
21 #include "libpaw.h" 22 23 module m_pawpsp 24 25 USE_DEFS 26 USE_MSG_HANDLING 27 USE_MPI_WRAPPERS 28 USE_MEMORY_PROFILING 29 30 use m_libpaw_libxc 31 #if defined LIBPAW_HAVE_FOX 32 use fox_sax 33 #endif 34 35 use m_libpaw_tools, only : libpaw_basename, libpaw_get_free_unit 36 37 use m_pawang, only: pawang_type 38 use m_pawtab, only: pawtab_type, wvlpaw_type, wvlpaw_allocate, wvlpaw_rholoc_free, & 39 & pawtab_free, wvlpaw_free, wvlpaw_rholoc_nullify, pawtab_bcast, & 40 & pawtab_set_flags, wvlpaw_allocate, wvlpaw_free, wvlpaw_rholoc_nullify, & 41 & wvlpaw_rholoc_free 42 use m_pawxmlps, only: rdpawpsxml_core, paw_setup_t, paw_setuploc, paw_setup_free 43 use m_pawrad, only: pawrad_type, pawrad_init, pawrad_free, pawrad_copy, & 44 & pawrad_bcast, pawrad_ifromr, simp_gen, nderiv_gen, bound_deriv, pawrad_deducer0, poisson 45 use m_paw_numeric, only: paw_splint, paw_spline, paw_smooth, paw_jbessel_4spline 46 use m_paw_atom, only: atompaw_shapebes, atompaw_vhnzc, atompaw_shpfun, & 47 & atompaw_dij0, atompaw_kij 48 use m_pawxc, only: pawxc, pawxcm 49 use m_paw_gaussfit, only: gaussfit_projector 50 51 implicit none 52 53 private 54 55 public:: pawpsp_calc_d5 !calculate up to the 5th derivative 56 public:: pawpsp_main !main routine to read psp 57 public:: pawpsp_nl !make paw projector form factors f_l(q) 58 public:: pawpsp_read !read psp from file 59 public:: pawpsp_read_header !read header of psp file 60 public:: pawpsp_read_corewf !read core wavefunction 61 public:: pawpsp_read_header_2 !reads pspversion, basis_size and lmn_size 62 public:: pawpsp_rw_atompaw !read and writes ATOMPAW psp with gaussian |p> 63 public:: pawpsp_wvl !wavelet and icoulomb>0 related operations 64 public:: pawpsp_wvl_calc !wavelet related operations 65 public:: pawpsp_7in !reads non-XML atomic data 66 public:: pawpsp_17in !reads XML atomic data 67 public:: pawpsp_calc !calculates atomic quantities from psp info 68 public:: pawpsp_read_header_xml !read header of psp file for XML 69 public:: pawpsp_read_pawheader !read header variables from XML objects 70 public:: pawpsp_bcast ! broadcast PAW psp data 71 public:: pawpsp_cg !compute sin FFT transform of a density 72 public:: pawpsp_lo !compute sin FFT transform of local potential 73 74 ! Private procedures 75 private:: pawpsp_wvl_sin2gauss !convert sin/cos to gaussians
m_pawpsp/pawpsp_17in [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_17in
FUNCTION
Initialize pspcod=17 ("PAW XML pseudopotentials"):
continue to read the corresponding file and compute the form factors
INPUTS
ipsp= id in the array of the currently read pseudo. ixc=exchange-correlation choice from main routine data file lmax=value of lmax mentioned at the second line of the psp file lnmax=max. number of (l,n) components over all type of psps angular momentum of nonlocal pseudopotential mmax=max number of pts in real space grid (already read in the psp file header) mqgrid_ff=dimension of q (or G) grid for nl form factors (array ffspl) mqgrid_vl=dimension of q (or G) grid for Vloc (array vlspl) pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) pspheads= header of the current pseudopotential qgrid_ff(psps%mqgrid_ff)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors qgrid_vl(psps%mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for Vloc xclevel= XC functional level zion=nominal valence of atom as specified in psp file znucl=atomic number of atom as specified in input file to main routine
OUTPUT
epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
ffspl(psps%mqgrid_ff,2,psps%lnmax)=form factor f_l(q) and second derivative
from spline fit for each angular momentum and each projector;
pawrad <type(pawrad_type)>=paw radial mesh and related data
pawtab <type(pawtab_type)>=paw tabulated starting data
vlspl(psps%mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
wvl_crmult,wvl_frmult= variables definining the fine and coarse grids in a wavelets calculation
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
xcccrc=XC core correction cutoff radius (bohr) from psp file
NOTES
Spin-orbit not yet implemented (to be done) Comments: * mesh_type= type of radial mesh mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n * radial shapefunction type shape_type=-1 ; gl(r)=numeric (read from psp file) shape_type= 1 ; gl(r)=k(r).r^l; k(r)=exp[-(r/sigma)**lambda] shape_type= 2 ; gl(r)=k(r).r^l; k(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2 if r<=rshp shape_type= 3 ; gl(r)=Alpha(1,l)*jl(q(1,l)*r)+Alpha(2,l)*jl(q(2,l)*r) for each l
PARENTS
m_pawpsp,pspatm
CHILDREN
SOURCE
2953 subroutine pawpsp_17in(epsatm,ffspl,icoulomb,ipsp,ixc,lmax,& 2954 & lnmax,mmax,mqgrid_ff,mqgrid_vl,pawpsp_header,pawrad,pawtab,& 2955 & pawxcdev, qgrid_ff,qgrid_vl,usewvl,usexcnhat_in,vlspl,xcccrc,& 2956 & xclevel,xc_denpos,zion,znucl) 2957 2958 2959 !This section has been created automatically by the script Abilint (TD). 2960 !Do not modify the following lines by hand. 2961 #undef ABI_FUNC 2962 #define ABI_FUNC 'pawpsp_17in' 2963 !End of the abilint section 2964 2965 implicit none 2966 2967 !Arguments ------------------------------------ 2968 !scalars 2969 integer,intent(in) :: ipsp,ixc,lmax,lnmax,mqgrid_ff,mqgrid_vl,pawxcdev,usexcnhat_in 2970 integer,intent(inout) ::mmax 2971 integer,intent(in) :: xclevel,icoulomb,usewvl 2972 real(dp),intent(in) :: xc_denpos,zion,znucl 2973 real(dp),intent(out) :: epsatm,xcccrc 2974 type(pawpsp_header_type),intent(in) :: pawpsp_header 2975 type(pawrad_type),intent(inout) :: pawrad 2976 type(pawtab_type),intent(inout) :: pawtab 2977 !arrays 2978 real(dp),intent(in) :: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl) 2979 real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax) 2980 real(dp),intent(out) :: vlspl(mqgrid_vl,2) 2981 2982 !Local variables ------------------------------ 2983 !scalars 2984 integer :: has_v_minushalf,ib,icoremesh,il,ilm,ilmn,ilmn0,iln,imainmesh,imsh,iprojmesh 2985 integer :: ir,iread1,ishpfmesh,ivalemesh,ivlocmesh,j0lmn,jlm,pngau 2986 integer :: jlmn,jln,klmn,msz,nmesh,nval,pspversion,shft,sz10,usexcnhat,vlocopt 2987 real(dp), parameter :: rmax_vloc=10.0_dp 2988 real(dp) :: fourpi,occ,rc,yp1,ypn 2989 logical :: save_core_msz 2990 character(len=500) :: msg 2991 type(pawrad_type) :: core_mesh,shpf_mesh,tproj_mesh,vale_mesh,vloc_mesh 2992 !arrays 2993 integer,allocatable :: mesh_shift(:),nprj(:) 2994 real(dp),allocatable :: kij(:),ncore(:) 2995 real(dp),allocatable :: shpf(:,:),tncore(:),tnvale(:),tproj(:,:),vhnzc(:),vlocr(:) 2996 real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:) 2997 type(pawrad_type),allocatable :: radmesh(:) 2998 2999 !************************************************************************ 3000 3001 !========================================================== 3002 !Destroy everything in pawtab but optional flags 3003 call pawtab_free(pawtab) 3004 !Destroy everything in pawrad 3005 call pawrad_free(pawrad) 3006 3007 !========================================================== 3008 !Initialize useful data 3009 3010 pawtab%usexcnhat=usexcnhat_in 3011 fourpi=4*acos(-1.d0) 3012 pspversion=pawpsp_header%pawver 3013 save_core_msz=(usewvl==1 .or. icoulomb .ne. 0) 3014 3015 !========================================================== 3016 !Initialize partial waves quantum numbers 3017 3018 pawtab%basis_size=pawpsp_header%basis_size 3019 LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size)) 3020 do ib=1,pawtab%basis_size 3021 pawtab%orbitals(ib)=paw_setuploc%valence_states%state(ib)%ll 3022 end do 3023 3024 !========================================================== 3025 !Initialize various dims and indexes 3026 3027 pawtab%lmn_size=pawpsp_header%lmn_size 3028 pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2 3029 pawtab%l_size=2*maxval(pawtab%orbitals)+1 3030 pawtab%ij_size=pawtab%basis_size*(pawtab%basis_size+1)/2 3031 3032 !indlmn calculation (indices for (l,m,n) basis) 3033 if (allocated(pawtab%indlmn)) then 3034 LIBPAW_DEALLOCATE(pawtab%indlmn) 3035 end if 3036 LIBPAW_ALLOCATE(pawtab%indlmn,(6,pawtab%lmn_size)) 3037 pawtab%indlmn(:,:)=0 3038 LIBPAW_BOUND1_ALLOCATE(nprj,BOUNDS(0,maxval(pawtab%orbitals))) 3039 ilmn=0;iln=0;nprj=0 3040 do ib=1,pawtab%basis_size 3041 il=pawtab%orbitals(ib) 3042 nprj(il)=nprj(il)+1 3043 iln=iln+1 3044 do ilm=1,2*il+1 3045 pawtab%indlmn(1,ilmn+ilm)=il 3046 pawtab%indlmn(2,ilmn+ilm)=ilm-(il+1) 3047 pawtab%indlmn(3,ilmn+ilm)=nprj(il) 3048 pawtab%indlmn(4,ilmn+ilm)=il*il+ilm 3049 pawtab%indlmn(5,ilmn+ilm)=iln 3050 pawtab%indlmn(6,ilmn+ilm)=1 3051 end do 3052 ilmn=ilmn+2*il+1 3053 end do 3054 LIBPAW_DEALLOCATE(nprj) 3055 !Are ilmn (found here) and pawtab%lmn_size compatibles ? 3056 if (ilmn/=pawtab%lmn_size) then 3057 write(msg, '(a,a,a,a,a)' )& 3058 & 'Calculated lmn size differs from',ch10,& 3059 & 'lmn_size read from pseudo !',ch10,& 3060 & 'Action: check your pseudopotential file.' 3061 MSG_ERROR(msg) 3062 end if 3063 3064 !========================================================== 3065 !Read and initialize radial meshes 3066 3067 nmesh=paw_setuploc%ngrid 3068 LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh)) 3069 LIBPAW_ALLOCATE(mesh_shift,(nmesh)) 3070 do imsh=1,nmesh 3071 radmesh(imsh)%mesh_type=-1 3072 radmesh(imsh)%rstep=zero 3073 radmesh(imsh)%lstep=zero 3074 mesh_shift(imsh)=0 3075 select case(trim(paw_setuploc%radial_grid(imsh)%eq)) 3076 case("r=a*exp(d*i)") 3077 mesh_shift(imsh)=1 3078 radmesh(imsh)%mesh_type=3 3079 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3080 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3081 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa 3082 radmesh(imsh)%lstep=paw_setuploc%radial_grid(imsh)%dd 3083 case("r=a*i/(1-b*i)") 3084 write(msg, '(3a)' )& 3085 & 'The grid r=a*i/(1-b*i) is not implemented in ABINIT !',ch10,& 3086 & 'Action: check your psp file.' 3087 MSG_ERROR(msg) 3088 case("r=a*i/(n-i)") 3089 mesh_shift(imsh)=0 3090 radmesh(imsh)%mesh_type=5 3091 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3092 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3093 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa 3094 radmesh(imsh)%lstep=dble(paw_setuploc%radial_grid(imsh)%nn) 3095 case("r=a*(exp(d*i)-1)") 3096 mesh_shift(imsh)=0 3097 radmesh(imsh)%mesh_type=2 3098 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3099 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3100 if(paw_setuploc%radial_grid(imsh)%istart==1)radmesh(imsh)%mesh_size=radmesh(imsh)%mesh_size+1 3101 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%aa 3102 radmesh(imsh)%lstep=paw_setuploc%radial_grid(imsh)%dd 3103 case("r=d*i") 3104 mesh_shift(imsh)=0 3105 radmesh(imsh)%mesh_type=1 3106 radmesh(imsh)%mesh_size=paw_setuploc%radial_grid(imsh)%iend & 3107 & -paw_setuploc%radial_grid(imsh)%istart+1+mesh_shift(imsh) 3108 if(paw_setuploc%radial_grid(imsh)%istart==1)radmesh(imsh)%mesh_size=radmesh(imsh)%mesh_size+1 3109 radmesh(imsh)%rstep=paw_setuploc%radial_grid(imsh)%dd 3110 case("r=(i/n+a)^5/a-a^4") 3111 write(msg, '(3a)' )& 3112 & 'The grid r=(i/n+a)^5/a-a^4 is not implemented in ABINIT !',ch10,& 3113 & 'Action: check your psp file.' 3114 MSG_ERROR(msg) 3115 end select 3116 end do 3117 3118 !Initialize radial meshes 3119 do imsh=1,nmesh 3120 call pawrad_init(radmesh(imsh)) 3121 end do 3122 3123 pawtab%rpaw=pawpsp_header%rpaw 3124 3125 !========================================================== 3126 !Here reading shapefunction parameters 3127 3128 pawtab%shape_type=pawpsp_header%shape_type 3129 pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99 3130 pawtab%rshp=pawpsp_header%rshp 3131 pawtab%shape_lambda=paw_setuploc%shape_function%lamb 3132 if(trim(paw_setuploc%shape_function%gtype)=="gauss")pawtab%shape_lambda=2 3133 pawtab%shape_sigma=paw_setuploc%shape_function%rc 3134 !If shapefunction type is gaussian, check exponent 3135 if (pawtab%shape_type==1) then 3136 if (pawtab%shape_lambda<2) then 3137 write(msg, '(3a)' )& 3138 & 'For a gaussian shape function, exponent lambda must be >1 !',ch10,& 3139 & 'Action: check your psp file.' 3140 MSG_ERROR(msg) 3141 end if 3142 end if 3143 3144 !If shapefunction type is Bessel, deduce here its parameters from rc 3145 if (pawtab%shape_type==3) then 3146 LIBPAW_ALLOCATE(pawtab%shape_alpha,(2,pawtab%l_size)) 3147 LIBPAW_ALLOCATE(pawtab%shape_q,(2,pawtab%l_size)) 3148 rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw 3149 do il=1,pawtab%l_size 3150 call atompaw_shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc) 3151 end do 3152 end if 3153 3154 !========================================================== 3155 !Mirror pseudopotential parameters to the output and log files 3156 3157 write(msg,'(a,i2)')' Pseudopotential format is: paw',pspversion 3158 call wrtout(ab_out,msg,'COLL') 3159 call wrtout(std_out, msg,'COLL') 3160 write(msg,'(2(a,i3),a,64i4)') & 3161 & ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',& 3162 & pawtab%lmn_size,'), orbitals=',pawtab%orbitals(1:pawtab%basis_size) 3163 call wrtout(ab_out,msg,'COLL') 3164 call wrtout(std_out, msg,'COLL') 3165 write(msg,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw 3166 call wrtout(ab_out,msg,'COLL') 3167 call wrtout(std_out, msg,'COLL') 3168 write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:' 3169 call wrtout(ab_out,msg,'COLL') 3170 call wrtout(std_out, msg,'COLL') 3171 3172 do imsh=1,nmesh 3173 if (radmesh(imsh)%mesh_type==1) & 3174 & write(msg,'(a,i1,a,i4,a,g12.5)') & 3175 & ' - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,& 3176 & ' , step=',radmesh(imsh)%rstep 3177 if (radmesh(imsh)%mesh_type==2) & 3178 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 3179 & ' - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,& 3180 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 3181 if (radmesh(imsh)%mesh_type==3) & 3182 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 3183 & ' - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,& 3184 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 3185 if (radmesh(imsh)%mesh_type==4) & 3186 & write(msg,'(a,i1,a,i4,a,g12.5)') & 3187 & ' - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,& 3188 & ' , AA=',radmesh(imsh)%rstep 3189 if (radmesh(imsh)%mesh_type==5) & 3190 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 3191 & ' - mesh ',imsh,': r(i)=-AA*i/(NN-i)), n=size=',radmesh(imsh)%mesh_size,& 3192 & ' , AA=',radmesh(imsh)%rstep,' NN=',radmesh(imsh)%lstep 3193 call wrtout(ab_out,msg,'COLL') 3194 call wrtout(std_out, msg,'COLL') 3195 end do 3196 if (pawtab%shape_type==-1) then 3197 write(msg,'(a)')& 3198 ' Shapefunction is NUMERIC type: directly read from atomic data file' 3199 call wrtout(ab_out,msg,'COLL') 3200 call wrtout(std_out, msg,'COLL') 3201 end if 3202 if (pawtab%shape_type==1) then 3203 write(msg,'(2a,a,f6.3,a,i3)')& 3204 & ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,& 3205 & ' with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda 3206 call wrtout(ab_out,msg,'COLL') 3207 call wrtout(std_out, msg,'COLL') 3208 end if 3209 if (pawtab%shape_type==2) then 3210 write(msg,'(a)')& 3211 ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2' 3212 call wrtout(ab_out,msg,'COLL') 3213 call wrtout(std_out, msg,'COLL') 3214 end if 3215 if (pawtab%shape_type==3) then 3216 write(msg,'(a)')& 3217 & ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)' 3218 call wrtout(ab_out,msg,'COLL') 3219 call wrtout(std_out, msg,'COLL') 3220 end if 3221 if (pawtab%rshp<1.d-8) then 3222 write(msg,'(a)') ' Radius for shape functions = sphere core radius' 3223 else 3224 write(msg,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp 3225 end if 3226 call wrtout(ab_out,msg,'COLL') 3227 call wrtout(std_out, msg,'COLL') 3228 3229 !========================================================== 3230 !Perfom tests 3231 3232 !Are lmax and orbitals compatibles ? 3233 if (lmax/=maxval(pawtab%orbitals)) then 3234 write(msg, '(a,a,a)' )& 3235 & 'lmax /= MAX(orbitals) !',ch10,& 3236 & 'Action: check your pseudopotential file.' 3237 MSG_ERROR(msg) 3238 end if 3239 3240 !Only mesh_type=1,2, 3 or 5 allowed 3241 do imsh=1,nmesh 3242 if (radmesh(imsh)%mesh_type>5) then 3243 write(msg, '(a,a,a)' )& 3244 & 'Only mesh types 1,2,3 or 5 allowed !',ch10,& 3245 & 'Action : check your pseudopotential or input file.' 3246 MSG_ERROR(msg) 3247 end if 3248 end do 3249 3250 !========================================================== 3251 !Read tabulated atomic data 3252 3253 !--------------------------------- 3254 !Read wave-functions (phi) 3255 3256 do ib=1,pawtab%basis_size 3257 3258 if (ib==1) then 3259 do imsh=1,nmesh 3260 if(trim(paw_setuploc%ae_partial_wave(1)%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3261 mmax=radmesh(imsh)%mesh_size 3262 call pawrad_init(pawrad,mesh_size=mmax,mesh_type=radmesh(imsh)%mesh_type, & 3263 & rstep=radmesh(imsh)%rstep,lstep=radmesh(imsh)%lstep,r_for_intg=pawtab%rpaw) 3264 pawtab%partialwave_mesh_size=pawrad%mesh_size 3265 pawtab%mesh_size=pawrad_ifromr(pawrad,pawtab%rpaw)+5 3266 pawtab%mesh_size=min(pawtab%mesh_size,pawrad%mesh_size) 3267 imainmesh=imsh 3268 exit 3269 end if 3270 end do 3271 LIBPAW_ALLOCATE(pawtab%phi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 3272 else if (trim(paw_setuploc%ae_partial_wave(ib)%grid)/=trim(paw_setuploc%radial_grid(imainmesh)%id)) then 3273 write(msg, '(a,a,a)' )& 3274 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 3275 & 'Action: check your pseudopotential file.' 3276 MSG_ERROR(msg) 3277 end if 3278 shft=mesh_shift(imainmesh) 3279 pawtab%phi(1+shft:pawtab%partialwave_mesh_size,ib)= & 3280 & paw_setuploc%ae_partial_wave(ib)%data(1:pawtab%partialwave_mesh_size-shft) & 3281 & *pawrad%rad(1+shft:pawtab%partialwave_mesh_size) 3282 if (shft==1) pawtab%phi(1,ib)=zero 3283 end do 3284 write(msg,'(a,i4)') & 3285 & ' mmax= ',mmax 3286 call wrtout(ab_out,msg,'COLL') 3287 call wrtout(std_out, msg,'COLL') 3288 3289 !--------------------------------- 3290 !Read pseudo wave-functions (tphi) 3291 3292 LIBPAW_ALLOCATE(pawtab%tphi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 3293 do ib=1,pawtab%basis_size 3294 3295 if(trim(paw_setuploc%pseudo_partial_wave(ib)%grid)/=trim(paw_setuploc%radial_grid(imainmesh)%id)) then 3296 write(msg, '(a,a,a)' )& 3297 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 3298 & 'Action: check your pseudopotential file.' 3299 MSG_ERROR(msg) 3300 end if 3301 shft=mesh_shift(imainmesh) 3302 pawtab%tphi(1+shft:pawtab%partialwave_mesh_size,ib)=& 3303 & paw_setuploc%pseudo_partial_wave(ib)%data(1:pawtab%partialwave_mesh_size-shft) & 3304 & *pawrad%rad(1+shft:pawtab%partialwave_mesh_size) 3305 if (shft==1) pawtab%tphi(1,ib)=zero 3306 end do 3307 write(msg,'(a,i1)') & 3308 & ' Radial grid used for partial waves is grid ',imainmesh 3309 call wrtout(ab_out,msg,'COLL') 3310 call wrtout(std_out, msg,'COLL') 3311 3312 !--------------------------------- 3313 !Read projectors (tproj) 3314 3315 if (allocated(paw_setuploc%projector_fit)) then 3316 call wvlpaw_allocate(pawtab%wvl) 3317 LIBPAW_ALLOCATE(pawtab%wvl%pngau,(pawtab%basis_size)) 3318 do ib=1,pawtab%basis_size 3319 pawtab%wvl%pngau(ib) = paw_setuploc%projector_fit(ib)%ngauss 3320 end do 3321 pawtab%wvl%ptotgau = sum(pawtab%wvl%pngau) * 2 3322 LIBPAW_ALLOCATE(pawtab%wvl%parg,(2,pawtab%wvl%ptotgau)) 3323 LIBPAW_ALLOCATE(pawtab%wvl%pfac,(2,pawtab%wvl%ptotgau)) 3324 pngau = 1 3325 do ib=1,pawtab%basis_size 3326 ! Complex gaussian 3327 pawtab%wvl%parg(:,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3328 & paw_setuploc%projector_fit(ib)%expos(:,1:pawtab%wvl%pngau(ib)) 3329 pawtab%wvl%pfac(:,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3330 & paw_setuploc%projector_fit(ib)%factors(:,1:pawtab%wvl%pngau(ib)) 3331 pngau = pngau + pawtab%wvl%pngau(ib) 3332 ! Conjugate gaussian 3333 pawtab%wvl%parg(1,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3334 & paw_setuploc%projector_fit(ib)%expos(1,1:pawtab%wvl%pngau(ib)) 3335 pawtab%wvl%parg(2,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3336 & -paw_setuploc%projector_fit(ib)%expos(2,1:pawtab%wvl%pngau(ib)) 3337 pawtab%wvl%pfac(1,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3338 & paw_setuploc%projector_fit(ib)%factors(1,1:pawtab%wvl%pngau(ib)) 3339 pawtab%wvl%pfac(2,pngau:pngau + pawtab%wvl%pngau(ib) - 1) = & 3340 & -paw_setuploc%projector_fit(ib)%factors(2,1:pawtab%wvl%pngau(ib)) 3341 pngau = pngau + pawtab%wvl%pngau(ib) 3342 pawtab%wvl%pngau(ib) = pawtab%wvl%pngau(ib) * 2 3343 end do 3344 pawtab%has_wvl=2 3345 else 3346 !Nullify wavelet objects for safety: 3347 pawtab%has_wvl=0 3348 call wvlpaw_free(pawtab%wvl) 3349 end if 3350 3351 do ib=1,pawtab%basis_size 3352 if (ib==1) then 3353 do imsh=1,nmesh 3354 if(trim(paw_setuploc%projector_function(1)%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3355 iprojmesh=imsh 3356 exit 3357 end if 3358 end do 3359 call pawrad_copy(radmesh(iprojmesh),tproj_mesh) 3360 LIBPAW_ALLOCATE(tproj,(tproj_mesh%mesh_size,pawtab%basis_size)) 3361 else if (trim(paw_setuploc%projector_function(ib)%grid)/=trim(paw_setuploc%radial_grid(iprojmesh)%id)) then 3362 write(msg, '(a,a,a)' )& 3363 & 'All tprojectors must be given on the same radial mesh !',ch10,& 3364 & 'Action: check your pseudopotential file.' 3365 MSG_ERROR(msg) 3366 end if 3367 shft=mesh_shift(iprojmesh) 3368 tproj(1+shft:tproj_mesh%mesh_size,ib)=paw_setuploc%projector_function(ib)%data(1:tproj_mesh%mesh_size-shft)& 3369 & *tproj_mesh%rad(1+shft:tproj_mesh%mesh_size) 3370 if (shft==1) tproj(1,ib)=zero 3371 end do 3372 write(msg,'(a,i1)') & 3373 & ' Radial grid used for projectors is grid ',iprojmesh 3374 call wrtout(ab_out,msg,'COLL') 3375 call wrtout(std_out, msg,'COLL') 3376 3377 3378 !--------------------------------- 3379 !Read core density (coredens) 3380 do imsh=1,nmesh 3381 if(trim(paw_setuploc%ae_core_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3382 icoremesh=imsh 3383 exit 3384 end if 3385 end do 3386 call pawrad_copy(radmesh(icoremesh),core_mesh) 3387 if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.& 3388 & (radmesh(icoremesh)%rstep /=pawrad%rstep) .or.& 3389 & (radmesh(icoremesh)%lstep /=pawrad%lstep)) then 3390 write(msg, '(a,a,a,a,a)' )& 3391 & 'Ncore must be given on a radial mesh with the same',ch10,& 3392 & 'type and step(s) than the main radial mesh (mesh for Phi) !',ch10,& 3393 & 'Action: check your pseudopotential file.' 3394 MSG_ERROR(msg) 3395 end if 3396 LIBPAW_ALLOCATE(ncore,(core_mesh%mesh_size)) 3397 shft=mesh_shift(icoremesh) 3398 ncore(1+shft:core_mesh%mesh_size)=paw_setuploc%ae_core_density%data(1:core_mesh%mesh_size-shft)/sqrt(fourpi) 3399 if (shft==1) call pawrad_deducer0(ncore,core_mesh%mesh_size,core_mesh) 3400 3401 !Construct and save VH[z_NC] if requested 3402 if (pawtab%has_vhnzc==1) then 3403 LIBPAW_ALLOCATE(pawtab%VHnZC,(pawtab%mesh_size)) 3404 LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size)) 3405 call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl) 3406 pawtab%VHnZC(1:pawtab%mesh_size)=vhnzc(1:pawtab%mesh_size) 3407 pawtab%has_vhnzc=2 3408 LIBPAW_DEALLOCATE(vhnzc) 3409 end if 3410 3411 pawtab%core_mesh_size=pawtab%mesh_size 3412 if(save_core_msz) pawtab%core_mesh_size=core_mesh%mesh_size 3413 LIBPAW_ALLOCATE(pawtab%coredens,(pawtab%core_mesh_size)) 3414 pawtab%rcore=core_mesh%rad(pawtab%core_mesh_size) 3415 pawtab%coredens(1:pawtab%core_mesh_size)=ncore(1:pawtab%core_mesh_size) 3416 3417 !--------------------------------- 3418 !Read pseudo core density (tcoredens) 3419 do imsh=1,nmesh 3420 if(trim(paw_setuploc%pseudo_core_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3421 iread1=imsh 3422 exit 3423 end if 3424 end do 3425 if (iread1/=icoremesh) then 3426 write(msg, '(a,a,a,a,a,a,a,a)' )& 3427 & 'Pseudized core density (tNcore) must be given',ch10,& 3428 & 'on the same radial mesh as core density (Ncore) !',ch10,& 3429 & 'Action: check your pseudopotential file.' 3430 MSG_ERROR(msg) 3431 end if 3432 LIBPAW_ALLOCATE(tncore,(core_mesh%mesh_size)) 3433 shft=mesh_shift(icoremesh) 3434 tncore(1+shft:core_mesh%mesh_size)=paw_setuploc%pseudo_core_density%data(1:core_mesh%mesh_size-shft)/sqrt(fourpi) 3435 if (shft==1) call pawrad_deducer0(tncore,core_mesh%mesh_size,core_mesh) 3436 if(save_core_msz) then 3437 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,6)) 3438 else 3439 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,1)) 3440 end if 3441 if (maxval(abs(tncore(:)))<tol6) then 3442 pawtab%usetcore=0 3443 pawtab%tcoredens(1:pawtab%core_mesh_size,:)=zero 3444 else 3445 pawtab%usetcore=1 3446 pawtab%tcoredens(1:pawtab%core_mesh_size,1)=tncore(1:pawtab%core_mesh_size) 3447 end if 3448 write(msg,'(a,i1)') & 3449 & ' Radial grid used for (t)core density is grid ',icoremesh 3450 call wrtout(ab_out,msg,'COLL') 3451 call wrtout(std_out, msg,'COLL') 3452 3453 3454 !--------------------------------- 3455 !Read local pseudopotential=Vh(tn_zc) or Vbare 3456 3457 if ((paw_setuploc%blochl_local_ionic_potential%tread).and.& 3458 & (pawtab%usexcnhat==-1.or.pawtab%usexcnhat==0.or.(pawtab%usexcnhat==1.and.& 3459 & ((.not.paw_setuploc%zero_potential%tread).or.(.not.paw_setuploc%kresse_joubert_local_ionic_potential%tread))))) then 3460 usexcnhat=0;vlocopt=2 3461 do imsh=1,nmesh 3462 if(trim(paw_setuploc%blochl_local_ionic_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3463 iread1=imsh 3464 exit 3465 end if 3466 end do 3467 ivlocmesh=iread1 3468 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 3469 LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 3470 shft=mesh_shift(ivlocmesh) 3471 vlocr(1+shft:vloc_mesh%mesh_size)=paw_setuploc%blochl_local_ionic_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3472 if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh) 3473 else if((paw_setuploc%kresse_joubert_local_ionic_potential%tread).and.& 3474 & (pawtab%usexcnhat==-1.or.pawtab%usexcnhat==1.or.(pawtab%usexcnhat==0.and.& 3475 & (.not.paw_setuploc%zero_potential%tread)))) then 3476 usexcnhat=1;vlocopt=1 3477 do imsh=1,nmesh 3478 if(trim(paw_setuploc%kresse_joubert_local_ionic_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3479 iread1=imsh 3480 exit 3481 end if 3482 end do 3483 ivlocmesh=iread1 3484 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 3485 LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 3486 shft=mesh_shift(ivlocmesh) 3487 vlocr(1+shft:vloc_mesh%mesh_size)= & 3488 & paw_setuploc%kresse_joubert_local_ionic_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3489 if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh) 3490 else if(paw_setuploc%zero_potential%tread) then 3491 usexcnhat=0;vlocopt=0 3492 do imsh=1,nmesh 3493 if(trim(paw_setuploc%zero_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3494 iread1=imsh 3495 exit 3496 end if 3497 end do 3498 ivlocmesh=iread1 3499 vloc_mesh%mesh_type=radmesh(ivlocmesh)%mesh_type 3500 vloc_mesh%rstep=radmesh(ivlocmesh)%rstep 3501 vloc_mesh%lstep=radmesh(ivlocmesh)%lstep 3502 vloc_mesh%mesh_size=pawrad_ifromr(radmesh(ivlocmesh),rmax_vloc) 3503 call pawrad_init(vloc_mesh) 3504 ! call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 3505 LIBPAW_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 3506 vlocr=zero 3507 shft=mesh_shift(ivlocmesh) 3508 vlocr(1+shft:vloc_mesh%mesh_size)=paw_setuploc%zero_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3509 if (shft==1) call pawrad_deducer0(vlocr,vloc_mesh%mesh_size,vloc_mesh) 3510 else 3511 write(msg, '(a,a,a,a,a)' )& 3512 & 'At least one local potential must be given',ch10,& 3513 & 'Action: check your pseudopotential file.' 3514 MSG_ERROR(msg) 3515 end if 3516 3517 write(msg,'(a,i1)') & 3518 & ' Radial grid used for Vloc is grid ',ivlocmesh 3519 call wrtout(ab_out,msg,'COLL') 3520 call wrtout(std_out, msg,'COLL') 3521 3522 !------------------------------------------------- 3523 !Read LDA-1/2 potential 3524 if (paw_setuploc%LDA_minus_half_potential%tread) then 3525 do imsh=1,nmesh 3526 if(trim(paw_setuploc%LDA_minus_half_potential%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3527 iread1=imsh 3528 exit 3529 end if 3530 end do 3531 if(iread1/=ivlocmesh) then 3532 write(msg, '(a)' )& 3533 & 'The LDA-1/2 potential must be given on the same grid as the local potential.' 3534 MSG_ERROR(msg) 3535 end if 3536 has_v_minushalf=1 3537 LIBPAW_ALLOCATE(pawtab%vminushalf,(vloc_mesh%mesh_size)) 3538 shft=mesh_shift(ivlocmesh) 3539 pawtab%vminus_mesh_size=vloc_mesh%mesh_size 3540 pawtab%vminushalf(1+shft:vloc_mesh%mesh_size)= & 3541 & paw_setuploc%LDA_minus_half_potential%data(1:vloc_mesh%mesh_size-shft)/sqrt(fourpi) 3542 if (shft==1) call pawrad_deducer0(pawtab%vminushalf,vloc_mesh%mesh_size,vloc_mesh) 3543 write(msg,'(a,i1)') & 3544 & ' Radial grid used for LDA-1/2 potential is grid ',ivlocmesh 3545 call wrtout(ab_out,msg,'COLL') 3546 call wrtout(std_out, msg,'COLL') 3547 else 3548 LIBPAW_ALLOCATE(pawtab%vminushalf,(0)) 3549 has_v_minushalf=0 3550 end if 3551 if(has_v_minushalf==0.and.pawtab%has_vminushalf==1) then 3552 write(msg, '(a)' )& 3553 & 'The LDA-1/2 potential must be given in the XML PAW datafile.' 3554 MSG_ERROR(msg) 3555 end if 3556 3557 !--------------------------------- 3558 !Eventually read "numeric" shapefunctions (if shape_type=-1) 3559 if (pawtab%shape_type==-1) then 3560 LIBPAW_ALLOCATE(pawtab%shapefunc,(pawtab%mesh_size,pawtab%l_size)) 3561 do imsh=1,nmesh 3562 if(trim(paw_setuploc%shape_function%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3563 iread1=imsh 3564 exit 3565 end if 3566 end do 3567 call pawrad_copy(radmesh(iread1),shpf_mesh) 3568 ishpfmesh=iread1 3569 LIBPAW_ALLOCATE(shpf,(shpf_mesh%mesh_size,pawtab%l_size)) 3570 shft=mesh_shift(ishpfmesh) 3571 shpf(1,1)=one 3572 do ir=2,shpf_mesh%mesh_size 3573 shpf(ir,1)=paw_setuploc%shape_function%data(ir-shft,1) 3574 end do 3575 sz10=size(paw_setuploc%shape_function%data,2) 3576 if(sz10>=2) then 3577 do il=2,pawtab%l_size 3578 shpf(1,il)=zero 3579 do ir=2,shpf_mesh%mesh_size 3580 shpf(ir,il)=paw_setuploc%shape_function%data(ir-shft,il) 3581 end do 3582 end do 3583 else 3584 do il=2,pawtab%l_size 3585 shpf(1,il)=zero 3586 do ir=2,shpf_mesh%mesh_size 3587 shpf(ir,il)=paw_setuploc%shape_function%data(ir-shft,1)*shpf_mesh%rad(ir)**(il-1) 3588 end do 3589 end do 3590 end if 3591 write(msg,'(a,i1)') & 3592 & ' Radial grid used for shape functions is grid ',iread1 3593 call wrtout(ab_out,msg,'COLL') 3594 call wrtout(std_out, msg,'COLL') 3595 3596 ! Has to spline shape functions if mesh is not the "main" mesh 3597 if (ishpfmesh/=imainmesh) then 3598 msz=shpf_mesh%mesh_size 3599 LIBPAW_ALLOCATE(work1,(msz)) 3600 LIBPAW_ALLOCATE(work2,(msz)) 3601 LIBPAW_ALLOCATE(work3,(msz)) 3602 LIBPAW_ALLOCATE(work4,(pawrad%mesh_size)) 3603 work3(1:msz)=shpf_mesh%rad(1:msz) 3604 work4(1:pawrad%mesh_size)=pawrad%rad(1:pawrad%mesh_size) 3605 do il=1,pawtab%l_size 3606 call bound_deriv(shpf(1:msz,il),shpf_mesh,msz,yp1,ypn) 3607 call paw_spline(work3,shpf(:,il),msz,yp1,ypn,work1) 3608 call paw_splint(msz,work3,shpf(:,il),work1,pawrad%mesh_size,work4,pawtab%shapefunc(:,il)) 3609 end do 3610 LIBPAW_DEALLOCATE(work1) 3611 LIBPAW_DEALLOCATE(work2) 3612 LIBPAW_DEALLOCATE(work3) 3613 LIBPAW_DEALLOCATE(work4) 3614 else 3615 pawtab%shapefunc(:,:)=shpf(:,:) 3616 end if 3617 LIBPAW_DEALLOCATE(shpf) 3618 end if 3619 3620 !--------------------------------- 3621 !Read pseudo valence density 3622 if (paw_setuploc%pseudo_valence_density%tread) then 3623 do imsh=1,nmesh 3624 if(trim(paw_setuploc%pseudo_valence_density%grid)==trim(paw_setuploc%radial_grid(imsh)%id)) then 3625 iread1=imsh 3626 exit 3627 end if 3628 end do 3629 ivalemesh=iread1 3630 call pawrad_copy(radmesh(iread1),vale_mesh) 3631 LIBPAW_ALLOCATE(tnvale,(vale_mesh%mesh_size)) 3632 shft=mesh_shift(ivalemesh) 3633 tnvale(1+shft:vale_mesh%mesh_size)=paw_setuploc%pseudo_valence_density%data(1:vale_mesh%mesh_size-shft)/sqrt(fourpi) 3634 if (shft==1) call pawrad_deducer0(tnvale,vale_mesh%mesh_size,vale_mesh) 3635 pawtab%has_tvale=1 3636 write(msg,'(a,i1)') & 3637 & ' Radial grid used for pseudo valence density is grid ',ivalemesh 3638 call wrtout(ab_out,msg,'COLL') 3639 call wrtout(std_out, msg,'COLL') 3640 else 3641 pawtab%has_tvale=0 3642 LIBPAW_ALLOCATE(tnvale,(0)) 3643 end if 3644 3645 !--------------------------------- 3646 !Read initial guess of rhoij (rhoij0) 3647 3648 LIBPAW_ALLOCATE(pawtab%rhoij0,(pawtab%lmn2_size)) 3649 pawtab%rhoij0=zero 3650 ilmn0=0 3651 do ib=1,pawtab%basis_size 3652 il=2*pawtab%orbitals(ib)+1 3653 occ=paw_setuploc%valence_states%state(ib)%ff 3654 if (occ<zero)occ=zero 3655 do ilmn=ilmn0+1,ilmn0+il 3656 pawtab%rhoij0(ilmn*(ilmn+1)/2)=occ/dble(il) 3657 end do 3658 ilmn0=ilmn0+il 3659 end do 3660 3661 !--------------------------------- 3662 !Read Kij terms (kij0) and deduce eventually Dij0 3663 3664 LIBPAW_ALLOCATE(kij,(pawtab%lmn2_size)) 3665 kij=zero 3666 nval=paw_setuploc%valence_states%nval 3667 do jlmn=1,pawtab%lmn_size 3668 j0lmn=jlmn*(jlmn-1)/2 3669 jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn) 3670 do ilmn=1,jlmn 3671 klmn=j0lmn+ilmn 3672 ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn) 3673 if (ilm==jlm) kij(klmn)=paw_setuploc%kinetic_energy_differences%data(jln+(iln-1)*nval) 3674 end do 3675 end do 3676 if (vlocopt>0) then 3677 LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size)) 3678 if (size(pawtab%vminushalf)>0.and.pawtab%has_vminushalf==1) then 3679 vlocr(1:vloc_mesh%mesh_size)=vlocr(1:vloc_mesh%mesh_size)+pawtab%vminushalf(1:vloc_mesh%mesh_size) 3680 end if 3681 call atompaw_dij0(pawtab%indlmn,kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,& 3682 & vloc_mesh,vlocr,znucl) 3683 end if 3684 3685 !Keep eventualy Kij in memory 3686 if (pawtab%has_kij==1.or.vlocopt==0) then 3687 LIBPAW_ALLOCATE(pawtab%kij,(pawtab%lmn2_size)) 3688 pawtab%kij(:)=kij(:) 3689 if (vlocopt> 0) pawtab%has_kij=2 3690 ! This -1 means that pawtab%kij will be freed later 3691 if (vlocopt==0) pawtab%has_kij=-1 3692 end if 3693 3694 LIBPAW_DEALLOCATE(kij) 3695 3696 !--------------------------------- 3697 !Read exact-exchange Fock terms for core-valence interactions (ex_cvij) 3698 if (paw_setuploc%exact_exchange_matrix%tread.eqv..true.) then 3699 pawtab%has_fock=2 3700 LIBPAW_ALLOCATE(pawtab%ex_cvij,(pawtab%lmn2_size)) 3701 pawtab%ex_cvij=zero 3702 nval=paw_setuploc%valence_states%nval 3703 do jlmn=1,pawtab%lmn_size 3704 j0lmn=jlmn*(jlmn-1)/2 3705 jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn) 3706 do ilmn=1,jlmn 3707 klmn=j0lmn+ilmn 3708 ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn) 3709 if (ilm==jlm) pawtab%ex_cvij(klmn)=paw_setuploc%exact_exchange_matrix%data(jln+(iln-1)*nval) 3710 end do 3711 end do 3712 pawtab%ex_cc=paw_setuploc%ex_cc 3713 end if 3714 3715 !========================================================== 3716 !Compute additional atomic data only depending on present DATASET 3717 3718 call pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,ixc,lnmax,& 3719 & mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,& 3720 & qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,& 3721 & vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl) 3722 3723 if(usewvl==1 .or. icoulomb > 0) then 3724 ! Calculate up to the 5th derivative of tcoredens 3725 call pawpsp_calc_d5(core_mesh,pawtab%core_mesh_size,pawtab%tcoredens) 3726 ! Other wvl related operations 3727 call pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr) 3728 else if (pawtab%has_wvl>0) then 3729 call wvlpaw_rholoc_nullify(pawtab%wvl%rholoc) 3730 end if 3731 3732 !========================================================== 3733 !Free temporary allocated space 3734 3735 call pawrad_free(radmesh) 3736 LIBPAW_DATATYPE_DEALLOCATE(radmesh) 3737 LIBPAW_DEALLOCATE(mesh_shift) 3738 3739 call pawrad_free(tproj_mesh) 3740 call pawrad_free(core_mesh) 3741 call pawrad_free(vloc_mesh) 3742 3743 LIBPAW_DEALLOCATE(vlocr) 3744 LIBPAW_DEALLOCATE(ncore) 3745 LIBPAW_DEALLOCATE(tncore) 3746 LIBPAW_DEALLOCATE(tproj) 3747 3748 if(pawtab%shape_type==-1) then 3749 call pawrad_free(shpf_mesh) 3750 end if 3751 if (paw_setuploc%pseudo_valence_density%tread) then 3752 call pawrad_free(vale_mesh) 3753 end if 3754 3755 LIBPAW_DEALLOCATE(tnvale) 3756 3757 end subroutine pawpsp_17in
m_pawpsp/pawpsp_7in [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_7in
FUNCTION
Initialize pspcod=7 ("PAW pseudopotentials"):
continue to read the corresponding file and compute the form factors
INPUTS
icoulomb==0 : usual reciprocal space computation =1 : free boundary conditions are used ipsp=id in the array of the currently read pseudo. ixc=exchange-correlation choice from main routine data file lloc=angular momentum choice of local pseudopotential lmax=value of lmax mentioned at the second line of the psp file pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) xclevel= XC functional level zion=nominal valence of atom as specified in psp file
OUTPUT
epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
ffspl(psps%mqgrid_ff,2,psps%lnmax)=form factor f_l(q) and second derivative
from spline fit for each angular momentum and each projector;
pawrad <type(pawrad_type)>=paw radial mesh and related data
pawtab <type(pawtab_type)>=paw tabulated starting data
vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
xcccrc=XC core correction cutoff radius (bohr) from psp file
NOTES
Spin-orbit not yet implemented (to be done)
PARENTS
m_pawpsp,pspatm
CHILDREN
SOURCE
3801 subroutine pawpsp_7in(epsatm,ffspl,icoulomb,ixc,& 3802 & lmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,& 3803 & pawrad,pawtab,pawxcdev,qgrid_ff,qgrid_vl,& 3804 & usewvl,usexcnhat_in,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl) 3805 3806 3807 !This section has been created automatically by the script Abilint (TD). 3808 !Do not modify the following lines by hand. 3809 #undef ABI_FUNC 3810 #define ABI_FUNC 'pawpsp_7in' 3811 !End of the abilint section 3812 3813 implicit none 3814 3815 !Arguments ------------------------------------ 3816 !scalars 3817 integer, intent(in):: icoulomb,ixc 3818 integer, intent(in):: lmax,lnmax,mmax 3819 integer, intent(in):: mqgrid_ff,mqgrid_vl,pawxcdev 3820 integer, intent(in):: usewvl,usexcnhat_in,xclevel 3821 real(dp), intent(in):: xc_denpos,zion,znucl 3822 real(dp), intent(out):: epsatm,xcccrc 3823 type(pawrad_type), intent(inout):: pawrad 3824 type(pawtab_type), intent(inout) :: pawtab 3825 !arrays 3826 real(dp),intent(in):: qgrid_ff(mqgrid_ff),qgrid_vl(mqgrid_vl) 3827 real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax) 3828 real(dp),intent(out) :: vlspl(mqgrid_vl,2) 3829 3830 !Local variables ------------------------------ 3831 !scalars 3832 integer :: imainmesh 3833 integer :: nmesh 3834 integer :: pspversion,usexcnhat,vlocopt 3835 logical :: save_core_msz 3836 type(pawrad_type) :: core_mesh,tproj_mesh,vale_mesh,vloc_mesh 3837 real(dp),pointer :: ncore(:),tncore(:),tnvale(:),tproj(:,:),vlocr(:) 3838 type(pawrad_type),pointer :: radmesh(:) 3839 3840 !************************************************************************ 3841 3842 !Destroy everything in pawtab but optional flags 3843 call pawtab_free(pawtab) 3844 !Destroy everything in pawrad 3845 call pawrad_free(pawrad) 3846 3847 save_core_msz=(usewvl==1 .or. icoulomb .ne. 0) 3848 nullify(ncore);nullify(tncore);nullify(tnvale) 3849 nullify(tproj);nullify(vlocr) 3850 nullify(radmesh) 3851 3852 call pawpsp_read(core_mesh,tmp_unit,imainmesh,lmax,& 3853 & ncore,nmesh,pawrad,pawtab,pspversion,radmesh,save_core_msz,& 3854 & tncore,tnvale,tproj,tproj_mesh,usexcnhat_in,usexcnhat,vale_mesh,vlocopt,vlocr,& 3855 & vloc_mesh,znucl) 3856 3857 call pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,ixc,lnmax,& 3858 & mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,& 3859 & qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,& 3860 & vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl) 3861 3862 if(usewvl==1 .or. icoulomb > 0) then 3863 ! Calculate up to the 5th derivative of tcoredens 3864 call pawpsp_calc_d5(core_mesh,pawtab%core_mesh_size,pawtab%tcoredens) 3865 ! Other wvl related operations 3866 call pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr) 3867 else if (pawtab%has_wvl>0) then 3868 call wvlpaw_rholoc_nullify(pawtab%wvl%rholoc) 3869 end if 3870 3871 !========================================================== 3872 !Free temporary allocated space 3873 call pawrad_free(radmesh) 3874 call pawrad_free(tproj_mesh) 3875 call pawrad_free(core_mesh) 3876 call pawrad_free(vloc_mesh) 3877 LIBPAW_DATATYPE_DEALLOCATE(radmesh) 3878 if (associated(vlocr)) then 3879 LIBPAW_POINTER_DEALLOCATE(vlocr) 3880 end if 3881 if (associated(ncore)) then 3882 LIBPAW_POINTER_DEALLOCATE(ncore) 3883 end if 3884 if (associated(tncore)) then 3885 LIBPAW_POINTER_DEALLOCATE(tncore) 3886 end if 3887 if (associated(tnvale)) then 3888 LIBPAW_POINTER_DEALLOCATE(tnvale) 3889 end if 3890 if (associated(tproj)) then 3891 LIBPAW_POINTER_DEALLOCATE(tproj) 3892 end if 3893 if (pspversion>=4) then 3894 call pawrad_free(vale_mesh) 3895 end if 3896 3897 end subroutine pawpsp_7in
m_pawpsp/pawpsp_bcast [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_bcast
FUNCTION
Communicate paw data to all processors
INPUTS
comm_mpi= communicator used to broadcast data lnmax= Max. number of (l,n) components over all type of psps mqgrid_ff= dimension of ffspl mqgrid_vl= dimension of vlspl
OUTPUT
epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and derivative
pawrad=<type pawrad_type>
pawtab=<type pawtab_type>
vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
xcccrc=XC core correction cutoff radius (bohr) from psp file
PARENTS
m_pawpsp,pspatm
CHILDREN
SOURCE
4656 subroutine pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc) 4657 4658 4659 !This section has been created automatically by the script Abilint (TD). 4660 !Do not modify the following lines by hand. 4661 #undef ABI_FUNC 4662 #define ABI_FUNC 'pawpsp_bcast' 4663 !End of the abilint section 4664 4665 implicit none 4666 4667 !Arguments ------------------------------------ 4668 integer,intent(in) :: comm_mpi 4669 real(dp),intent(inout) :: epsatm,xcccrc 4670 real(dp),intent(inout) :: ffspl(:,:,:),vlspl(:,:) 4671 type(pawrad_type),intent(inout) :: pawrad 4672 type(pawtab_type),intent(inout) :: pawtab 4673 4674 !Local variables------------------------------- 4675 integer :: ierr,ii,me,nn_dpr 4676 integer :: siz_ffspl,siz1_ffspl,siz2_ffspl,siz3_ffspl,siz_vlspl,siz1_vlspl,siz2_vlspl 4677 integer,allocatable :: list_int(:) 4678 real(dp),allocatable :: list_dpr(:) 4679 4680 !************************************************************************* 4681 4682 me=xmpi_comm_rank(comm_mpi) 4683 4684 !Broadcast pawrad 4685 call pawrad_bcast(pawrad,comm_mpi) 4686 4687 !Broadcast pawtab (only data read from file) 4688 call pawtab_bcast(pawtab,comm_mpi,only_from_file=.true.) 4689 4690 !Broadcast the sizes of the arrays 4691 LIBPAW_ALLOCATE(list_int,(5)) 4692 if (me==0) then 4693 siz1_vlspl=size(vlspl,1); list_int(1)=siz1_vlspl 4694 siz2_vlspl=size(vlspl,2); list_int(2)=siz2_vlspl 4695 siz1_ffspl=size(ffspl,1); list_int(3)=siz1_ffspl 4696 siz2_ffspl=size(ffspl,2); list_int(4)=siz2_ffspl 4697 siz3_ffspl=size(ffspl,3); list_int(5)=siz3_ffspl 4698 end if 4699 call xmpi_bcast(list_int,0,comm_mpi,ierr) 4700 if (me/=0) then 4701 siz1_vlspl=list_int(1) 4702 siz2_vlspl=list_int(2) 4703 siz1_ffspl=list_int(3) 4704 siz2_ffspl=list_int(4) 4705 siz3_ffspl=list_int(5) 4706 end if 4707 siz_vlspl=siz1_vlspl*siz2_vlspl 4708 siz_ffspl=siz1_ffspl*siz2_ffspl*siz3_ffspl 4709 LIBPAW_DEALLOCATE(list_int) 4710 4711 !Broadcast the reals 4712 nn_dpr=2+siz_vlspl+siz_ffspl 4713 LIBPAW_ALLOCATE(list_dpr,(nn_dpr)) 4714 if (me==0) then 4715 ii=1 4716 list_dpr(ii)=epsatm ;ii=ii+1 4717 list_dpr(ii)=xcccrc ;ii=ii+1 4718 list_dpr(ii:ii+siz_vlspl-1)=reshape(vlspl,(/siz_vlspl/)) ;ii=ii+siz_vlspl 4719 list_dpr(ii:ii+siz_ffspl-1)=reshape(ffspl,(/siz_ffspl/)) ;ii=ii+siz_ffspl 4720 end if 4721 call xmpi_bcast(list_dpr,0,comm_mpi,ierr) 4722 if (me/=0) then 4723 ii=1 4724 epsatm=list_dpr(ii) ;ii=ii+1 4725 xcccrc=list_dpr(ii) ;ii=ii+1 4726 vlspl=reshape(list_dpr(ii:ii+siz_vlspl-1),(/siz1_vlspl,siz2_vlspl/)) 4727 ii=ii+siz_vlspl 4728 ffspl=reshape(list_dpr(ii:ii+siz_ffspl-1),(/siz1_ffspl,siz2_ffspl,siz3_ffspl/)) 4729 ii=ii+siz_ffspl 4730 end if 4731 LIBPAW_DEALLOCATE(list_dpr) 4732 4733 end subroutine pawpsp_bcast
m_pawpsp/pawpsp_calc [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_calc
FUNCTION
Performs tests and compute data related to pspcod=7 or 17 ("PAW pseudopotentials")
INPUTS
core_mesh<type(pawrad_type)>= radial mesh for the core density imainmesh= serial number of the main mesh ixc=exchange-correlation choice from main routine data file lnmax=max. number of (l,n) components over all type of psps angular momentum of nonlocal pseudopotential mqgrid_ff=dimension of q (or G) grid for nl form factors (array ffspl) mqgrid_vl=dimension of q (or G) grid for Vloc (array vlspl) ncore(core_mesh%mesh_size)= core density nmesh= number of radial meshes pawrad <type(pawrad_type)>=paw radial mesh and related data pawxcdev=choice of XC development (0=no dev. (use of angular mesh) ; 1 or 2=dev. on moments) pspversion= version of the atompaw code used to generate paw data. qgrid_ff(mqgrid_ff)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors qgrid_vl(mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for Vloc radmesh(nmesh)<type(pawrad_type)>=paw radial meshes and related data tncore(core_mesh%mesh_size)= pseudo core density tproj(tproj_mesh%mesh_size)= non-local projectors in real space tproj_mesh<type(pawrad_type)>= radial mesh for the projectors usexcnhat=0 if compensation charge density is not included in XC terms 1 if compensation charge density is included in XC terms vale_mesh<type(pawrad_type)>= radial mesh for the valence density xc_denpos= lowest allowed density (usually for the computation of the XC functionals) vlocopt= option for the local potential.(0=Vbare, 1=VH(tnzc) with hat in XC, 2=VH(tnzc) w/o hat in XC) vlocr(vloc_mesh%mesh_size)= local potential according to vlocopt. xclevel= XC functional level zion=nominal valence of atom as specified in psp file znucl=atomic number of atom as specified in input file to main routine
OUTPUT
epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree)
ffspl(mqgrid_ff,2,lnmax)=form factor f_l(q) and second derivative
from spline fit for each angular momentum and each projector;
vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
xc_denpos= lowest allowed density (usually for the computation of the XC functionals)
xcccrc=XC core correction cutoff radius (bohr) from psp file
SIDE EFFECTS
pawtab <type(pawtab_type)>=paw tabulated starting data tnvale(vale_mesh%mesh_size)= pseudo valence density (+ nhat in output) vloc_mesh<type(pawrad_type)>= radial mesh for the local potential
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
1759 subroutine pawpsp_calc(core_mesh,epsatm,ffspl,imainmesh,ixc,lnmax,& 1760 & mmax,mqgrid_ff,mqgrid_vl,ncore,nmesh,pawrad,pawtab,pawxcdev,pspversion,& 1761 & qgrid_ff,qgrid_vl,radmesh,tncore,tnvale,tproj,tproj_mesh,usexcnhat,vale_mesh,& 1762 & vloc_mesh,vlocopt,vlocr,vlspl,xcccrc,xclevel,xc_denpos,zion,znucl) 1763 1764 1765 !This section has been created automatically by the script Abilint (TD). 1766 !Do not modify the following lines by hand. 1767 #undef ABI_FUNC 1768 #define ABI_FUNC 'pawpsp_calc' 1769 !End of the abilint section 1770 1771 implicit none 1772 1773 !Arguments ------------------------------------ 1774 !scalars 1775 integer,intent(in) :: imainmesh,ixc,lnmax,mqgrid_ff,mqgrid_vl 1776 integer,intent(in) :: nmesh,pawxcdev,pspversion,usexcnhat,vlocopt 1777 integer,intent(in) ::mmax 1778 integer,intent(in) :: xclevel 1779 real(dp),intent(in) :: xc_denpos,zion,znucl 1780 real(dp),intent(out) :: epsatm,xcccrc 1781 type(pawrad_type),intent(in) :: core_mesh,tproj_mesh,vale_mesh 1782 type(pawrad_type),intent(inout) ::pawrad,vloc_mesh 1783 type(pawtab_type),intent(inout) :: pawtab 1784 !arrays 1785 real(dp),intent(in) :: ncore(core_mesh%mesh_size),tncore(core_mesh%mesh_size) 1786 real(dp),intent(in) :: qgrid_vl(mqgrid_vl),qgrid_ff(mqgrid_ff) 1787 real(dp),intent(inout) :: ffspl(mqgrid_ff,2,lnmax) 1788 real(dp),intent(out) :: vlspl(mqgrid_vl,2) 1789 real(dp),intent(inout) :: tnvale(vale_mesh%mesh_size*pawtab%has_tvale) 1790 real(dp),intent(inout) :: tproj(tproj_mesh%mesh_size,pawtab%basis_size) 1791 real(dp),intent(inout) :: vlocr(vloc_mesh%mesh_size) 1792 type(pawrad_type),intent(in) :: radmesh(nmesh) 1793 1794 !Local variables ------------------------------ 1795 !scalars 1796 integer,parameter :: reduced_mshsz=2501 1797 integer :: ib,il,ilm,ilmn,iln,ir,isnotzero,itest 1798 integer :: j0lmn,jlm,jlmn,jln,klmn,msz,msz1,msz_tmp,mst_tmp,nspden 1799 logical :: has_dij0,non_magnetic_xc,reduced_ncor,reduced_nval,reduced_vloc,testval 1800 real(dp),parameter :: reduced_rstep=0.00025_dp,rm_vloc=20.0_dp 1801 real(dp) :: d2nvdq0,intg,intvh,lstep_tmp,qcore,qq,rstep_tmp,yp1,ypn 1802 character(len=500) :: msg 1803 type(pawang_type) :: pawang_tmp 1804 type(pawrad_type) :: rcore_mesh,rvale_mesh,rvloc_mesh,tproj_mesh_new 1805 !arrays 1806 real(dp) :: tmp_qgrid(1),tmp_q2vq(1) 1807 real(dp),allocatable :: ncorwk(:),nhat(:),nhatwk(:),nwk(:),r2k(:) 1808 real(dp),allocatable :: rtncor(:),rtnval(:),rvlocr(:) 1809 real(dp),allocatable :: vbare(:),vh(:),vhnzc(:) 1810 real(dp),allocatable :: vxc1(:),vxc2(:),work1(:),work1b(:),work2(:),work3(:),work4(:) 1811 logical :: tmp_lmselect(1) 1812 1813 ! ************************************************************************* 1814 1815 !========================================================== 1816 !Perfom tests on meshes 1817 1818 ! initialise logical 1819 non_magnetic_xc=.false. 1820 1821 !Are radial meshes for Phi and Vloc compatibles ? 1822 ! if (vloc_mesh%rmax<pawrad%rmax) then 1823 ! write(msg, '(a,a,a)' )& 1824 !& 'Rmax for Vloc < Rmax !',ch10,& 1825 !& 'Action : check your pseudopotential (increase Vloc meshSize).' 1826 ! MSG_ERROR(msg) 1827 ! end if 1828 1829 !Are mmax and mesh_size for partial waves compatibles ? 1830 if (mmax/=pawtab%partialwave_mesh_size) then 1831 write(msg, '(a,a,a)' )& 1832 & 'mmax /= phi_mesh_size in psp file !',ch10,& 1833 & 'Action: check your pseudopotential file.' 1834 MSG_ERROR(msg) 1835 end if 1836 1837 !Are radial meshes for (t)Ncore and Phi compatibles ? 1838 if (core_mesh%mesh_size<pawtab%mesh_size) then 1839 write(msg, '(a,a,a,a,a)' )& 1840 & 'Mesh size for core density must be equal or larger',ch10,& 1841 & 'than mesh size for PAW augmentation regions !',ch10,& 1842 & 'Action : check your pseudopotential (increase Ncore meshSize).' 1843 MSG_ERROR(msg) 1844 end if 1845 1846 !Are radial meshes for (t)Nvale and Phi compatibles ? 1847 if ((pawtab%has_tvale==1).and.(vale_mesh%rmax<pawrad%rmax)) then 1848 write(msg, '(a,a,a)' )& 1849 & 'Rmax for tNvale < Rmax for Phi !',ch10,& 1850 & 'Action : check your pseudopotential (increase tNvale meshSize).' 1851 MSG_ERROR(msg) 1852 end if 1853 1854 !Is PAW radius included inside radial mesh ? 1855 if (pawtab%rpaw>pawrad%rmax+tol8) then 1856 write(msg, '(a,a,a)' )& 1857 & 'Radius of PAW sphere is outside the radial mesh !',ch10,& 1858 & 'Action: check your pseudopotential file.' 1859 MSG_ERROR(msg) 1860 end if 1861 1862 !Max. radius of mesh for Vloc has to be "small" in order to avoid numeric noise ? 1863 if (vloc_mesh%rmax>rm_vloc) then 1864 msz_tmp=pawrad_ifromr(vloc_mesh,rm_vloc);mst_tmp=vloc_mesh%mesh_type 1865 rstep_tmp=vloc_mesh%rstep;lstep_tmp=vloc_mesh%lstep 1866 call pawrad_free(vloc_mesh) 1867 call pawrad_init(vloc_mesh,mesh_size=msz_tmp,mesh_type=mst_tmp,& 1868 & rstep=rstep_tmp,lstep=lstep_tmp,r_for_intg=rm_vloc) 1869 write(msg, '(a,i4,a)' ) ' Mesh size for Vloc has been set to ', & 1870 & vloc_mesh%mesh_size,' to avoid numerical noise.' 1871 call wrtout(std_out,msg,'COLL') 1872 call wrtout(ab_out,msg,'COLL') 1873 end if 1874 1875 !This test has been disable... MT 2006-25-10 1876 !For Simpson rule, it is better to have odd mesh sizes 1877 !itest=0 1878 !do imsh=1,nmesh 1879 !if (mod(radmesh(imsh)%mesh_size,2)==0.and.radmesh(imsh)%mesh_type==1) itest=1 1880 !end do 1881 !if (itest==1) then 1882 ! write(msg, '(5a)' ) & 1883 !& 'Regular radial meshes should have odd number of points ',ch10,& 1884 !& 'for better accuracy of integration sheme (Simpson rule).',ch10,& 1885 !& 'Althought it''s not compulsory, you should change mesh sizes in psp file.' 1886 ! MSG_WARNING(msg) 1887 !end if 1888 1889 !Test the compatibilty between Rpaw and mesh for (t)Phi 1890 if (pspversion>=3) then 1891 itest=pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw) 1892 ! This test has been disable... MT 2015-02-12 1893 ! if (itest+2>radmesh(imainmesh)%mesh_size) then 1894 ! write(msg, '(9a)' ) & 1895 !& 'Atomic data could produce inaccurate results:',ch10,& 1896 !& 'Wavefunctions and pseudo-wavefunctions should',ch10,& 1897 !& 'be given on a radial mesh larger than the PAW',ch10,& 1898 !& 'spheres (at least 2 additional points) !',ch10,& 1899 !& 'Action: check your pseudopotential file.' 1900 ! MSG_WARNING(msg) 1901 ! end if 1902 if (abs(pawtab%rpaw-radmesh(imainmesh)%rad(itest))<tol8) itest=itest-1 1903 ib=0;isnotzero=0 1904 do while ((isnotzero==0).and.(ib<pawtab%basis_size)) 1905 ib=ib+1;ir=itest 1906 do while ((isnotzero==0).and.(ir<pawtab%mesh_size)) 1907 ir=ir+1;if (abs(pawtab%phi(ir,ib)-pawtab%tphi(ir,ib))>tol8) isnotzero=1 1908 end do 1909 end do 1910 if (isnotzero>0) then 1911 write(msg, '(7a)' )& 1912 & 'Atomic data are inconsistent:',ch10,& 1913 & 'For r>=r_paw, pseudo wavefunctions are not',ch10,& 1914 & 'equal to wave functions (Phi(r)/=tPhi(r)) !',ch10,& 1915 & 'Action: check your pseudopotential file.' 1916 MSG_ERROR(msg) 1917 end if 1918 else 1919 ! For compatibility reasons set PAW radius at the end of mesh (older versions) 1920 if (pawtab%rpaw/=pawrad%rmax) then 1921 msz_tmp=pawrad%mesh_size;mst_tmp=pawrad%mesh_type 1922 rstep_tmp=pawrad%rstep;lstep_tmp=pawrad%lstep 1923 call pawrad_free(pawrad) 1924 call pawrad_init(pawrad,mesh_size=msz_tmp,mesh_type=mst_tmp,rstep=rstep_tmp,lstep=lstep_tmp) 1925 pawtab%rpaw=pawrad%rmax 1926 end if 1927 end if 1928 !If Vloc is a "Vbare" potential, it has to be localized inside PAW spheres 1929 if (vlocopt==0.and.(vloc_mesh%rmax>pawtab%rpaw+tol10)) then 1930 if(vlocr(pawrad_ifromr(vloc_mesh,pawtab%rpaw))>tol10) then 1931 write(msg, '(7a)' )& 1932 & 'Atomic data are inconsistent:',ch10,& 1933 & 'Local potential is a "Vbare" potential',ch10,& 1934 & 'and is not localized inside PAW sphere !',ch10,& 1935 & 'Vbare is set to zero if r>rpaw.' 1936 MSG_WARNING(msg) 1937 do ir=pawrad_ifromr(vloc_mesh,pawtab%rpaw),vloc_mesh%mesh_size 1938 vlocr(ir)=zero 1939 end do 1940 end if 1941 end if 1942 1943 !========================================================== 1944 !Initializations 1945 1946 has_dij0=(allocated(pawtab%dij0)) 1947 1948 !Allocate/initialize some dummy variables 1949 tmp_lmselect(1)=.true. 1950 if (pawxcdev==0) then 1951 pawang_tmp%l_size_max=1;pawang_tmp%angl_size=1;pawang_tmp%ylm_size=1 1952 pawang_tmp%use_ls_ylm=0;pawang_tmp%gnt_option=0;pawang_tmp%ngnt=0;pawang_tmp%nsym=0 1953 LIBPAW_ALLOCATE(pawang_tmp%angwgth,(1)) 1954 pawang_tmp%angwgth(1)=one 1955 LIBPAW_ALLOCATE(pawang_tmp%anginit,(3,1)) 1956 pawang_tmp%anginit(1,1)=one 1957 pawang_tmp%anginit(2:3,1)=zero 1958 LIBPAW_ALLOCATE(pawang_tmp%ylmr,(1,1)) 1959 pawang_tmp%ylmr(1,1)=1._dp/sqrt(four_pi) 1960 LIBPAW_ALLOCATE(pawang_tmp%ylmrgr,(3,1,1)) 1961 pawang_tmp%ylmrgr(1:3,1,1)=zero 1962 end if 1963 1964 !========================================================== 1965 !Compute ffspl(q) (and derivatives) 1966 1967 ffspl=zero 1968 if (mqgrid_ff>0) then 1969 call pawpsp_nl(ffspl,pawtab%indlmn,pawtab%lmn_size,lnmax,mqgrid_ff,qgrid_ff,& 1970 & tproj_mesh,tproj) 1971 end if 1972 1973 !========================================================== 1974 !Compute eventually compensation charge radius (i.e. radius for shape functions) 1975 1976 if (pawtab%shape_type>0.and.pawtab%rshp<1.d-8) then 1977 pawtab%rshp=pawtab%rpaw 1978 else if (pawtab%shape_type==-1) then 1979 ir=pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw)+1;isnotzero=0 1980 do while ((isnotzero==0).and.(ir>1)) 1981 ir=ir-1;il=0 1982 do while ((isnotzero==0).and.(il<pawtab%l_size)) 1983 il=il+1;if (pawtab%shapefunc(ir,il)>tol16) isnotzero=1 1984 end do 1985 end do 1986 ir=min(ir+1,pawrad_ifromr(radmesh(imainmesh),pawtab%rpaw)) 1987 pawtab%rshp=radmesh(imainmesh)%rad(ir) 1988 do il=1,pawtab%l_size 1989 if (pawtab%shapefunc(ir,il)>tol6) then 1990 write(msg, '(a,a,a)' )& 1991 & 'Shape function is not zero at PAW radius !',ch10,& 1992 & 'Action: check your pseudopotential file.' 1993 MSG_ERROR(msg) 1994 end if 1995 end do 1996 end if 1997 1998 !========================================================== 1999 !Compute compensation charge density (nhat) 2000 !Add it to pseudo valence density 2001 2002 if (pawtab%has_tvale==1) then 2003 msz=vale_mesh%mesh_size 2004 LIBPAW_ALLOCATE(nhat,(msz)) 2005 ! A-Has to compute norm of nhat (Int[n-tild_n]) 2006 testval=(abs(tnvale(msz))<tol9) 2007 ! A1-If tnvale is not given with enough points, 2008 ! try to compute it from rhoij0 and tphi 2009 if (.not.testval) then 2010 msz1=pawtab%mesh_size 2011 ! Compute n and tild_n from phi and tphi 2012 LIBPAW_ALLOCATE(work1,(msz1)) 2013 LIBPAW_ALLOCATE(work2,(msz1)) 2014 work1=zero 2015 work2=zero 2016 do jlmn=1,pawtab%lmn_size 2017 j0lmn=jlmn*(jlmn-1)/2;jln=pawtab%indlmn(5,jlmn) 2018 do ilmn=1,jlmn 2019 klmn=j0lmn+ilmn;iln=pawtab%indlmn(5,ilmn) 2020 yp1=two;if (ilmn==jlmn) yp1=one 2021 work1(1:msz1)=work1(1:msz1)+yp1*pawtab%rhoij0(klmn) & 2022 & *pawtab% phi(1:msz1,iln)*pawtab% phi(1:msz1,jln) 2023 work2(1:msz1)=work2(1:msz1)+yp1*pawtab%rhoij0(klmn) & 2024 & *pawtab%tphi(1:msz1,iln)*pawtab%tphi(1:msz1,jln) 2025 end do 2026 end do 2027 ! Spline tnvale onto pawrad if needed 2028 LIBPAW_ALLOCATE(nwk,(msz1)) 2029 if ((vale_mesh%mesh_type/=pawrad%mesh_type).or.(vale_mesh%rstep/=pawrad%rstep).or.& 2030 & (vale_mesh%lstep/=pawrad%lstep)) then 2031 LIBPAW_ALLOCATE(work3,(vale_mesh%mesh_size)) 2032 call bound_deriv(tnvale(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn) 2033 call paw_spline(vale_mesh%rad,tnvale,vale_mesh%mesh_size,yp1,ypn,work3) 2034 call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,tnvale,work3,msz1,pawrad%rad(1:msz1),nwk(1:msz1)) 2035 LIBPAW_DEALLOCATE(work3) 2036 else 2037 nwk(1:msz1)=tnvale(1:msz1) 2038 end if 2039 ! Compare tild_n and tnvale (inside aug. region) 2040 if (maxval(abs((nwk(1:msz1)*four_pi*pawrad%rad(1:msz1)**2)-work2(1:msz1)))<tol6) then 2041 ! If equality then compute Int[n-tild_n] 2042 work1=work1-work2 2043 call simp_gen(qq,work1,pawrad) 2044 qq=qq/four_pi 2045 else 2046 ! If not equality, will use tnvale 2047 testval=.true. 2048 write(msg, '(3a)' ) & 2049 & 'Valence density is not given with enough points',ch10,& 2050 & 'in psp file. Some charge estimations will be coarse.' 2051 MSG_WARNING(msg) 2052 end if 2053 LIBPAW_DEALLOCATE(nwk) 2054 LIBPAW_DEALLOCATE(work1) 2055 LIBPAW_DEALLOCATE(work2) 2056 end if 2057 ! A2-If tnvale is given with enough points, use it 2058 if (testval) then 2059 nhat(1:msz)=tnvale(1:msz)*vale_mesh%rad(1:msz)**2 2060 call simp_gen(qq,nhat,vale_mesh) 2061 qq=zion/four_pi-qq 2062 end if 2063 ! B-Compute nhat and add it to pseudo valence density 2064 call atompaw_shpfun(0,vale_mesh,intg,pawtab,nhat) 2065 nhat(1:msz)=qq*nhat(1:msz) 2066 tnvale(1:msz)=tnvale(1:msz)+nhat(1:msz) 2067 end if 2068 2069 !========================================================== 2070 !If Vloc potential is in "Vbare" format, translate it into VH(tnzc) format 2071 2072 if (vlocopt==0) then 2073 write(msg,'(a)') ' Local potential is in "Vbare" format... ' 2074 call wrtout(ab_out,msg,'COLL') 2075 call wrtout(std_out, msg,'COLL') 2076 msz=core_mesh%mesh_size 2077 LIBPAW_ALLOCATE(r2k,(msz)) 2078 call atompaw_shpfun(0,core_mesh,intg,pawtab,r2k) 2079 r2k(1:msz)=r2k(1:msz)*core_mesh%rad(1:msz)**2 2080 ! Compute VH[4pi.r2.n(r)=4pi.r2.tncore(r)+(Qcore-Z).r2.k(r)] 2081 LIBPAW_ALLOCATE(nwk,(core_mesh%mesh_size)) 2082 LIBPAW_ALLOCATE(vh,(core_mesh%mesh_size)) 2083 if (core_mesh%mesh_type==5) then 2084 nwk(1:msz)=tncore(1:msz)*four_pi*core_mesh%rad(1:msz)**2 2085 call simp_gen(qcore,nwk,core_mesh) 2086 qcore=znucl-zion-qcore 2087 else 2088 nwk(1:msz)=(ncore(1:msz)-tncore(1:msz))*four_pi*core_mesh%rad(1:msz)**2 2089 ib=1 2090 do ir=msz,2,-1 2091 if(abs(nwk(ir))<tol14)ib=ir 2092 end do 2093 call simp_gen(qcore,nwk,core_mesh,r_for_intg=core_mesh%rad(ib)) 2094 nwk(1:msz)=tncore(1:msz)*four_pi*core_mesh%rad(1:msz)**2 2095 end if 2096 nwk(1:msz)=nwk(1:msz)+r2k(1:msz)*(qcore-znucl) 2097 call poisson(nwk,0,core_mesh,vh) 2098 vh(2:msz)=vh(2:msz)/core_mesh%rad(2:msz) 2099 call pawrad_deducer0(vh,msz,core_mesh) 2100 2101 LIBPAW_DEALLOCATE(nwk) 2102 ! Eventually spline Vbare 2103 LIBPAW_ALLOCATE(vbare,(core_mesh%mesh_size)) 2104 if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2105 & (core_mesh%rstep /=vloc_mesh%rstep) .or.& 2106 & (core_mesh%lstep /=vloc_mesh%lstep)) then 2107 msz=core_mesh%mesh_size;if (vloc_mesh%rmax<core_mesh%rmax) msz=pawrad_ifromr(core_mesh,vloc_mesh%rmax) 2108 call bound_deriv(vlocr(1:vloc_mesh%mesh_size),vloc_mesh,vloc_mesh%mesh_size,yp1,ypn) 2109 LIBPAW_ALLOCATE(work1,(vloc_mesh%mesh_size)) 2110 LIBPAW_ALLOCATE(work2,(vloc_mesh%mesh_size)) 2111 call paw_spline(vloc_mesh%rad,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1) 2112 call paw_splint(vloc_mesh%mesh_size,vloc_mesh%rad,vlocr,work1,msz,core_mesh%rad(1:msz),vbare) 2113 LIBPAW_DEALLOCATE(work1) 2114 LIBPAW_DEALLOCATE(work2) 2115 else 2116 msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size) 2117 vbare(1:msz)=vlocr(1:msz) 2118 end if 2119 ! Build VH(tnzc) from Vbare 2120 vlocr(1:msz)=vbare(1:msz)+vh(1:msz) 2121 if(vloc_mesh%mesh_size>msz)then 2122 vlocr(msz+1:vloc_mesh%mesh_size)=vh(msz)*vloc_mesh%rad(msz)/vloc_mesh%rad(msz+1:vloc_mesh%mesh_size) 2123 end if 2124 LIBPAW_DEALLOCATE(vbare) 2125 LIBPAW_DEALLOCATE(vh) 2126 2127 ! Compute <tPhi_i|VH(tnzc)|tPhi_j> and int[VH(tnzc)*Qijhat(r)dr] parts of Dij0 2128 ! Note: it is possible as core_mesh and radmesh(imainmesh) have the same steps 2129 if (has_dij0) then 2130 msz=radmesh(imainmesh)%mesh_size 2131 LIBPAW_ALLOCATE(work1,(msz)) 2132 work1(1:msz)=vlocr(1:msz)*r2k(1:msz) 2133 call simp_gen(intvh,work1,radmesh(imainmesh)) 2134 do jlmn=1,pawtab%lmn_size 2135 j0lmn=jlmn*(jlmn-1)/2;jlm=pawtab%indlmn(4,jlmn);jln=pawtab%indlmn(5,jlmn) 2136 do ilmn=1,jlmn 2137 klmn=j0lmn+ilmn;ilm=pawtab%indlmn(4,ilmn);iln=pawtab%indlmn(5,ilmn) 2138 if (jlm==ilm) then 2139 work1(1:msz)=pawtab%tphi(1:msz,iln)*pawtab%tphi(1:msz,jln)*(vlocr(1:msz)-intvh) & 2140 & -pawtab%phi (1:msz,iln)*pawtab%phi (1:msz,jln)*intvh 2141 call simp_gen(intg,work1,radmesh(imainmesh)) 2142 pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg 2143 end if 2144 end do 2145 end do 2146 LIBPAW_DEALLOCATE(work1) 2147 end if 2148 LIBPAW_DEALLOCATE(r2k) 2149 end if 2150 2151 !========================================================== 2152 !If usexcnhat in psp file is different from usexcnhat chosen 2153 !by user, convert VH(tnzc) and Dij0 2154 2155 if (pawtab%usexcnhat==-1) then 2156 pawtab%usexcnhat=usexcnhat 2157 else if (usexcnhat/=pawtab%usexcnhat) then 2158 if (pawtab%has_tvale==0) then 2159 write(msg, '(5a)' ) & 2160 & 'It is only possible to modify the use of compensation charge density',ch10,& 2161 & 'for a file format greater or equal than paw4 !',ch10,& 2162 & 'Action: use usexcnhat=-1 in input file or change psp file format.' 2163 MSG_WARNING(msg) 2164 else 2165 msz=vloc_mesh%mesh_size 2166 ! Retrieve tvale and nhat onto vloc mesh 2167 LIBPAW_ALLOCATE(nwk,(msz)) 2168 LIBPAW_ALLOCATE(ncorwk,(msz)) 2169 LIBPAW_ALLOCATE(nhatwk,(msz)) 2170 nwk=zero;ncorwk=zero;nhatwk=zero 2171 if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2172 & (core_mesh%rstep /=vloc_mesh%rstep) .or.& 2173 & (core_mesh%lstep /=vloc_mesh%lstep)) then 2174 LIBPAW_ALLOCATE(work1,(core_mesh%mesh_size)) 2175 msz1=msz;if (core_mesh%rmax<vloc_mesh%rmax) msz1=pawrad_ifromr(vloc_mesh,core_mesh%rmax) 2176 call bound_deriv(tncore(1:core_mesh%mesh_size),core_mesh,core_mesh%mesh_size,yp1,ypn) 2177 call paw_spline(core_mesh%rad,tncore,core_mesh%mesh_size,yp1,ypn,work1) 2178 call paw_splint(core_mesh%mesh_size,core_mesh%rad,tncore,work1,msz1,vloc_mesh%rad(1:msz1),ncorwk(1:msz1)) 2179 LIBPAW_DEALLOCATE(work1) 2180 else 2181 msz1=min(core_mesh%mesh_size,msz) 2182 ncorwk(1:msz1)=tncore(1:msz1) 2183 end if 2184 if ((vale_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2185 & (vale_mesh%rstep /=vloc_mesh%rstep) .or.& 2186 & (vale_mesh%lstep /=vloc_mesh%lstep)) then 2187 LIBPAW_ALLOCATE(work1,(vale_mesh%mesh_size)) 2188 msz1=msz;if (vale_mesh%rmax<vloc_mesh%rmax) msz1=pawrad_ifromr(vloc_mesh,vale_mesh%rmax) 2189 call bound_deriv(tnvale(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn) 2190 call paw_spline(vale_mesh%rad,tnvale,vale_mesh%mesh_size,yp1,ypn,work1) 2191 call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,tnvale,work1,msz1,vloc_mesh%rad(1:msz1),nwk(1:msz1)) 2192 call bound_deriv(nhat(1:vale_mesh%mesh_size),vale_mesh,vale_mesh%mesh_size,yp1,ypn) 2193 call paw_spline(vale_mesh%rad,nhat,vale_mesh%mesh_size,yp1,ypn,work1) 2194 call paw_splint(vale_mesh%mesh_size,vale_mesh%rad,nhat,work1,msz1,vloc_mesh%rad(1:msz1),nhatwk(1:msz1)) 2195 LIBPAW_DEALLOCATE(work1) 2196 else 2197 msz1=min(vale_mesh%mesh_size,msz) 2198 nwk (1:msz1)=tnvale(1:msz1) 2199 nhatwk(1:msz1)=nhat (1:msz1) 2200 end if 2201 nwk=nwk-nhatwk 2202 nwk=sqrt(four_pi)*nwk;nhatwk=sqrt(four_pi)*nhatwk ! 0th-order moment of densities 2203 2204 ! Compute Vxc without nhat (vxc1) and with nhat (vxc2) 2205 nspden=1 2206 #if defined LIBPAW_HAVE_LIBXC 2207 if (ixc<0) nspden=libxc_functionals_nspin() 2208 #endif 2209 if (ixc<0) then 2210 LIBPAW_ALLOCATE(vxc1,(msz*nspden)) 2211 LIBPAW_ALLOCATE(vxc2,(msz*nspden)) 2212 LIBPAW_ALLOCATE(work1,(msz)) 2213 LIBPAW_ALLOCATE(work1b,(msz)) 2214 LIBPAW_ALLOCATE(work2,(msz*nspden)) 2215 LIBPAW_ALLOCATE(work3,(msz*nspden)) 2216 work2(1:msz)=nwk 2217 work3(1:msz)=nhatwk 2218 if (nspden==2) work2(msz+1:2*msz)=half*nwk 2219 if (nspden==2) work3(msz+1:2*msz)=half*nhatwk 2220 if (pawxcdev/=0) then 2221 call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,work3,0,non_magnetic_xc,msz,nspden,5,& 2222 & pawang_tmp,vloc_mesh,pawxcdev,work2,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2223 call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,work3,0,non_magnetic_xc,msz,nspden,5,& 2224 & pawang_tmp,vloc_mesh,pawxcdev,work2,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2225 vxc1=vxc1/sqrt(four_pi);vxc2=vxc2/sqrt(four_pi) ! Deduce Vxc from its first moment 2226 else 2227 call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,work3,0,0,non_magnetic_xc,msz,nspden,5,& 2228 & pawang_tmp,vloc_mesh,work2,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2229 call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,work3,0,0,non_magnetic_xc,msz,nspden,5,& 2230 & pawang_tmp,vloc_mesh,work2,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2231 end if 2232 LIBPAW_DEALLOCATE(nwk) 2233 LIBPAW_DEALLOCATE(ncorwk) 2234 LIBPAW_DEALLOCATE(nhatwk) 2235 LIBPAW_DEALLOCATE(work1) 2236 LIBPAW_DEALLOCATE(work1b) 2237 LIBPAW_DEALLOCATE(work2) 2238 LIBPAW_DEALLOCATE(work3) 2239 else 2240 LIBPAW_ALLOCATE(vxc1,(msz)) 2241 LIBPAW_ALLOCATE(vxc2,(msz)) 2242 LIBPAW_ALLOCATE(work1,(msz)) 2243 LIBPAW_ALLOCATE(work1b,(msz)) 2244 if (pawxcdev/=0) then 2245 call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,nhatwk,0,non_magnetic_xc,msz,1,5,& 2246 & pawang_tmp,vloc_mesh,pawxcdev,nwk,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2247 call pawxcm(ncorwk,yp1,ypn,0,ixc,work1,1,tmp_lmselect,nhatwk,0,non_magnetic_xc,msz,1,5,& 2248 & pawang_tmp,vloc_mesh,pawxcdev,nwk,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2249 vxc1=vxc1/sqrt(four_pi);vxc2=vxc2/sqrt(four_pi) ! Deduce Vxc from its first moment 2250 else 2251 call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,nhatwk,0,0,non_magnetic_xc,msz,1,5,& 2252 & pawang_tmp,vloc_mesh,nwk,pawtab%usetcore,0,vxc1,xclevel,xc_denpos) 2253 call pawxc(ncorwk,yp1,ypn,ixc,work1,work1b,1,tmp_lmselect,nhatwk,0,0,non_magnetic_xc,msz,1,5,& 2254 & pawang_tmp,vloc_mesh,nwk,pawtab%usetcore,2,vxc2,xclevel,xc_denpos) 2255 end if 2256 LIBPAW_DEALLOCATE(nwk) 2257 LIBPAW_DEALLOCATE(ncorwk) 2258 LIBPAW_DEALLOCATE(nhatwk) 2259 LIBPAW_DEALLOCATE(work1) 2260 LIBPAW_DEALLOCATE(work1b) 2261 endif 2262 ! Compute difference of XC potentials 2263 if (usexcnhat==0.and.pawtab%usexcnhat/=0) vxc1(1:msz)=vxc2(1:msz)-vxc1(1:msz) 2264 if (usexcnhat/=0.and.pawtab%usexcnhat==0) vxc1(1:msz)=vxc1(1:msz)-vxc2(1:msz) 2265 ! Modify VH(tnzc) 2266 vlocr(1:msz)=vlocr(1:msz)-vxc1(1:msz) 2267 if (has_dij0) then 2268 ! Modify Dij0 2269 LIBPAW_ALLOCATE(work2,(pawtab%lmn2_size)) 2270 call atompaw_kij(pawtab%indlmn,work2,pawtab%lmn_size,ncore,0,0,pawtab,pawrad,& 2271 & core_mesh,vloc_mesh,vxc1(1:msz),znucl) 2272 pawtab%dij0=work2 2273 LIBPAW_DEALLOCATE(work2) 2274 end if 2275 LIBPAW_DEALLOCATE(vxc1) 2276 LIBPAW_DEALLOCATE(vxc2) 2277 end if ! has_tvale/=0 2278 end if 2279 if (pawtab%usexcnhat==0) then 2280 write(msg,'(a)') & 2281 & ' Compensation charge density is not taken into account in XC energy/potential' 2282 call wrtout(ab_out,msg,'COLL') 2283 call wrtout(std_out, msg,'COLL') 2284 end if 2285 if (pawtab%usexcnhat==1) then 2286 write(msg,'(a)') & 2287 & ' Compensation charge density is taken into account in XC energy/potential' 2288 call wrtout(ab_out,msg,'COLL') 2289 call wrtout(std_out, msg,'COLL') 2290 end if 2291 2292 !========================================================== 2293 ! Calculate the coefficient beta = \int { vH[nZc](r) - vloc(r) } 4pi r^2 dr 2294 ! 2295 LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size)) 2296 LIBPAW_ALLOCATE(nwk,(core_mesh%mesh_size)) 2297 ! get vH[nZc] 2298 call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl) 2299 2300 !Transpose vlocr mesh into core mesh 2301 nwk(:)=zero 2302 if ((core_mesh%mesh_type/=vloc_mesh%mesh_type).or.& 2303 & (core_mesh%rstep /=vloc_mesh%rstep) .or.& 2304 & (core_mesh%lstep /=vloc_mesh%lstep)) then 2305 msz=core_mesh%mesh_size;if (vloc_mesh%rmax<core_mesh%rmax) msz=pawrad_ifromr(core_mesh,vloc_mesh%rmax) 2306 call bound_deriv(vlocr(1:vloc_mesh%mesh_size),vloc_mesh,vloc_mesh%mesh_size,yp1,ypn) 2307 LIBPAW_ALLOCATE(work1,(vloc_mesh%mesh_size)) 2308 LIBPAW_ALLOCATE(work2,(vloc_mesh%mesh_size)) 2309 call paw_spline(vloc_mesh%rad,vlocr,vloc_mesh%mesh_size,yp1,ypn,work1) 2310 call paw_splint(vloc_mesh%mesh_size,vloc_mesh%rad,vlocr,work1,msz,core_mesh%rad(1:msz),nwk) 2311 LIBPAW_DEALLOCATE(work1) 2312 LIBPAW_DEALLOCATE(work2) 2313 else 2314 msz=min(core_mesh%mesh_size,vloc_mesh%mesh_size) 2315 nwk(1:msz)=vlocr(1:msz) 2316 end if 2317 2318 !Difference 2319 nwk(1:msz)=vhnzc(1:msz)-nwk(1:msz) 2320 if (msz<core_mesh%mesh_size) nwk(msz+1:core_mesh%mesh_size)=zero 2321 2322 !Perform the spherical integration 2323 nwk(1:msz)=nwk(1:msz)*four_pi*core_mesh%rad(1:msz)**2 2324 2325 call simp_gen(pawtab%beta,nwk,core_mesh) 2326 2327 LIBPAW_DEALLOCATE(vhnzc) 2328 LIBPAW_DEALLOCATE(nwk) 2329 2330 write(msg,'(a,e18.6)') & 2331 & ' beta integral value: ',pawtab%beta 2332 call wrtout(std_out,msg,'COLL') 2333 2334 2335 !========================================================== 2336 !Try to optimize CPU time: 2337 !If Vloc mesh size is big, spline Vloc into a smaller log. mesh 2338 2339 reduced_vloc=(vloc_mesh%mesh_size>int(reduced_mshsz)) 2340 if (reduced_vloc) then 2341 msz=vloc_mesh%mesh_size 2342 lstep_tmp=log(0.9999999_dp*vloc_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2343 call pawrad_init(rvloc_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2344 & rstep=reduced_rstep,lstep=lstep_tmp) 2345 LIBPAW_ALLOCATE(rvlocr,(reduced_mshsz)) 2346 call bound_deriv(vlocr(1:msz),vloc_mesh,msz,yp1,ypn) 2347 LIBPAW_ALLOCATE(work1,(msz)) 2348 LIBPAW_ALLOCATE(work2,(msz)) 2349 LIBPAW_ALLOCATE(work3,(msz)) 2350 work3(1:msz)=vloc_mesh%rad(1:msz) 2351 call paw_spline(work3,vlocr,msz,yp1,ypn,work1) 2352 call paw_splint(msz,work3,vlocr,work1,reduced_mshsz,rvloc_mesh%rad,rvlocr) 2353 LIBPAW_DEALLOCATE(work1) 2354 LIBPAW_DEALLOCATE(work2) 2355 LIBPAW_DEALLOCATE(work3) 2356 end if 2357 2358 !Keep VH(tnZc) eventually in memory 2359 if (pawtab%has_vhtnzc==1) then 2360 LIBPAW_ALLOCATE(pawtab%vhtnzc,(pawtab%mesh_size)) 2361 if ((reduced_vloc).and.(rvloc_mesh%mesh_type==pawrad%mesh_type)& 2362 & .and.(rvloc_mesh%rstep==pawrad%rstep).and.(rvloc_mesh%lstep==pawrad%lstep)) then 2363 pawtab%vhtnzc(1:pawtab%mesh_size)=rvlocr(1:pawtab%mesh_size) 2364 pawtab%has_vhtnzc=2 2365 else if ((vloc_mesh%mesh_type==pawrad%mesh_type)& 2366 & .and.(vloc_mesh%rstep==pawrad%rstep).and.(vloc_mesh%lstep==pawrad%lstep)) then 2367 pawtab%vhtnzc(1:pawtab%mesh_size)=vlocr(1:pawtab%mesh_size) 2368 pawtab%has_vhtnzc=2 2369 else 2370 msg = 'Vloc mesh is not right !' 2371 MSG_ERROR(msg) 2372 end if 2373 end if 2374 2375 !========================================================== 2376 !Try to optimize CPU time: 2377 !If ncore mesh size is big, spline tncore into a smaller log. mesh 2378 2379 reduced_ncor=(core_mesh%mesh_size>int(reduced_mshsz)).and.(pawtab%usetcore/=0) 2380 if (reduced_ncor) then 2381 msz=core_mesh%mesh_size 2382 lstep_tmp=log(0.9999999_dp*core_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2383 call pawrad_init(rcore_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2384 & rstep=reduced_rstep,lstep=lstep_tmp) 2385 LIBPAW_ALLOCATE(rtncor,(reduced_mshsz)) 2386 call bound_deriv(tncore(1:msz),core_mesh,msz,yp1,ypn) 2387 LIBPAW_ALLOCATE(work1,(msz)) 2388 LIBPAW_ALLOCATE(work2,(msz)) 2389 LIBPAW_ALLOCATE(work3,(msz)) 2390 work3(1:msz)=core_mesh%rad(1:msz) 2391 call paw_spline(work3,tncore,msz,yp1,ypn,work1) 2392 call paw_splint(msz,work3,tncore,work1,reduced_mshsz,rcore_mesh%rad,rtncor) 2393 LIBPAW_DEALLOCATE(work1) 2394 LIBPAW_DEALLOCATE(work2) 2395 LIBPAW_DEALLOCATE(work3) 2396 end if 2397 2398 !========================================================== 2399 !Try to optimize CPU time: 2400 !If vale mesh size is big, spline tnvale into a smaller log. mesh 2401 2402 if (pawtab%has_tvale==1) then 2403 reduced_nval=(vale_mesh%mesh_size>int(reduced_mshsz)) 2404 if (reduced_nval) then 2405 msz=vale_mesh%mesh_size 2406 lstep_tmp=log(0.9999999_dp*vale_mesh%rmax/reduced_rstep)/dble(reduced_mshsz-2) 2407 call pawrad_init(rvale_mesh,mesh_size=reduced_mshsz,mesh_type=3,& 2408 & rstep=reduced_rstep,lstep=lstep_tmp) 2409 LIBPAW_ALLOCATE(rtnval,(reduced_mshsz)) 2410 call bound_deriv(tnvale(1:msz),vale_mesh,msz,yp1,ypn) 2411 LIBPAW_ALLOCATE(work1,(msz)) 2412 LIBPAW_ALLOCATE(work2,(msz)) 2413 LIBPAW_ALLOCATE(work3,(msz)) 2414 work3(1:msz)=vale_mesh%rad(1:msz) 2415 call paw_spline(work3,tnvale,msz,yp1,ypn,work1) 2416 call paw_splint(msz,work3,tnvale,work1,reduced_mshsz,rvale_mesh%rad,rtnval) 2417 LIBPAW_DEALLOCATE(work1) 2418 LIBPAW_DEALLOCATE(work2) 2419 LIBPAW_DEALLOCATE(work3) 2420 end if 2421 else 2422 reduced_nval=.false. 2423 end if 2424 2425 !========================================================== 2426 !Compute Vlspl(q) (and second derivative) from Vloc(r) 2427 2428 !Compute Vlspl(q)=q^2.Vloc(q) from vloc(r) 2429 if(mqgrid_vl>0) then 2430 if (reduced_vloc) then 2431 call pawpsp_lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),rvloc_mesh,rvlocr,yp1,ypn,zion) 2432 else 2433 call pawpsp_lo(epsatm,mqgrid_vl,qgrid_vl,vlspl(:,1),vloc_mesh,vlocr,yp1,ypn,zion) 2434 end if 2435 ! Compute second derivative of Vlspl(q) 2436 call paw_spline(qgrid_vl,vlspl(:,1),mqgrid_vl,yp1,ypn,vlspl(:,2)) 2437 else 2438 ! Only to compute epsatm 2439 epsatm=zero 2440 if (reduced_vloc) then 2441 call pawpsp_lo(epsatm,1,tmp_qgrid,tmp_q2vq,rvloc_mesh,rvlocr,yp1,ypn,zion) 2442 else 2443 call pawpsp_lo(epsatm,1,tmp_qgrid,tmp_q2vq,vloc_mesh,vlocr,yp1,ypn,zion) 2444 end if 2445 end if 2446 2447 !========================================================== 2448 !Compute tcorespl(q) (and second derivative) from tNcore(r) 2449 2450 pawtab%mqgrid=mqgrid_vl 2451 xcccrc=core_mesh%rmax 2452 LIBPAW_ALLOCATE(pawtab%tcorespl,(pawtab%mqgrid,2)) 2453 2454 if(mqgrid_vl>0.and.pawtab%usetcore/=0) then 2455 ! Compute tcorespl(q)=tNc(q) from tNcore(r) 2456 if (reduced_ncor) then 2457 call pawpsp_cg(pawtab%dncdq0,pawtab%d2ncdq0,mqgrid_vl,qgrid_vl,pawtab%tcorespl(:,1),rcore_mesh,rtncor,yp1,ypn) 2458 else 2459 call pawpsp_cg(pawtab%dncdq0,pawtab%d2ncdq0,mqgrid_vl,qgrid_vl,pawtab%tcorespl(:,1),core_mesh,tncore,yp1,ypn) 2460 end if 2461 ! Compute second derivative of tcorespl(q) 2462 call paw_spline(qgrid_vl,pawtab%tcorespl(:,1),mqgrid_vl,yp1,ypn,pawtab%tcorespl(:,2)) 2463 else 2464 pawtab%tcorespl=zero 2465 pawtab%dncdq0=zero 2466 pawtab%d2ncdq0=zero 2467 end if 2468 2469 !========================================================== 2470 !Compute tvalespl(q) (and second derivative) from tNvale(r) 2471 2472 if (pawtab%has_tvale/=0.and.mqgrid_vl>0) then 2473 LIBPAW_ALLOCATE(pawtab%tvalespl,(pawtab%mqgrid,2)) 2474 if (reduced_nval) then 2475 call pawpsp_cg(pawtab%dnvdq0,d2nvdq0,mqgrid_vl,qgrid_vl,pawtab%tvalespl(:,1),rvale_mesh,rtnval,yp1,ypn) 2476 pawtab%tnvale_mesh_size=rvale_mesh%mesh_size 2477 else 2478 call pawpsp_cg(pawtab%dnvdq0,d2nvdq0,mqgrid_vl,qgrid_vl,pawtab%tvalespl(:,1),vale_mesh,tnvale,yp1,ypn) 2479 pawtab%tnvale_mesh_size=vale_mesh%mesh_size 2480 end if 2481 ! Compute second derivative of tvalespl(q) 2482 call paw_spline(qgrid_vl,pawtab%tvalespl(:,1),mqgrid_vl,yp1,ypn,pawtab%tvalespl(:,2)) 2483 else 2484 pawtab%dnvdq0=zero 2485 pawtab%tnvale_mesh_size=0 2486 end if 2487 2488 !================================================== 2489 !Compute Ex-correlation energy for the core density 2490 2491 nspden=1 2492 #if defined LIBPAW_HAVE_LIBXC 2493 if (ixc<0) nspden=libxc_functionals_nspin() 2494 #endif 2495 2496 LIBPAW_ALLOCATE(work1,(core_mesh%mesh_size*nspden)) 2497 LIBPAW_ALLOCATE(work1b,(core_mesh%mesh_size*nspden)) 2498 LIBPAW_ALLOCATE(work2,(core_mesh%mesh_size*nspden)) 2499 LIBPAW_ALLOCATE(work3,(1)) 2500 LIBPAW_ALLOCATE(work4,(core_mesh%mesh_size)) 2501 work1(:)=zero 2502 if (pawxcdev/=0) then 2503 call pawxcm(ncore,pawtab%exccore,yp1,0,ixc,work4,1,tmp_lmselect,work3,0,non_magnetic_xc,core_mesh%mesh_size,& 2504 & nspden,4,pawang_tmp,core_mesh,pawxcdev,work1,1,0,work2,xclevel,xc_denpos) 2505 else 2506 call pawxc(ncore,pawtab%exccore,yp1,ixc,work4,work1b,1,tmp_lmselect,work3,0,0,non_magnetic_xc,core_mesh%mesh_size,& 2507 & nspden,4,pawang_tmp,core_mesh,work1,1,0,work2,xclevel,xc_denpos) 2508 end if 2509 LIBPAW_DEALLOCATE(work1) 2510 LIBPAW_DEALLOCATE(work1b) 2511 LIBPAW_DEALLOCATE(work2) 2512 LIBPAW_DEALLOCATE(work3) 2513 LIBPAW_DEALLOCATE(work4) 2514 2515 !================================================== 2516 !Compute atomic contribution to Dij (Dij0) 2517 !if not already in memory 2518 if ((.not.has_dij0).and.(pawtab%has_kij==2.or.pawtab%has_kij==-1)) then 2519 LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size)) 2520 if (reduced_vloc) then 2521 call atompaw_dij0(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,& 2522 & rvloc_mesh,rvlocr,znucl) 2523 else 2524 call atompaw_dij0(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,pawtab,pawrad,core_mesh,& 2525 & vloc_mesh,vlocr,znucl) 2526 end if 2527 has_dij0=.true. 2528 end if 2529 !================================================== 2530 !Compute kinetic operator contribution to Dij 2531 2532 if (pawtab%has_kij==1.and.has_dij0) then 2533 LIBPAW_ALLOCATE(pawtab%kij,(pawtab%lmn2_size)) 2534 call atompaw_kij(pawtab%indlmn,pawtab%kij,pawtab%lmn_size,ncore,0,1,pawtab,pawrad,core_mesh,& 2535 & vloc_mesh,vlocr,znucl) 2536 pawtab%has_kij=2 2537 end if 2538 2539 !pawtab%has_kij=-1 means that kij does not have to be kept in memory 2540 if (pawtab%has_kij==-1) then 2541 LIBPAW_DEALLOCATE(pawtab%kij) 2542 pawtab%has_kij=0 2543 end if 2544 2545 !========================================================== 2546 !If projectors have to be kept in memory, we need 2547 !them on the main radial mesh (so, spline them if necessary) 2548 2549 if (pawtab%has_tproj>0) then 2550 if ((tproj_mesh%mesh_type/=pawrad%mesh_type).or.& 2551 & (tproj_mesh%rstep /=pawrad%rstep).or.& 2552 & (tproj_mesh%lstep /=pawrad%lstep)) then 2553 ir=pawrad_ifromr(pawrad,tproj_mesh%rmax) 2554 call pawrad_init(tproj_mesh_new,mesh_size=ir,mesh_type=pawrad%mesh_type,& 2555 & rstep=pawrad%rstep,lstep=pawrad%lstep) 2556 LIBPAW_ALLOCATE(pawtab%tproj,(tproj_mesh_new%mesh_size,pawtab%basis_size)) 2557 LIBPAW_ALLOCATE(work1,(tproj_mesh%mesh_size)) 2558 do ib=1,pawtab%basis_size 2559 call bound_deriv(tproj(:,ib),tproj_mesh,tproj_mesh%mesh_size,yp1,ypn) 2560 call paw_spline(tproj_mesh%rad,tproj(:,ib),tproj_mesh%mesh_size,yp1,ypn,work1) 2561 call paw_splint(tproj_mesh%mesh_size,tproj_mesh%rad,tproj(:,ib),work1,& 2562 & tproj_mesh_new%mesh_size,tproj_mesh_new%rad,pawtab%tproj(:,ib)) 2563 end do 2564 LIBPAW_DEALLOCATE(work1) 2565 call pawrad_free(tproj_mesh_new) 2566 else 2567 LIBPAW_ALLOCATE(pawtab%tproj,(tproj_mesh%mesh_size,pawtab%basis_size)) 2568 pawtab%tproj(:,:)=tproj(:,:) 2569 end if 2570 pawtab%has_tproj=2 2571 end if 2572 2573 !========================================================== 2574 !Free temporary allocated space 2575 2576 if (pawtab%has_tvale==1) then 2577 LIBPAW_DEALLOCATE(nhat) 2578 end if 2579 if (reduced_vloc) then 2580 call pawrad_free(rvloc_mesh) 2581 LIBPAW_DEALLOCATE(rvlocr) 2582 end if 2583 if (reduced_ncor) then 2584 call pawrad_free(rcore_mesh) 2585 LIBPAW_DEALLOCATE(rtncor) 2586 end if 2587 if (reduced_nval) then 2588 call pawrad_free(rvale_mesh) 2589 LIBPAW_DEALLOCATE(rtnval) 2590 end if 2591 if (pawxcdev==0) then 2592 LIBPAW_DEALLOCATE(pawang_tmp%angwgth) 2593 LIBPAW_DEALLOCATE(pawang_tmp%anginit) 2594 LIBPAW_DEALLOCATE(pawang_tmp%ylmr) 2595 LIBPAW_DEALLOCATE(pawang_tmp%ylmrgr) 2596 end if 2597 2598 end subroutine pawpsp_calc
m_pawpsp/pawpsp_calc_d5 [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_calc_d5
FUNCTION
Compute the first to the 5th derivatives of a given function in a pawrad mesh
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
2626 subroutine pawpsp_calc_d5(mesh,mesh_size,tcoredens) 2627 2628 2629 !This section has been created automatically by the script Abilint (TD). 2630 !Do not modify the following lines by hand. 2631 #undef ABI_FUNC 2632 #define ABI_FUNC 'pawpsp_calc_d5' 2633 !End of the abilint section 2634 2635 implicit none 2636 2637 !Arguments ------------------------------------ 2638 integer,intent(in) :: mesh_size 2639 type(pawrad_type),intent(in) :: mesh 2640 real(dp),intent(inout) :: tcoredens(mesh_size,6) 2641 2642 !Local variables------------------------------- 2643 integer,parameter :: it=1 !number of steps for smoothing function 2644 logical,parameter :: use_smooth=.true. 2645 2646 ! ************************************************************************* 2647 2648 !calculate first derivative from density, 2649 !and store it 2650 call nderiv_gen(tcoredens(:,2),tcoredens(:,1),mesh) 2651 2652 !get second derivative from density, and store it 2653 call paw_spline(mesh%rad,tcoredens(:,1),mesh_size,& 2654 & zero,zero,tcoredens(:,3)) 2655 2656 !smooth functions, to avoid numerical instabilities 2657 if(use_smooth) then 2658 call paw_smooth(tcoredens(:,2),mesh_size,it) 2659 call paw_smooth(tcoredens(:,3),mesh_size,it) 2660 end if 2661 2662 !get third derivative from first derivative: 2663 call paw_spline(mesh%rad,tcoredens(:,2),mesh_size,& 2664 & zero,zero,tcoredens(:,4)) 2665 2666 !get fourth derivative from second derivative: 2667 call paw_spline(mesh%rad,tcoredens(:,3),mesh_size,& 2668 & zero,zero,tcoredens(:,5)) 2669 2670 !smooth 3rd and 4th order derivatives 2671 if(use_smooth) then 2672 call paw_smooth(tcoredens(:,4),mesh_size,it) 2673 call paw_smooth(tcoredens(:,5),mesh_size,it) 2674 end if 2675 2676 !get fifth derivative from third derivative: 2677 call paw_spline(mesh%rad,tcoredens(:,4),mesh_size,& 2678 & zero,zero,tcoredens(:,6)) 2679 2680 !smooth 5th order derivative 2681 if(use_smooth) then 2682 call paw_smooth(tcoredens(:,6),mesh_size,it) 2683 end if 2684 2685 end subroutine pawpsp_calc_d5
m_pawpsp/pawpsp_cg [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_cg
FUNCTION
Compute sine transform to transform from n(r) to n(q). Computes integrals on (generalized) grid using corrected trapezoidal integration.
INPUTS
mqgrid=number of grid points in q from 0 to qmax. qgrid(mqgrid)=q grid values (bohr**-1). radmesh <type(pawrad_type)>=data containing radial grid information nr(:)=n(r) on radial grid.
OUTPUT
dnqdq0= 1/q dn(q)/dq for q=0 d2nqdq0 = Gives contribution of d2(tNcore(q))/d2q for q=0 compute \int{(16/15)*pi^5*n(r)*r^6* dr} nq(mqgrid)= n(q) = 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 n(r))dr]. yp1,ypn=derivatives of n(q) wrt q at q=0 and q=qmax (needed for spline fitter).
PARENTS
dfpt_eltfrxc,m_pawpsp,m_psps
CHILDREN
SOURCE
492 subroutine pawpsp_cg(dnqdq0,d2nqdq0,mqgrid,qgrid,nq,radmesh,nr,yp1,ypn) 493 494 495 !This section has been created automatically by the script Abilint (TD). 496 !Do not modify the following lines by hand. 497 #undef ABI_FUNC 498 #define ABI_FUNC 'pawpsp_cg' 499 !End of the abilint section 500 501 implicit none 502 503 !Arguments---------------------------------------------------------- 504 !scalars 505 integer,intent(in) :: mqgrid 506 real(dp),intent(out) :: dnqdq0,d2nqdq0,yp1,ypn 507 type(pawrad_type),intent(in) :: radmesh 508 !arrays 509 real(dp),intent(in) :: nr(:) 510 real(dp),intent(in) :: qgrid(mqgrid) 511 real(dp),intent(out) :: nq(mqgrid) 512 513 !Local variables------------------------------- 514 !scalars 515 integer :: iq,ir,mesh_size 516 real(dp) :: aexp,arg,bexp,dn,r0tor1,r1torm,rm,rmtoin 517 logical :: begin_r0 518 !character(len=500) :: msg 519 !arrays 520 real(dp),allocatable :: ff(:),rnr(:) 521 522 ! ************************************************************************* 523 524 mesh_size=min(size(nr),radmesh%mesh_size) 525 LIBPAW_ALLOCATE(ff,(mesh_size)) 526 LIBPAW_ALLOCATE(rnr,(mesh_size)) 527 ff=zero;rnr=zero 528 529 do ir=1,mesh_size 530 rnr(ir)=radmesh%rad(ir)*nr(ir) 531 end do 532 533 !Is mesh beginning with r=0 ? 534 begin_r0=(radmesh%rad(1)<1.d-20) 535 536 !Adjustment of an exponentional at r_max (n_exp(r)=aexp*Exp[-bexp*r]) 537 rm=radmesh%rad(mesh_size) 538 dn=one/(12._dp*radmesh%stepint*radmesh%radfact(mesh_size)) & 539 & *( 3._dp*nr(mesh_size-4) & 540 & -16._dp*nr(mesh_size-3) & 541 & +36._dp*nr(mesh_size-2) & 542 & -48._dp*nr(mesh_size-1) & 543 & +25._dp*nr(mesh_size)) 544 if (dn<0._dp.and. & 545 & abs(radmesh%rad(mesh_size)*nr(mesh_size))>1.d-20) then 546 bexp=-dn/nr(mesh_size) 547 if (bexp * rm > 50._dp) then 548 ! This solves the problem with the weird core charge used in v4[62] in which bexp x rm ~= 10^3 549 !write(msg,"(a,es16.8)")"Tooooo large bexp * rm: ", bexp*rm, ", setting aexp to 0" 550 !MSG_WARNING(msg) 551 bexp=0.001_dp;aexp=zero 552 else 553 aexp=nr(mesh_size)*exp(bexp*rm) 554 if (abs(aexp)<1.d-20) then 555 bexp=0.001_dp;aexp=zero 556 end if 557 end if 558 else 559 bexp=0.001_dp;aexp=zero 560 end if 561 562 !=========================================== 563 !=== Compute n(q) for q=0 separately 564 !=========================================== 565 566 !Integral from 0 to r1 (only if r1<>0) 567 r0tor1=zero 568 if (.not.begin_r0) r0tor1=(rnr(1)*radmesh%rad(1)**2)/3.d0 569 570 !Integral from r1 to rmax 571 do ir=1,mesh_size 572 if (abs(rnr(ir))>1.d-20) ff(ir)=rnr(ir)*radmesh%rad(ir) 573 end do 574 call simp_gen(r1torm,ff,radmesh) 575 576 !Integral from rmax to infinity 577 !This part is approximated using an exponential density aexp*Exp[-bexp*r] 578 !(formulae obtained with mathematica) 579 rmtoin=aexp*exp(-bexp*rm)/bexp**3*(two+two*bexp*rm+bexp*bexp*rm*rm) 580 581 !Some of the three parts 582 nq(1)=four_pi*(r0tor1+r1torm+rmtoin) 583 584 !=========================================== 585 !=== Compute n(q) for other q''s 586 !=========================================== 587 588 !Loop over q values 589 do iq=2,mqgrid 590 arg=two_pi*qgrid(iq) 591 592 ! Integral from 0 to r1 (only if r1<>0) 593 r0tor1=zero;if (.not.begin_r0) & 594 & r0tor1=nr(1)*(sin(arg*radmesh%rad(1))/arg/arg& 595 & -radmesh%rad(1)*cos(arg*radmesh%rad(1))/arg) 596 597 ! Integral from r1 to rmax 598 do ir=1,mesh_size 599 if (abs(rnr(ir))>1.d-20) ff(ir)=sin(arg*radmesh%rad(ir))*rnr(ir) 600 end do 601 call simp_gen(r1torm,ff,radmesh) 602 603 ! Integral from rmax to infinity 604 ! This part is approximated using an exponential density aexp*Exp[-bexp*r] 605 ! (formulae obtained with mathematica) 606 rmtoin=aexp*exp(-bexp*rm)/(arg**2+bexp**2)**2 & 607 & *(arg*(two*bexp+arg**2*rm+bexp**2*rm)*cos(arg*rm) & 608 & +(arg**2*(bexp*rm-one)+bexp**2*(bexp*rm+one))*sin(arg*rm)) 609 610 ! Store q^2 v(q) 611 nq(iq)=two/qgrid(iq)*(r0tor1+r1torm+rmtoin) 612 end do 613 614 !=========================================== 615 !=== Compute derivatives of n(q) 616 !=== at ends of interval 617 !=========================================== 618 619 !yp(0)=zero 620 yp1=zero 621 622 !yp(qmax)=$ 2\int_0^\infty[(-\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r) r n(r) dr]$ 623 arg=two_pi*qgrid(mqgrid) 624 625 !Integral from 0 to r1 (only if r1<>0) 626 r0tor1=zero;if (.not.begin_r0) & 627 & r0tor1=two_pi*nr(1)*(3.d0*radmesh%rad(1)/arg /arg*cos(arg*radmesh%rad(1))+ & 628 & (radmesh%rad(1)**2/arg-3.0d0/arg**3)*sin(arg*radmesh%rad(1))) 629 630 !Integral from r1 to rmax 631 do ir=1,mesh_size 632 if (abs(rnr(ir))>1.d-20) ff(ir)=(two_pi*radmesh%rad(ir)*cos(arg*radmesh%rad(ir)) & 633 & - sin(arg*radmesh%rad(ir))/qgrid(mqgrid)) *rnr(ir) 634 end do 635 call simp_gen(r1torm,ff,radmesh) 636 637 !Integral from rmax to infinity 638 !This part is approximated using an exponential density aexp*Exp[-bexp*r] 639 !(formulae obtained with mathematica) 640 rmtoin=-one/(qgrid(mqgrid)*(arg**2+bexp**2)**3) & 641 & *aexp*exp(-bexp*rm) & 642 & *((arg**5*rm-two_pi*arg**4*qgrid(mqgrid)*rm*(bexp*rm-two) & 643 & +two*arg**3*bexp*(bexp*rm+one)+arg*bexp**3*(bexp*rm+two) & 644 & -four_pi*arg**2*bexp*qgrid(mqgrid)*(bexp**2*rm**2-three) & 645 & -two_pi*bexp**3*qgrid(mqgrid)*(bexp**2*rm**2+two*bexp*rm+two))*cos(arg*rm) & 646 & +(two*arg**2*bexp**3*rm+two_pi*arg**5*qgrid(mqgrid)*rm**2 & 647 & +arg**4*(bexp*rm-one)+bexp**4*(bexp*rm+one) & 648 & +four_pi*arg**3*qgrid(mqgrid)*(bexp**2*rm**2+two*bexp*rm-one) & 649 & +two_pi*arg*bexp**2*qgrid(mqgrid)*(bexp**2*rm**2+four*bexp*rm+6._dp))*sin(arg*rm)) 650 651 !Some of the three parts 652 ypn=two/qgrid(mqgrid)*(r0tor1+r1torm+rmtoin) 653 654 !=========================================== 655 !=== Compute 1/q dn(q)/dq at q=0 656 !=========================================== 657 658 !Integral from 0 to r1 (only if r1<>0) 659 r0tor1=zero 660 if (.not.begin_r0) r0tor1=(rnr(1)*radmesh%rad(1)**4)/5.d0 661 662 !Integral from r1 to rmax 663 do ir=1,mesh_size 664 if (abs(rnr(ir))>1.d-20) ff(ir)=rnr(ir)*radmesh%rad(ir)**3 665 end do 666 call simp_gen(r1torm,ff,radmesh) 667 668 !Integral from rmax to infinity 669 !This part is approximated using an exponential density aexp*Exp[-bexp*r] 670 !(formulae obtained with mathematica) 671 rmtoin=aexp*exp(-bexp*rm)/bexp**5 & 672 & *(24._dp+24._dp*bexp*rm+12._dp*bexp**2*rm**2+four*bexp**3*rm**3+bexp**4*rm**4) 673 674 !Some of the three parts 675 dnqdq0=-(2.d0/3.d0)*two_pi**3*(r0tor1+r1torm+rmtoin) 676 677 LIBPAW_DEALLOCATE(ff) 678 LIBPAW_DEALLOCATE(rnr) 679 680 d2nqdq0 = 1_dp 681 682 end subroutine pawpsp_cg
m_pawpsp/pawpsp_header_type [ Types ]
[ Top ] [ m_pawpsp ] [ Types ]
NAME
pawpsp_header_type
FUNCTION
For PAW, header related data
SOURCE
89 type, public :: pawpsp_header_type 90 91 !Integer scalars 92 integer :: basis_size ! Number of elements of the wf basis ((l,n) quantum numbers) 93 integer :: l_size ! Maximum value of l+1 leading to a non zero Gaunt coefficient 94 integer :: lmn_size ! Number of elements of the paw basis 95 integer :: mesh_size ! Dimension of (main) radial mesh 96 integer :: pawver ! Version number of paw psp format 97 integer :: shape_type ! Type of shape function 98 real(dp) :: rpaw ! Radius for paw spheres 99 real(dp) :: rshp ! Cut-off radius of shape function 100 101 end type pawpsp_header_type
m_pawpsp/pawpsp_lo [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_lo
FUNCTION
Compute sine transform to transform from V(r) to q^2 V(q). Computes integrals on (generalized) grid using corrected trapezoidal integration.
INPUTS
mqgrid=number of grid points in q from 0 to qmax. qgrid(mqgrid)=q grid values (bohr**-1). radmesh <type(pawrad_type)>=data containing radial grid information vloc(:)=V(r) on radial grid. zion=nominal valence charge of atom.
OUTPUT
epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
q2vq(mqgrid)
=q^2 V(q)
= -\frac{Zv}{\pi}
+ q^2 4\pi\int[(\frac{\sin(2\pi q r)}{2\pi q r})(r^2 V(r)+r Zv)dr].
yp1,ypn=derivatives of q^2 V(q) wrt q at q=0 and q=qmax (needed for spline fitter).
PARENTS
m_pawpsp
CHILDREN
SOURCE
315 subroutine pawpsp_lo(epsatm,mqgrid,qgrid,q2vq,radmesh,vloc,yp1,ypn,zion) 316 317 318 !This section has been created automatically by the script Abilint (TD). 319 !Do not modify the following lines by hand. 320 #undef ABI_FUNC 321 #define ABI_FUNC 'pawpsp_lo' 322 !End of the abilint section 323 324 implicit none 325 326 !Arguments---------------------------------------------------------- 327 !scalars 328 integer,intent(in) :: mqgrid 329 real(dp),intent(in) :: zion 330 real(dp),intent(out) :: epsatm,yp1,ypn 331 type(pawrad_type),intent(in) :: radmesh 332 !arrays 333 real(dp),intent(in) :: qgrid(mqgrid) 334 real(dp),intent(in) :: vloc(:) 335 real(dp),intent(out) :: q2vq(mqgrid) 336 337 !Local variables ------------------------------ 338 !scalars 339 integer :: iq,ir,irmax,mesh_size 340 real(dp) :: arg,r0tor1,r1torm,rmtoin 341 logical :: begin_r0 342 !arrays 343 real(dp),allocatable :: ff(:),rvpz(:) 344 345 !************************************************************************ 346 347 mesh_size=size(vloc) 348 irmax=pawrad_ifromr(radmesh,min(20._dp,radmesh%rmax)) 349 irmax=min(irmax,mesh_size) 350 351 !Particular case of a zero potential 352 if (maxval(abs(vloc(1:irmax)))<=1.e-20_dp) then 353 q2vq=zero;yp1=zero;ypn=zero;epsatm=zero 354 return 355 end if 356 357 LIBPAW_ALLOCATE(ff,(mesh_size)) 358 LIBPAW_ALLOCATE(rvpz,(mesh_size)) 359 ff=zero;rvpz=zero 360 361 !Is mesh beginning with r=0 ? 362 begin_r0=(radmesh%rad(1)<1.e-20_dp) 363 364 !Store r.V+Z 365 do ir=1,irmax 366 rvpz(ir)=radmesh%rad(ir)*vloc(ir)+zion 367 end do 368 369 !=========================================== 370 !=== Compute q^2 v(q) for q=0 separately 371 !=========================================== 372 373 !Integral from 0 to r1 (only if r1<>0) 374 r0tor1=zero;if (.not.begin_r0) & 375 & r0tor1=(zion*0.5_dp+radmesh%rad(1)*vloc(1)/3._dp)*radmesh%rad(1)**2 376 377 !Integral from r1 to rmax 378 do ir=1,irmax 379 if (abs(rvpz(ir))>1.e-20_dp) then 380 ff(ir)=radmesh%rad(ir)*rvpz(ir) 381 end if 382 end do 383 384 call simp_gen(r1torm,ff,radmesh) 385 386 !Integral from rmax to infinity 387 !This part is neglected... might be improved. 388 rmtoin=zero 389 390 !Some of the three parts 391 epsatm=four_pi*(r0tor1+r1torm+rmtoin) 392 393 q2vq(1)=-zion/pi 394 395 !=========================================== 396 !=== Compute q^2 v(q) for other q''s 397 !=========================================== 398 399 !Loop over q values 400 do iq=2,mqgrid 401 arg=two_pi*qgrid(iq) 402 403 ! Integral from 0 to r1 (only if r1<>0) 404 r0tor1=zero;if (.not.begin_r0) & 405 & r0tor1=( vloc(1)/arg*sin(arg*radmesh%rad(1)) & 406 & -rvpz(1) *cos(arg*radmesh%rad(1)) +zion )/pi 407 408 ! Integral from r1 to rmax 409 do ir=1,irmax 410 if (abs(rvpz(ir))>1.e-20_dp) ff(ir)=sin(arg*radmesh%rad(ir))*rvpz(ir) 411 end do 412 call simp_gen(r1torm,ff,radmesh) 413 414 ! Integral from rmax to infinity 415 ! This part is neglected... might be improved. 416 rmtoin=zero 417 418 ! Store q^2 v(q) 419 q2vq(iq)=-zion/pi + two*qgrid(iq)*(r0tor1+r1torm+rmtoin) 420 end do 421 422 !=========================================== 423 !=== Compute derivatives of q^2 v(q) 424 !=== at ends of interval 425 !=========================================== 426 427 !yp(0)=zero 428 yp1=zero 429 430 !yp(qmax)=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$ 431 arg=two_pi*qgrid(mqgrid) 432 433 !Integral from 0 to r1 (only if r1<>0) 434 r0tor1=zero;if (.not.begin_r0) & 435 & r0tor1=zion*radmesh%rad(1) *sin(arg*radmesh%rad(1)) & 436 & +three*radmesh%rad(1)*vloc(1)/arg *cos(arg*radmesh%rad(1)) & 437 & +(radmesh%rad(1)**2-one/arg**2)*vloc(1)*sin(arg*radmesh%rad(1)) 438 439 !Integral from r1 to rmax 440 do ir=1,irmax 441 if (abs(rvpz(ir))>1.e-20_dp) ff(ir)=( arg*radmesh%rad(ir)*cos(arg*radmesh%rad(ir)) & 442 & + sin(arg*radmesh%rad(ir))) *rvpz(ir) 443 end do 444 call simp_gen(r1torm,ff,radmesh) 445 446 !Integral from rmax to infinity 447 !This part is neglected... might be improved. 448 rmtoin=zero 449 450 !Some of the three parts 451 ypn=two*(r0tor1+r1torm+rmtoin) 452 453 LIBPAW_DEALLOCATE(ff) 454 LIBPAW_DEALLOCATE(rvpz) 455 456 end subroutine pawpsp_lo
m_pawpsp/pawpsp_main [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_main
FUNCTION
Reads a PAW dataset (atomic data)
INPUTS
filpsp=name of the file containing the PAW dataset
usewvl=1 if we use a wavelet basis, 0 other wise (plane waves)
icoulomb=1 if we use a Poisson routine with wavelets, 0 otherwise
ixc=index of the XC correlation functional
xclevel=type of XC functional (1=LDA, 2=GGA, ...)
pawxcdev=order of the developement of the PAW on-site terms
(0: full calculation, 1: order 1, 2:order 2)
usexcnhat=flag controlling the use of compensation charge (nhat) in XC potential
qgrid_ff=size of the mesh for the sin FFT transform of the non-local projectors (form factors)
(plane waves only, 0 otherwise)
qgrid_vl=size of the mesh for the sin FFT transform of the local potential
(plane waves only, 0 otherwise)
ffspl=sin FFT transform of the non-local projectors (form factors) (plane waves only)
vlspl=sin FFT transform of the local potential (plane waves only)
epsatm=$ 4\pi\int[r^2 (V(r)+\frac{Zv}{r}dr]$.
xcccrc=XC core correction cutoff radius (bohr)
zionpsp=valence of atom as specified in input file
znuclpsp=atomic number of atom as specified in input file
===== Optional arguments for wvl =====
wvl_ngauss
===== Other optional arguments =====
psxml=datastructure containing a XMP PAW dataset
comm_mpi=MPI communicator
xc_denpos=tolerance on density/potential for the calculation of XC potential
(if density<xc_denpos, density=zero)
OUTPUT
pawrad <type(pawrad_type)>=data containing PAW radial grid information pawtab <type(pawtab_type)>=data containing the PAW dataset (partial waves...)
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
4787 subroutine pawpsp_main( & 4788 & pawrad,pawtab,& 4789 & filpsp,usewvl,icoulomb,ixc,xclevel,pawxcdev,usexcnhat,& 4790 & qgrid_ff,qgrid_vl,ffspl,vlspl,epsatm,xcccrc,zionpsp,znuclpsp,& 4791 & wvl_ngauss,psxml,comm_mpi,xc_denpos) 4792 4793 4794 !This section has been created automatically by the script Abilint (TD). 4795 !Do not modify the following lines by hand. 4796 #undef ABI_FUNC 4797 #define ABI_FUNC 'pawpsp_main' 4798 !End of the abilint section 4799 4800 implicit none 4801 4802 !Arguments ------------------------------------ 4803 !scalars 4804 integer,intent(in) :: icoulomb,ixc 4805 integer,intent(in) :: pawxcdev,usewvl,usexcnhat,xclevel 4806 integer,optional,intent(in) :: comm_mpi 4807 real(dp),intent(in):: zionpsp,znuclpsp 4808 real(dp),optional,intent(in) :: xc_denpos 4809 real(dp),intent(out) :: epsatm,xcccrc 4810 character(len=fnlen),intent(in):: filpsp ! name of the psp file 4811 type(pawrad_type),intent(inout) :: pawrad 4812 type(pawtab_type),intent(inout) :: pawtab 4813 type(paw_setup_t),optional,intent(in) :: psxml 4814 !arrays 4815 integer,optional,intent(in) :: wvl_ngauss(2) 4816 real(dp),intent(in) :: qgrid_ff(:),qgrid_vl(:) 4817 real(dp),intent(inout) :: ffspl(:,:,:) 4818 real(dp),intent(out) :: vlspl(:,:) 4819 4820 !Local variables------------------------------- 4821 integer :: has_tproj,has_wvl,ipsp,lmax,lloc,lnmax,mmax,me,mqgrid_ff,mqgrid_vl 4822 integer :: pspcod,pspxc,usexml 4823 real(dp),parameter :: xc_denpos_default=tol14 4824 real(dp) :: my_xc_denpos,r2well,zion,znucl 4825 character(len=500) :: msg 4826 type(pawpsp_header_type) :: pawpsp_header 4827 !arrays 4828 4829 ! ************************************************************************* 4830 4831 !Check consistency of parameters 4832 if (icoulomb/= 0.or.usewvl==1) then 4833 if (.not.present(wvl_ngauss)) then 4834 msg='usewvl==1 or icoulomb/=0: a mandatory argument is missing!' 4835 MSG_BUG(msg) 4836 end if 4837 end if 4838 4839 mqgrid_ff=size(qgrid_ff) 4840 mqgrid_vl=size(qgrid_vl) 4841 lnmax=size(ffspl,3) 4842 if (size(ffspl,1)/=mqgrid_ff.or.size(ffspl,2)/=2) then 4843 msg='invalid sizes for ffspl!' 4844 MSG_BUG(msg) 4845 end if 4846 if (size(vlspl,1)/=mqgrid_vl.or.size(vlspl,2)/=2) then 4847 msg='invalid sizes for vlspl!' 4848 MSG_BUG(msg) 4849 end if 4850 4851 my_xc_denpos=xc_denpos_default;if (present(xc_denpos)) my_xc_denpos=xc_denpos 4852 pawtab%usexcnhat=usexcnhat 4853 me=0;if (present(comm_mpi))me=xmpi_comm_rank(comm_mpi) 4854 4855 has_wvl=0; if (usewvl==1.or.icoulomb/=0) has_wvl=1 4856 has_tproj=0; if (usewvl==1) has_tproj=1 4857 call pawtab_set_flags(pawtab,has_tvale=1,has_wvl=has_wvl,has_tproj=has_tproj) 4858 4859 if(me==0) then 4860 write(msg, '(a,t38,a)' )'- pspatm: opening atomic psp file',trim(filpsp) 4861 call wrtout(ab_out, msg,'COLL') 4862 call wrtout(std_out, msg,'COLL') 4863 4864 ! This checks if file is xml or UPF 4865 ! It sets usexml as well 4866 call pawpsp_check_xml_upf(filpsp) 4867 4868 ! ---------------------------------------------------------------------------- 4869 if (usexml /= 1) then 4870 ! Open the atomic data file, and read the three first lines 4871 open (unit=tmp_unit,file=filpsp,form='formatted',status='old') 4872 rewind (unit=tmp_unit) 4873 ! Read first 3 lines of psp file: 4874 call pawpsp_read_header(tmp_unit,lloc,lmax,mmax,pspcod,& 4875 & pspxc,r2well,zion,znucl) 4876 4877 else if (usexml == 1 .and. present(psxml)) then 4878 write(msg,'(a,a)') & 4879 & '- pawpsp : Reading pseudopotential header in XML form from ', trim(filpsp) 4880 call wrtout(ab_out,msg,'COLL') 4881 call wrtout(std_out, msg,'COLL') 4882 4883 ! Return header information 4884 call pawpsp_read_header_xml(lloc,lmax,pspcod,& 4885 & pspxc,psxml,r2well,zion,znucl) 4886 ! Fill in pawpsp_header object: 4887 call pawpsp_read_pawheader(pawpsp_header%basis_size,& 4888 & lmax,pawpsp_header%lmn_size,& 4889 & pawpsp_header%l_size,pawpsp_header%mesh_size,& 4890 & pawpsp_header%pawver,psxml,& 4891 & pawpsp_header%rpaw,pawpsp_header%rshp,pawpsp_header%shape_type) 4892 end if 4893 4894 ! Check data for consistency against main routine input 4895 call pawpsp_consistency() 4896 4897 ! Read rest of the PSP file 4898 if (pspcod==7) then 4899 ! ABINIT proprietary format 4900 call pawpsp_7in(epsatm,ffspl,icoulomb,ixc,& 4901 & lmax,lnmax,mmax,mqgrid_ff,mqgrid_vl,& 4902 & pawrad,pawtab,pawxcdev,qgrid_ff,qgrid_vl,& 4903 & usewvl,usexcnhat,vlspl,xcccrc,xclevel,my_xc_denpos,zion,znucl) 4904 4905 else if (pspcod==17)then 4906 ! XML format 4907 ipsp=1 4908 call pawpsp_17in(epsatm,ffspl,icoulomb,ipsp,ixc,lmax,& 4909 & lnmax,mmax,mqgrid_ff,mqgrid_vl,pawpsp_header,pawrad,pawtab,& 4910 & pawxcdev,qgrid_ff,qgrid_vl,usewvl,usexcnhat,vlspl,xcccrc,& 4911 & xclevel,my_xc_denpos,zion,znucl) 4912 4913 end if 4914 end if!me==0 4915 4916 close(unit=tmp_unit) 4917 4918 write(msg,'(3a)') ' pawpsp: atomic psp has been read ',& 4919 & ' and splines computed',ch10 4920 call wrtout(ab_out,msg,'COLL') 4921 call wrtout(std_out, msg,'COLL') 4922 4923 !Communicate PAW objects 4924 if(present(comm_mpi)) then 4925 if(xmpi_comm_size(comm_mpi)>1) then 4926 call pawpsp_bcast(comm_mpi,epsatm,ffspl,pawrad,pawtab,vlspl,xcccrc) 4927 end if 4928 end if 4929 4930 !WVL+PAW: 4931 if(icoulomb/=0.or.usewvl==1) then 4932 if(present(comm_mpi))then 4933 call pawpsp_wvl(filpsp,pawrad,pawtab,usewvl,wvl_ngauss,comm_mpi) 4934 else 4935 call pawpsp_wvl(filpsp,pawrad,pawtab,usewvl,wvl_ngauss) 4936 end if 4937 end if 4938 4939 contains
m_pawpsp/pawpsp_nl [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_nl
FUNCTION
Make paw projector form factors f_l(q) for each l
INPUTS
indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn lmnmax=max number of (l,m,n) components lnmax=max number of (l,n) components mqgrid=number of grid points for q grid qgrid(mqgrid)=values at which form factors are returned radmesh <type(pawrad_type)>=data containing radial grid information wfll(:,lnmax)=paw projector on radial grid
OUTPUT
ffspl(mqgrid,2,lnmax)= form factor f_l(q) and second derivative
NOTES
u_l(r) is the paw projector (input as wfll);
j_l(q) is a spherical Bessel function;
f_l(q) = $ \int_0^{rmax}[j_l(2\pi q r) u_l(r) r dr]$
PARENTS
m_pawpsp,pawinit
CHILDREN
SOURCE
141 subroutine pawpsp_nl(ffspl,indlmn,lmnmax,lnmax,mqgrid,qgrid,radmesh,wfll) 142 143 144 !This section has been created automatically by the script Abilint (TD). 145 !Do not modify the following lines by hand. 146 #undef ABI_FUNC 147 #define ABI_FUNC 'pawpsp_nl' 148 !End of the abilint section 149 150 implicit none 151 152 !Arguments ------------------------------------ 153 !scalars 154 integer,intent(in) :: lmnmax,lnmax,mqgrid 155 type(pawrad_type),intent(in) :: radmesh 156 !arrays 157 integer,intent(in) :: indlmn(6,lmnmax) 158 real(dp),intent(in) :: qgrid(mqgrid) 159 real(dp),intent(in) :: wfll(:,:) 160 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) 161 162 !Local variables------------------------------- 163 !scalars 164 integer :: ilmn,iln,iln0,iq,ir,ll,meshsz,mmax 165 real(dp),parameter :: eps=tol14**4,TOLJ=0.001_dp 166 real(dp) :: arg,argn,bes 167 real(dp) :: besp,qr 168 real(dp) :: yp1,ypn 169 character(len=100) :: msg 170 type(pawrad_type) :: tmpmesh 171 !arrays 172 real(dp),allocatable :: ff(:),gg(:),rr(:),rr2(:),rr2wf(:),rrwf(:),work(:) 173 174 !************************************************************************* 175 176 !Is mesh beginning with r=0 ? 177 if (radmesh%rad(1)>tol10) then 178 msg='Radial mesh cannot begin with r<>0!' 179 MSG_BUG(msg) 180 end if 181 182 meshsz=size(wfll,1) 183 if (meshsz>radmesh%mesh_size) then 184 msg='wrong size for wfll!' 185 MSG_BUG(msg) 186 end if 187 188 !Init. temporary arrays and variables 189 LIBPAW_ALLOCATE(ff,(meshsz)) 190 LIBPAW_ALLOCATE(gg,(meshsz)) 191 LIBPAW_ALLOCATE(rr,(meshsz)) 192 LIBPAW_ALLOCATE(rr2,(meshsz)) 193 LIBPAW_ALLOCATE(rrwf,(meshsz)) 194 LIBPAW_ALLOCATE(rr2wf,(meshsz)) 195 LIBPAW_ALLOCATE(work,(mqgrid)) 196 rr(1:meshsz) =radmesh%rad(1:meshsz) 197 rr2(1:meshsz)=two_pi*rr(1:meshsz)*rr(1:meshsz) 198 argn=two_pi*qgrid(mqgrid) 199 mmax=meshsz 200 201 !Loop on (l,n) projectors 202 iln0=0 203 do ilmn=1,lmnmax 204 iln=indlmn(5,ilmn) 205 if(iln>iln0) then 206 iln0=iln;ll=indlmn(1,ilmn) 207 208 ir=meshsz 209 do while (abs(wfll(ir,iln))<eps) 210 ir=ir-1 211 end do 212 mmax=min(ir+1,meshsz) 213 if (mmax/=radmesh%int_meshsz) then 214 call pawrad_init(tmpmesh,mesh_size=meshsz,mesh_type=radmesh%mesh_type, & 215 & rstep=radmesh%rstep,lstep=radmesh%lstep,r_for_intg=rr(mmax)) 216 else 217 call pawrad_copy(radmesh,tmpmesh) 218 end if 219 220 rrwf(:) =rr (:)*wfll(:,iln) 221 rr2wf(:)=rr2(:)*wfll(:,iln) 222 223 ! 1-Compute f_l(0<q<qmax) 224 if (mqgrid>2) then 225 do iq=2,mqgrid-1 226 arg=two_pi*qgrid(iq) 227 do ir=1,mmax 228 qr=arg*rr(ir) 229 call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ) 230 ff(ir)=bes*rrwf(ir) 231 end do 232 call simp_gen(ffspl(iq,1,iln),ff,tmpmesh) 233 end do 234 end if 235 236 ! 2-Compute f_l(q=0) and first derivative 237 ffspl(1,1,iln)=zero;yp1=zero 238 if (ll==0) then 239 call simp_gen(ffspl(1,1,iln),rrwf,tmpmesh) 240 end if 241 if (ll==1) then 242 call simp_gen(yp1,rr2wf,tmpmesh) 243 yp1=yp1*third 244 end if 245 246 ! 3-Compute f_l(q=qmax) and first derivative 247 if (mqgrid>1) then 248 ! if (ll==0.or.ll==1) then 249 do ir=1,mmax 250 qr=argn*rr(ir) 251 call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ) 252 ff(ir)=bes*rrwf(ir) 253 gg(ir)=besp*rr2wf(ir) 254 end do 255 call simp_gen(ffspl(mqgrid,1,iln),ff,tmpmesh) 256 call simp_gen(ypn,gg,tmpmesh) 257 else 258 ypn=yp1 259 end if 260 261 ! 4-Compute second derivative of f_l(q) 262 call paw_spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln)) 263 264 call pawrad_free(tmpmesh) 265 266 ! End loop on (l,n) projectors 267 end if 268 end do 269 270 LIBPAW_DEALLOCATE(ff) 271 LIBPAW_DEALLOCATE(gg) 272 LIBPAW_DEALLOCATE(rr) 273 LIBPAW_DEALLOCATE(rr2) 274 LIBPAW_DEALLOCATE(rrwf) 275 LIBPAW_DEALLOCATE(rr2wf) 276 LIBPAW_DEALLOCATE(work) 277 278 end subroutine pawpsp_nl
m_pawpsp/pawpsp_read [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_read
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
File format of formatted PAW psp input (the 3 first lines have already been read in calling -pspatm- routine) : (1) title (character) line (2) psps%znuclpsp(ipsp), zion, pspdat (3) pspcod, pspxc, lmax, lloc, mmax, r2well (4) psp_version, creatorID (5) basis_size, lmn_size (6) orbitals (for l=1 to basis_size) (7) number_of_meshes For imsh=1 to number_of_meshes (8) mesh_index, mesh_type ,mesh_size, rad_step[, log_step] (9) r_cut(SPH) (10) shape_type, r_shape[, shapefunction arguments] For iln=1 to basis_size (11) comment(character) (12) radial mesh index for phi (13) phi(r) (for ir=1 to phi_meshsz) For iln=1 to basis_size (14) comment(character) (15) radial mesh index for tphi (16) tphi(r) (for ir=1 to phi_mesh_size) For iln=1 to basis_size (17) comment(character) (18) radial mesh index for tproj (19) tproj(r) (for ir=1 to proj_mesh_size) (20) comment(character) (21) radial mesh index for core_density (22) core_density (for ir=1 to core_mesh_size) (23) comment(character) (24) radial mesh index for pseudo_core_density (25) tcore_density (for ir=1 to core_mesh_size) (26) comment(character) (27) Dij0 (for ij=1 to lmn_size*(lmn_size+1)/2) (28) comment(character) (29) Rhoij0 (for ij=1 to lmn_size*(lmn_size+1)/2) (30) comment(character) (31) radial mesh index for Vloc, format of Vloc (0=Vbare, 1=VH(tnzc), 2=VH(tnzc) without nhat in XC) (32) Vloc(r) (for ir=1 to vloc_mesh_size) ===== Following lines only if shape_type=-1 ===== For il=1 to 2*max(orbitals)+1 (33) comment(character) (34) radial mesh index for shapefunc (35) shapefunc(r)*gnorm(l)*r**l (for ir=1 to shape_mesh_size) (36) comment(character) (37) radial mesh index for pseudo_valence_density (38) tvale(r) (for ir=1 to vale_mesh_size) Comments: * psp_version= ID of PAW_psp version 4 characters string of the form 'pawn' (with n varying) * creatorID= ID of psp generator creatorid=1xyz : psp generated from Holzwarth AtomPAW generator version x.yz creatorid=2xyz : psp generated from Vanderbilt ultra-soft generator version x.yz creatorid=-1: psp for tests (for developpers only) * mesh_type= type of radial mesh mesh_type=1 (regular grid): rad(i)=(i-1)*AA mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n * radial shapefunction type shape_type=-1 ; gl(r)=numeric (read from psp file) shape_type= 1 ; gl(r)=k(r).r^l; k(r)=exp[-(r/sigma)**lambda] shape_type= 2 ; gl(r)=k(r).r^l; k(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2 if r<=rshp shape_type= 3 ; gl(r)=Alpha(1,l)*jl(q(1,l)*r)+Alpha(2,l)*jl(q(2,l)*r) for each l
PARENTS
m_pawpsp
CHILDREN
SOURCE
772 subroutine pawpsp_read(core_mesh,funit,imainmesh,lmax,& 773 & ncore,nmesh,pawrad,pawtab,pspversion,radmesh,save_core_msz,& 774 & tncore,tnvale,tproj,tproj_mesh,usexcnhat_in,usexcnhat_out,vale_mesh,& 775 & vlocopt,vlocr,vloc_mesh,znucl) 776 777 778 !This section has been created automatically by the script Abilint (TD). 779 !Do not modify the following lines by hand. 780 #undef ABI_FUNC 781 #define ABI_FUNC 'pawpsp_read' 782 !End of the abilint section 783 784 implicit none 785 786 !Arguments ------------------------------------ 787 integer,intent(in):: funit,lmax,usexcnhat_in 788 integer,intent(out) :: imainmesh,pspversion,usexcnhat_out,vlocopt 789 logical,intent(in) :: save_core_msz 790 real(dp),intent(in):: znucl 791 !arrays 792 real(dp),pointer :: ncore(:),tncore(:),tnvale(:),tproj(:,:),vlocr(:) 793 type(pawrad_type),intent(inout) :: pawrad 794 type(pawrad_type),intent(out)::core_mesh,tproj_mesh,vale_mesh,vloc_mesh 795 type(pawrad_type),pointer :: radmesh(:) 796 type(pawtab_type),intent(inout) :: pawtab 797 integer,intent(out)::nmesh 798 799 !Local variables------------------------------- 800 integer :: creatorid,imsh 801 integer :: icoremesh,ishpfmesh,ivalemesh,ivlocmesh 802 integer :: ib,il,ilm,ilmn,iln,iprojmesh 803 integer :: ii,ir,iread1,iread2,jj 804 integer :: msz,pngau_,ptotgau_ 805 real(dp):: rc,rread1,rread2 806 real(dp) :: yp1,ypn 807 !arrays 808 integer,allocatable :: nprj(:) 809 real(dp),allocatable :: shpf(:,:),val(:),vhnzc(:) 810 real(dp),allocatable :: work1(:),work2(:),work3(:),work4(:) 811 character :: blank=' ',numb=' ' 812 character(len=80) :: pspline 813 character(len=500) :: msg,submsg 814 logical :: read_gauss=.false. 815 type(pawrad_type)::shpf_mesh 816 817 ! ************************************************************************* 818 819 !========================================================== 820 !Read lines 4 to 11 of the header 821 822 !This is important for BigDFT in standalone mode 823 call pawpsp_read_header_2(funit,pspversion,pawtab%basis_size,pawtab%lmn_size) 824 825 !Check pspversion for wvl-paw 826 if(pspversion<4 .and. pawtab%has_wvl>0)then 827 write(msg, '(a,i2,a,a)' )& 828 & 'In reading atomic psp file, finds pspversion=',pspversion,ch10,& 829 & 'For WVL-PAW, pspversion >= 4 is required.' 830 MSG_BUG(msg) 831 end if 832 833 834 !Have to maintain compatibility with Abinit v4.2.x 835 if (pspversion==1) then 836 LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size)) 837 read(funit,*) (pawtab%orbitals(ib), ib=1,pawtab%basis_size) 838 pawtab%l_size=2*maxval(pawtab%orbitals)+1 839 nmesh=3 840 LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh)) 841 read(funit,'(a80)') pspline 842 radmesh(1)%lstep=zero 843 read(unit=pspline,fmt=*,err=10,end=10) radmesh(1)%mesh_type,& 844 & radmesh(1)%rstep,radmesh(1)%lstep 845 10 read(funit,*) pawtab%rpaw 846 read(funit,*) radmesh(1)%mesh_size,radmesh(2)%mesh_size,& 847 & radmesh(3)%mesh_size 848 read(funit,'(a80)') pspline 849 pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99 850 read(unit=pspline,fmt=*,err=11,end=11) pawtab%shape_type,& 851 & pawtab%shape_lambda,pawtab%shape_sigma 852 11 read(funit,*) creatorid 853 if (pawtab%shape_type==3) pawtab%shape_type=-1 854 radmesh(2)%mesh_type=radmesh(1)%mesh_type 855 radmesh(3)%mesh_type=radmesh(1)%mesh_type 856 radmesh(2)%rstep=radmesh(1)%rstep 857 radmesh(3)%rstep=radmesh(1)%rstep 858 radmesh(2)%lstep=radmesh(1)%lstep 859 radmesh(3)%lstep=radmesh(1)%lstep 860 else 861 862 ! Here psp file for Abinit 4.3+ 863 LIBPAW_ALLOCATE(pawtab%orbitals,(pawtab%basis_size)) 864 read(funit,*) (pawtab%orbitals(ib), ib=1,pawtab%basis_size) 865 pawtab%l_size=2*maxval(pawtab%orbitals)+1 866 read(funit,*) nmesh 867 LIBPAW_DATATYPE_ALLOCATE(radmesh,(nmesh)) 868 do imsh=1,nmesh 869 rread2=zero 870 read(funit,'(a80)') pspline 871 read(unit=pspline,fmt=*,err=20,end=20) ii,iread1,iread2,rread1,rread2 872 20 continue 873 if (ii<=nmesh) then 874 radmesh(ii)%mesh_type=iread1 875 radmesh(ii)%mesh_size=iread2 876 radmesh(ii)%rstep=rread1 877 radmesh(ii)%lstep=rread2 878 else 879 write(msg, '(3a)' )& 880 & 'Index of mesh out of range !',ch10,& 881 & 'Action : check your pseudopotential file.' 882 MSG_ERROR(msg) 883 end if 884 end do 885 read(funit,*) pawtab%rpaw 886 read(funit,'(a80)') pspline 887 read(unit=pspline,fmt=*) pawtab%shape_type 888 pawtab%shape_lambda=-1;pawtab%shape_sigma=1.d99 889 end if 890 891 !Initialize radial meshes 892 do imsh=1,nmesh 893 call pawrad_init(radmesh(imsh)) 894 end do 895 896 !========================================================== 897 !Initialize various dims and indexes 898 899 pawtab%l_size=2*maxval(pawtab%orbitals)+1 900 pawtab%lmn2_size=pawtab%lmn_size*(pawtab%lmn_size+1)/2 901 pawtab%ij_size=pawtab%basis_size*(pawtab%basis_size+1)/2 902 pawtab%usexcnhat=usexcnhat_in 903 904 !indlmn calculation (indices for (l,m,n) basis) 905 if (allocated(pawtab%indlmn)) then 906 LIBPAW_DEALLOCATE(pawtab%indlmn) 907 end if 908 LIBPAW_ALLOCATE(pawtab%indlmn,(6,pawtab%lmn_size)) 909 LIBPAW_BOUND1_ALLOCATE(nprj,BOUNDS(0,maxval(pawtab%orbitals))) 910 pawtab%indlmn(:,:)=0 911 ilmn=0;iln=0;nprj=0 912 do ib=1,pawtab%basis_size 913 il=pawtab%orbitals(ib) 914 nprj(il)=nprj(il)+1 915 iln=iln+1 916 do ilm=1,2*il+1 917 pawtab%indlmn(1,ilmn+ilm)=il 918 pawtab%indlmn(2,ilmn+ilm)=ilm-(il+1) 919 pawtab%indlmn(3,ilmn+ilm)=nprj(il) 920 pawtab%indlmn(4,ilmn+ilm)=il*il+ilm 921 pawtab%indlmn(5,ilmn+ilm)=iln 922 pawtab%indlmn(6,ilmn+ilm)=1 923 end do 924 ilmn=ilmn+2*il+1 925 end do 926 LIBPAW_DEALLOCATE(nprj) 927 !Are ilmn (found here) and pawtab%lmn_size compatibles ? 928 if (ilmn/=pawtab%lmn_size) then 929 write(msg, '(a,a,a,a,a)' )& 930 & 'Calculated lmn size differs from',ch10,& 931 & 'lmn_size read from pseudo !',ch10,& 932 & 'Action: check your pseudopotential file.' 933 MSG_ERROR(msg) 934 end if 935 936 !========================================================== 937 !Here reading shapefunction parameters 938 939 !Shapefunction parameters for Abinit 4.3...4.5 940 if (pspversion==2) then 941 if (pawtab%shape_type==1) read(unit=pspline,fmt=*) ii,pawtab%shape_lambda,pawtab%shape_sigma 942 if (pawtab%shape_type==3) pawtab%shape_type=-1 943 pawtab%rshp=zero 944 945 !Shapefunction parameters for Abinit 4.6+ 946 else if (pspversion>=3) then 947 pawtab%rshp=zero 948 if (pawtab%shape_type==-1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp 949 if (pawtab%shape_type== 1) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp, & 950 & pawtab%shape_lambda,pawtab%shape_sigma 951 if (pawtab%shape_type== 2) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp 952 if (pawtab%shape_type== 3) read(unit=pspline,fmt=*,err=21,end=21) ii,pawtab%rshp 953 end if 954 21 continue 955 !If shapefunction type is gaussian, check exponent 956 if (pawtab%shape_type==1) then 957 if (pawtab%shape_lambda<2) then 958 write(msg, '(3a)' )& 959 & 'For a gaussian shape function, exponent lambda must be >1 !',ch10,& 960 & 'Action: check your psp file.' 961 MSG_ERROR(msg) 962 end if 963 end if 964 !If shapefunction type is Bessel, deduce here its parameters from rc 965 if (pawtab%shape_type==3) then 966 LIBPAW_ALLOCATE(pawtab%shape_alpha,(2,pawtab%l_size)) 967 LIBPAW_ALLOCATE(pawtab%shape_q,(2,pawtab%l_size)) 968 rc=pawtab%rshp;if (rc<1.d-8) rc=pawtab%rpaw 969 do il=1,pawtab%l_size 970 call atompaw_shapebes(pawtab%shape_alpha(1:2,il),pawtab%shape_q(1:2,il),il-1,rc) 971 end do 972 end if 973 974 !========================================================== 975 !Mirror pseudopotential parameters to the output and log files 976 977 write(msg,'(a,i1)')' Pseudopotential format is: paw',pspversion 978 call wrtout(ab_out,msg,'COLL') 979 call wrtout(std_out, msg,'COLL') 980 write(msg,'(2(a,i3),a,64i4)') & 981 & ' basis_size (lnmax)=',pawtab%basis_size,' (lmn_size=',& 982 & pawtab%lmn_size,'), orbitals=',pawtab%orbitals(1:pawtab%basis_size) 983 call wrtout(ab_out,msg,'COLL') 984 call wrtout(std_out, msg,'COLL') 985 write(msg,'(a,f11.8)')' Spheres core radius: rc_sph=',pawtab%rpaw 986 call wrtout(ab_out,msg,'COLL') 987 call wrtout(std_out, msg,'COLL') 988 write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:' 989 call wrtout(ab_out,msg,'COLL') 990 call wrtout(std_out, msg,'COLL') 991 do imsh=1,nmesh 992 if (radmesh(imsh)%mesh_type==1) & 993 & write(msg,'(a,i1,a,i4,a,g12.5)') & 994 & ' - mesh ',imsh,': r(i)=step*(i-1), size=',radmesh(imsh)%mesh_size,& 995 & ' , step=',radmesh(imsh)%rstep 996 if (radmesh(imsh)%mesh_type==2) & 997 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 998 & ' - mesh ',imsh,': r(i)=AA*[exp(BB*(i-1))-1], size=',radmesh(imsh)%mesh_size,& 999 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 1000 if (radmesh(imsh)%mesh_type==3) & 1001 & write(msg,'(a,i1,a,i4,2(a,g12.5))') & 1002 & ' - mesh ',imsh,': r(i)=AA*exp(BB*(i-2)), size=',radmesh(imsh)%mesh_size,& 1003 & ' , AA=',radmesh(imsh)%rstep,' BB=',radmesh(imsh)%lstep 1004 if (radmesh(imsh)%mesh_type==4) & 1005 & write(msg,'(a,i1,a,i4,a,g12.5)') & 1006 & ' - mesh ',imsh,': r(i)=-AA*ln(1-(i-1)/n), n=size=',radmesh(imsh)%mesh_size,& 1007 & ' , AA=',radmesh(imsh)%rstep 1008 call wrtout(ab_out,msg,'COLL') 1009 call wrtout(std_out, msg,'COLL') 1010 end do 1011 if (pawtab%shape_type==-1) then 1012 write(msg,'(a)')& 1013 ' Shapefunction is NUMERIC type: directly read from atomic data file' 1014 call wrtout(ab_out,msg,'COLL') 1015 call wrtout(std_out, msg,'COLL') 1016 end if 1017 if (pawtab%shape_type==1) then 1018 write(msg,'(2a,a,f6.3,a,i3)')& 1019 & ' Shapefunction is EXP type: shapef(r)=exp(-(r/sigma)**lambda)',ch10,& 1020 & ' with sigma=',pawtab%shape_sigma,' and lambda=',pawtab%shape_lambda 1021 call wrtout(ab_out,msg,'COLL') 1022 call wrtout(std_out, msg,'COLL') 1023 end if 1024 if (pawtab%shape_type==2) then 1025 write(msg,'(a)')& 1026 ' Shapefunction is SIN type: shapef(r)=[sin(pi*r/rshp)/(pi*r/rshp)]**2' 1027 call wrtout(ab_out,msg,'COLL') 1028 call wrtout(std_out, msg,'COLL') 1029 end if 1030 if (pawtab%shape_type==3) then 1031 write(msg,'(a)')& 1032 & ' Shapefunction is BESSEL type: shapef(r,l)=aa(1,l)*jl(q(1,l)*r)+aa(2,l)*jl(q(2,l)*r)' 1033 call wrtout(ab_out,msg,'COLL') 1034 call wrtout(std_out, msg,'COLL') 1035 end if 1036 if (pawtab%rshp<1.d-8) then 1037 write(msg,'(a)') ' Radius for shape functions = sphere core radius' 1038 else 1039 write(msg,'(a,f11.8)') ' Radius for shape functions = ',pawtab%rshp 1040 end if 1041 call wrtout(ab_out,msg,'COLL') 1042 call wrtout(std_out, msg,'COLL') 1043 1044 !========================================================== 1045 !Perfom tests 1046 1047 !Are lmax and orbitals compatibles ? 1048 if (lmax/=maxval(pawtab%orbitals)) then 1049 write(msg, '(a,a,a)' )& 1050 & 'lmax /= MAX(orbitals) !',ch10,& 1051 & 'Action: check your pseudopotential file.' 1052 MSG_ERROR(msg) 1053 end if 1054 1055 !Only mesh_type=1,2, 3 or 4 allowed 1056 do imsh=1,nmesh 1057 if (radmesh(imsh)%mesh_type>4) then 1058 write(msg, '(a,a,a)' )& 1059 & 'Only mesh types 1,2, 3 or 4 allowed !',ch10,& 1060 & 'Action : check your pseudopotential or input file.' 1061 MSG_ERROR(msg) 1062 end if 1063 end do 1064 1065 !========================================================== 1066 !Read tabulated atomic data 1067 1068 !--------------------------------- 1069 !Read wave-functions (phi) 1070 do ib=1,pawtab%basis_size 1071 read (funit,*) 1072 if (pspversion==1) iread1=1 1073 if (pspversion>1) read (funit,*) iread1 1074 if (ib==1) then 1075 call pawrad_free(pawrad) 1076 call pawrad_init(pawrad,mesh_size=radmesh(iread1)%mesh_size,mesh_type=radmesh(iread1)%mesh_type,& 1077 & rstep=radmesh(iread1)%rstep,lstep=radmesh(iread1)%lstep,r_for_intg=pawtab%rpaw) 1078 pawtab%partialwave_mesh_size=pawrad%mesh_size 1079 pawtab%mesh_size=pawrad_ifromr(pawrad,pawtab%rpaw)+5 1080 pawtab%mesh_size=min(pawtab%mesh_size,pawrad%mesh_size) 1081 if (pawtab%mesh_size>pawrad%mesh_size-2) pawtab%mesh_size=pawrad%mesh_size 1082 imainmesh=iread1 1083 LIBPAW_ALLOCATE(pawtab%phi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 1084 else if (iread1/=imainmesh) then 1085 write(msg, '(a,a,a)' )& 1086 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 1087 & 'Action: check your pseudopotential file.' 1088 MSG_ERROR(msg) 1089 end if 1090 read (funit,*) (pawtab%phi(ir,ib),ir=1,pawtab%partialwave_mesh_size) 1091 end do 1092 1093 !--------------------------------- 1094 !Read pseudo wave-functions (tphi) 1095 LIBPAW_ALLOCATE(pawtab%tphi,(pawtab%partialwave_mesh_size,pawtab%basis_size)) 1096 do ib=1,pawtab%basis_size 1097 read (funit,*) 1098 if (pspversion==1) iread1=1 1099 if (pspversion>1) read (funit,*) iread1 1100 if (iread1/=imainmesh) then 1101 write(msg, '(a,a,a)' )& 1102 & 'All Phi and tPhi must be given on the same radial mesh !',ch10,& 1103 & 'Action: check your pseudopotential file.' 1104 MSG_ERROR(msg) 1105 end if 1106 read (funit,*) (pawtab%tphi(ir,ib),ir=1,pawtab%partialwave_mesh_size) 1107 end do 1108 write(msg,'(a,i1)') & 1109 & ' Radial grid used for partial waves is grid ',imainmesh 1110 call wrtout(ab_out,msg,'COLL') 1111 call wrtout(std_out, msg,'COLL') 1112 1113 !--------------------------------- 1114 !Read projectors (tproj) 1115 do ib=1,pawtab%basis_size 1116 read (funit,*) 1117 if (pspversion==1) iread1=2 1118 if (pspversion>1) read (funit,*) iread1 1119 if (ib==1) then 1120 iprojmesh=iread1 1121 call pawrad_copy(radmesh(iprojmesh),tproj_mesh) 1122 LIBPAW_POINTER_ALLOCATE(tproj,(tproj_mesh%mesh_size,pawtab%basis_size)) 1123 else if (iread1/=iprojmesh) then 1124 write(msg, '(a,a,a)' )& 1125 & 'All tprojectors must be given on the same radial mesh !',ch10,& 1126 & 'Action: check your pseudopotential file.' 1127 MSG_ERROR(msg) 1128 end if 1129 ! read projectors from a mesh 1130 read (funit,*) (tproj(ir,ib),ir=1,tproj_mesh%mesh_size) 1131 end do 1132 write(msg,'(a,i2)') & 1133 & ' Radial grid used for projectors is grid ',iprojmesh 1134 call wrtout(ab_out,msg,'COLL') 1135 call wrtout(std_out, msg,'COLL') 1136 1137 !--------------------------------- 1138 !Read gaussian projectors for wavelets 1139 ! -- only if pawtab%has_wvl flag is on 1140 ! -- if not, we skip the lines 1141 read(funit,'(a80)') pspline 1142 if(index(trim(pspline),'GAUSSIAN')/=0) read_gauss=.true. 1143 if (read_gauss) then 1144 if (pawtab%has_wvl>0) then 1145 call wvlpaw_allocate(pawtab%wvl) 1146 jj=0 1147 do ib=1,pawtab%basis_size 1148 if(ib/=1) read(funit,*) pspline 1149 ! read Gaussian coefficients 1150 read(funit,*) pngau_, ptotgau_ !total number of gaussians 1151 if(ib==1) then 1152 pawtab%wvl%ptotgau=ptotgau_ 1153 LIBPAW_ALLOCATE(pawtab%wvl%pngau,(pawtab%basis_size)) 1154 LIBPAW_ALLOCATE(pawtab%wvl%parg,(2,pawtab%wvl%ptotgau)) 1155 LIBPAW_ALLOCATE(pawtab%wvl%pfac,(2,pawtab%wvl%ptotgau)) 1156 else 1157 if(pawtab%wvl%ptotgau/=ptotgau_) then 1158 write(msg,'(3a)')& 1159 & 'Total number of gaussians, should be the same for all projectors !',ch10,& 1160 & 'Action: check your pseudopotential file.' 1161 MSG_ERROR(msg) 1162 end if 1163 end if !ib==1 1164 read(funit,*)(pawtab%wvl%parg(:,ii),ii=jj+1,jj+pngau_) 1165 read(funit,*)(pawtab%wvl%pfac(:,ii),ii=jj+1,jj+pngau_) 1166 pawtab%wvl%pngau(ib)=pngau_ 1167 jj=jj+pngau_ 1168 end do 1169 pawtab%has_wvl=2 1170 else 1171 ! If pawtab%has_wvl=0, we skip the lines 1172 do ib=1,pawtab%basis_size 1173 if(ib/=1) read(funit,*) 1174 read(funit,*) pngau_, ptotgau_ 1175 LIBPAW_ALLOCATE(val, (pngau_ *2)) 1176 read(funit,*) val 1177 read(funit,*) val 1178 LIBPAW_DEALLOCATE(val) 1179 end do 1180 end if 1181 end if 1182 1183 !--------------------------------- 1184 !Read core density (coredens) 1185 if(read_gauss) read (funit,*) !if not read_gauss, this line was already read 1186 if (pspversion==1) iread1=1 1187 if (pspversion>1) read (funit,*) iread1 1188 icoremesh=iread1 1189 call pawrad_copy(radmesh(icoremesh),core_mesh) 1190 if ((radmesh(icoremesh)%mesh_type/=pawrad%mesh_type).or.& 1191 & (radmesh(icoremesh)%rstep /=pawrad%rstep) .or.& 1192 & (radmesh(icoremesh)%lstep /=pawrad%lstep)) then 1193 write(msg, '(a,a,a,a,a)' )& 1194 & 'Ncore must be given on a radial mesh with the same',ch10,& 1195 & 'type and step(s) than the main radial mesh (mesh for Phi) !',ch10,& 1196 & 'Action: check your pseudopotential file.' 1197 MSG_ERROR(msg) 1198 end if 1199 LIBPAW_POINTER_ALLOCATE(ncore,(core_mesh%mesh_size)) 1200 read (funit,*) (ncore(ir),ir=1,core_mesh%mesh_size) 1201 1202 !Construct and save VH[z_NC] if requested 1203 if (pawtab%has_vhnzc==1) then 1204 LIBPAW_ALLOCATE(pawtab%VHnZC,(pawtab%mesh_size)) 1205 LIBPAW_ALLOCATE(vhnzc,(core_mesh%mesh_size)) 1206 call atompaw_vhnzc(ncore,core_mesh,vhnzc,znucl) 1207 pawtab%VHnZC(1:pawtab%mesh_size)=vhnzc(1:pawtab%mesh_size) 1208 pawtab%has_vhnzc=2 1209 LIBPAW_DEALLOCATE(vhnzc) 1210 end if 1211 1212 pawtab%core_mesh_size=pawrad%mesh_size 1213 if(save_core_msz) pawtab%core_mesh_size=core_mesh%mesh_size 1214 LIBPAW_ALLOCATE(pawtab%coredens,(pawtab%core_mesh_size)) 1215 pawtab%rcore=core_mesh%rad(pawtab%core_mesh_size) 1216 pawtab%coredens(1:pawtab%core_mesh_size)=ncore(1:pawtab%core_mesh_size) 1217 1218 !--------------------------------- 1219 !Read pseudo core density (tcoredens) 1220 if(save_core_msz) then 1221 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,6)) 1222 else 1223 LIBPAW_ALLOCATE(pawtab%tcoredens,(pawtab%core_mesh_size,1)) 1224 end if 1225 pawtab%tcoredens=zero 1226 read (funit,*) 1227 if (pspversion==1) iread1=1 1228 if (pspversion>1) read (funit,*) iread1 1229 if (iread1/=icoremesh) then 1230 write(msg, '(a,a,a,a,a,a,a,a)' )& 1231 & 'Pseudized core density (tNcore) must be given',ch10,& 1232 & 'on the same radial mesh as core density (Ncore) !',ch10,& 1233 & 'Action: check your pseudopotential file.' 1234 MSG_ERROR(msg) 1235 end if 1236 LIBPAW_POINTER_ALLOCATE(tncore,(core_mesh%mesh_size)) 1237 read (funit,*) (tncore(ir),ir=1,core_mesh%mesh_size) 1238 if (maxval(abs(tncore(:)))<tol6) then 1239 pawtab%usetcore=0 1240 else 1241 pawtab%usetcore=1 1242 pawtab%tcoredens(1:pawtab%core_mesh_size,1)=tncore(1:pawtab%core_mesh_size) 1243 end if 1244 write(msg,'(a,i1)') & 1245 & ' Radial grid used for (t)core density is grid ',icoremesh 1246 call wrtout(ab_out,msg,'COLL') 1247 call wrtout(std_out, msg,'COLL') 1248 1249 !--------------------------------- 1250 !Read frozen part of Dij terms (dij0) 1251 LIBPAW_ALLOCATE(pawtab%dij0,(pawtab%lmn2_size)) 1252 read (funit,*) 1253 read (funit,*) (pawtab%dij0(ib),ib=1,pawtab%lmn2_size) 1254 1255 !--------------------------------- 1256 !Read initial guess of rhoij (rhoij0) 1257 LIBPAW_ALLOCATE(pawtab%rhoij0,(pawtab%lmn2_size)) 1258 read (funit,*) 1259 read (funit,*) (pawtab%rhoij0(ib),ib=1,pawtab%lmn2_size) 1260 1261 !--------------------------------- 1262 !Read local pseudopotential=Vh(tn_zc) or Vbare 1263 read (funit,*) 1264 if (pspversion==1) ivlocmesh=3 1265 vlocopt=1 1266 if (pspversion==2) then 1267 read (funit,*) ivlocmesh 1268 else if (pspversion>2) then 1269 ! read (funit,fmt=*,err=30,end=30) ivlocmesh,vlocopt 1270 msg=blank 1271 read (funit,fmt='(a)') msg 1272 read (msg,fmt=*) ivlocmesh 1273 write(numb,'(i1)')ivlocmesh 1274 ii=index(msg,numb) 1275 if(len_trim(trim(msg(ii+1:)))/=0)then 1276 submsg=trim(msg(ii+1:)) 1277 if(len_trim(submsg)/=0)then 1278 do ii=1,len_trim(submsg) 1279 numb=submsg(ii:ii) 1280 if(numb==blank)cycle 1281 jj=index('0123456789',numb) 1282 if(jj<1 .or. jj>10)exit 1283 vlocopt=jj-1 1284 end do 1285 end if 1286 end if 1287 end if 1288 usexcnhat_out=0;if (vlocopt==1) usexcnhat_out=1 1289 call pawrad_copy(radmesh(ivlocmesh),vloc_mesh) 1290 LIBPAW_POINTER_ALLOCATE(vlocr,(vloc_mesh%mesh_size)) 1291 read (funit,*) (vlocr(ir),ir=1,vloc_mesh%mesh_size) 1292 write(msg,'(a,i1)') & 1293 & ' Radial grid used for Vloc is grid ',ivlocmesh 1294 call wrtout(ab_out,msg,'COLL') 1295 call wrtout(std_out, msg,'COLL') 1296 1297 !--------------------------------- 1298 !Eventually read "numeric" shapefunctions (if shape_type=-1) 1299 if (pawtab%shape_type==-1) then 1300 LIBPAW_ALLOCATE(pawtab%shapefunc,(pawtab%mesh_size,pawtab%l_size)) 1301 do il=1,pawtab%l_size 1302 read (funit,*) 1303 if (pspversion==1) iread1=1 1304 if (pspversion>1) read (funit,*) iread1 1305 if (il==1) then 1306 call pawrad_copy(radmesh(iread1),shpf_mesh) 1307 ishpfmesh=iread1 1308 LIBPAW_ALLOCATE(shpf,(shpf_mesh%mesh_size,pawtab%l_size)) 1309 else if (iread1/=ishpfmesh) then 1310 write(msg, '(a,a,a)' )& 1311 & 'All shape functions must be given on the same radial mesh !',ch10,& 1312 & 'Action: check your pseudopotential file.' 1313 MSG_ERROR(msg) 1314 end if 1315 read (funit,*) (shpf(ir,il),ir=1,shpf_mesh%mesh_size) 1316 end do 1317 write(msg,'(a,i1)') & 1318 & ' Radial grid used for shape functions is grid ',iread1 1319 call wrtout(ab_out,msg,'COLL') 1320 call wrtout(std_out, msg,'COLL') 1321 1322 ! Has to spline shape functions if mesh is not the "main" mesh 1323 if (ishpfmesh/=imainmesh) then 1324 msz=shpf_mesh%mesh_size 1325 LIBPAW_ALLOCATE(work1,(msz)) 1326 LIBPAW_ALLOCATE(work2,(msz)) 1327 LIBPAW_ALLOCATE(work3,(msz)) 1328 LIBPAW_ALLOCATE(work4,(pawtab%mesh_size)) 1329 work3(1:pawtab%mesh_size)=shpf_mesh%rad(1:pawtab%mesh_size) 1330 work4(1:pawtab%mesh_size)=pawrad%rad(1:pawtab%mesh_size) 1331 do il=1,pawtab%l_size 1332 call bound_deriv(shpf(1:msz,il),shpf_mesh,msz,yp1,ypn) 1333 call paw_spline(work3,shpf(:,il),msz,yp1,ypn,work1) 1334 call paw_splint(msz,work3,shpf(:,il),work1,pawtab%mesh_size,work4,pawtab%shapefunc(:,il)) 1335 end do 1336 LIBPAW_DEALLOCATE(work1) 1337 LIBPAW_DEALLOCATE(work2) 1338 LIBPAW_DEALLOCATE(work3) 1339 LIBPAW_DEALLOCATE(work4) 1340 else 1341 pawtab%shapefunc(:,:)=shpf(:,:) 1342 end if 1343 LIBPAW_DEALLOCATE(shpf) 1344 call pawrad_free(shpf_mesh) 1345 end if 1346 1347 !--------------------------------- 1348 !Read pseudo valence density (if psp version >=4) 1349 if (pspversion>=4) then 1350 read (funit,*) 1351 read (funit,*) iread1 1352 ivalemesh=iread1 1353 call pawrad_copy(radmesh(iread1),vale_mesh) 1354 LIBPAW_POINTER_ALLOCATE(tnvale,(vale_mesh%mesh_size)) 1355 read (funit,*) (tnvale(ir),ir=1,vale_mesh%mesh_size) 1356 pawtab%has_tvale=1 1357 write(msg,'(a,i1)') & 1358 & ' Radial grid used for pseudo valence density is grid ',ivalemesh 1359 call wrtout(ab_out,msg,'COLL') 1360 call wrtout(std_out, msg,'COLL') 1361 else 1362 pawtab%has_tvale=0 1363 LIBPAW_POINTER_ALLOCATE(tnvale,(0)) 1364 end if 1365 1366 end subroutine pawpsp_read
m_pawpsp/pawpsp_read_corewf [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_read_corewf
FUNCTION
INPUTS
[filename]= (optional) core WF file name
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
optics_paw_core,posdoppler
CHILDREN
SOURCE
1392 subroutine pawpsp_read_corewf(energy_cor,indlmn_core,lcor,lmncmax,ncor,nphicor,radmesh,phi_cor,& 1393 & filename) ! optional argument 1394 1395 1396 !This section has been created automatically by the script Abilint (TD). 1397 !Do not modify the following lines by hand. 1398 #undef ABI_FUNC 1399 #define ABI_FUNC 'pawpsp_read_corewf' 1400 !End of the abilint section 1401 1402 implicit none 1403 1404 !Arguments ------------------------------------ 1405 integer,intent(out) :: lmncmax,nphicor 1406 character(len=*),optional :: filename 1407 !arrays 1408 integer,allocatable,intent(inout) :: indlmn_core(:,:),lcor(:),ncor(:) 1409 real(dp),allocatable,intent(inout) :: phi_cor(:,:),energy_cor(:) 1410 type(pawrad_type),intent(in) :: radmesh 1411 1412 !Local variables------------------------------- 1413 integer :: ib,i1,i2,il,ilm,ilmn,iln,ios,jln,nmesh,npts,unt 1414 real(dp) :: noccor,r1,r2 1415 logical :: ex,oldformat,usexml 1416 character(len=8) :: dum,dum1,dum2,dum3,dum4 1417 character(len=80) :: fline 1418 character(len=500) :: filename_,msg 1419 !arrays 1420 integer,allocatable :: meshtp(:),meshsz(:) 1421 real(dp),allocatable :: rad(:),radstp(:),work(:) 1422 real(dp),allocatable :: logstp(:),phitmp(:) 1423 type(pawrad_type) :: tmpmesh 1424 1425 ! ************************************************************************ 1426 1427 !Check for core WF file existence and XML format 1428 usexml=.false.;oldformat=.false. 1429 if (present(filename)) then 1430 ! Core WF file given as optional argument 1431 filename_=filename;ex=.false. 1432 inquire(file=trim(filename_),iostat=ios,exist=ex) 1433 if (ios/=0) then 1434 write(msg,'(2a)') 'INQUIRE returns an error for file ',trim(filename_) 1435 MSG_ERROR(msg) 1436 end if 1437 if (.not.ex) then 1438 write(msg,'(3a)') 'This file does not exist: ',trim(filename_),'!' 1439 MSG_ERROR(msg) 1440 end if 1441 unt = libpaw_get_free_unit() 1442 open(unit=unt,file=trim(filename_),form='formatted',status='old', action="read") 1443 read(unt,*) fline 1444 close(unt) 1445 usexml=(fline(1:5)=='<?xml') 1446 else 1447 ! Core WF file: new format 1448 filename_='corewf.abinit';ex=.false. 1449 inquire(file=trim(filename_),iostat=ios,exist=ex) 1450 if (ios/=0) then 1451 write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!' 1452 MSG_ERROR(msg) 1453 end if 1454 if (.not.ex) then 1455 ! Core WF file: new format XML 1456 filename_='corewf.abinit.xml';ex=.false. 1457 inquire(file=trim(filename_),iostat=ios,exist=ex) 1458 if (ios/=0) then 1459 write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!' 1460 MSG_ERROR(msg) 1461 end if 1462 usexml=ex 1463 if (.not.ex) then 1464 ! Core WF file: old format 1465 filename_='corewf.dat';ex=.false. 1466 inquire(file=trim(filename_),iostat=ios,exist=ex) 1467 if (ios/=0) then 1468 write(msg,'(3a)') 'INQUIRE returns an error for file ',trim(filename_),'!' 1469 MSG_ERROR(msg) 1470 end if 1471 oldformat=ex 1472 if (.not.ex) then 1473 ! No core WF file found 1474 write(msg, '(3a)' )& 1475 & 'Checks for existence of files corewf.abinit[.xml] or corewf.dat',ch10,& 1476 & 'but INQUIRE finds file does not exist!' 1477 MSG_ERROR(msg) 1478 end if 1479 end if 1480 end if 1481 end if 1482 1483 !Core WF file is in new XML format 1484 if ((.not.oldformat).and.(usexml)) then 1485 call rdpawpsxml_core(energy_cor,trim(filename_),lcor,ncor,nphicor,radmesh,phi_cor) 1486 endif 1487 1488 !Core WF file is in new (proprietary) format 1489 if ((.not.oldformat).and.(.not.usexml)) then 1490 unt = libpaw_get_free_unit() 1491 open(unt,file=trim(filename_),form='formatted',action="read") 1492 read(unt,*) ! skip title 1493 read(unt,*) ! skip method,nspinor,nsppol 1494 read(unt,*) ! skip zatom,zcore,pspdat 1495 read(unt,*) ! skip pspcod,pspxc,lmax 1496 read(unt,*) ! skip pspfmt,creatorID 1497 read(unt,*) nphicor 1498 read(unt,*) ! skip orbitals 1499 read(unt,*) nmesh 1500 LIBPAW_ALLOCATE(meshsz,(nmesh)) 1501 LIBPAW_ALLOCATE(meshtp,(nmesh)) 1502 LIBPAW_ALLOCATE(radstp,(nmesh)) 1503 LIBPAW_ALLOCATE(logstp,(nmesh)) 1504 do iln=1,nmesh 1505 r2=zero;read(unt,'(a80)') fline 1506 read(unit=fline,fmt=*,err=20,end=20) ib,i1,i2,r1,r2 1507 20 continue 1508 if (ib<=nmesh) then 1509 meshtp(ib)=i1;meshsz(ib)=i2 1510 radstp(ib)=r1;logstp(ib)=r2 1511 end if 1512 end do 1513 read(unt,*) ! skip rmax(core) 1514 LIBPAW_ALLOCATE(ncor,(nphicor)) 1515 LIBPAW_ALLOCATE(lcor,(nphicor)) 1516 LIBPAW_ALLOCATE(energy_cor,(nphicor)) 1517 LIBPAW_ALLOCATE(phi_cor,(radmesh%mesh_size,nphicor)) 1518 do iln=1,nphicor 1519 read(unt,*) ! skip comment 1520 read(unt,*) i1 1521 read(unt,*) ncor(iln),lcor(iln) 1522 read(unt,*) energy_cor(iln) 1523 energy_cor(iln)=energy_cor(iln)*half ! For consistency reasons (in the legacy coreWF format, ernegies are in Ry) 1524 LIBPAW_ALLOCATE(phitmp,(meshsz(i1))) 1525 read(unt,*) phitmp 1526 if ((radmesh%mesh_type/=meshtp(i1)) & 1527 & .or.(radmesh%rstep/=radstp(i1)) & 1528 & .or.(radmesh%lstep/=logstp(i1))) then 1529 call pawrad_init(tmpmesh,mesh_size=meshsz(i1),mesh_type=meshtp(i1),rstep=radstp(i1),lstep=logstp(i1)) 1530 npts=radmesh%mesh_size 1531 if (tmpmesh%rmax<radmesh%rmax+tol8) npts=pawrad_ifromr(radmesh,tmpmesh%rmax)-1 1532 LIBPAW_ALLOCATE(work,(meshsz(i1))) 1533 call bound_deriv(phitmp,tmpmesh,meshsz(i1),r1,r2) 1534 call paw_spline(tmpmesh%rad,phitmp,meshsz(i1),r1,r2,work) 1535 call paw_splint(meshsz(i1),tmpmesh%rad,phitmp,work,npts,radmesh%rad(1:npts),phi_cor(1:npts,iln)) 1536 if (npts<radmesh%mesh_size) phi_cor(npts+1:radmesh%mesh_size,iln)=zero 1537 LIBPAW_DEALLOCATE(work) 1538 call pawrad_free(tmpmesh) 1539 else 1540 npts=min(meshsz(i1),radmesh%mesh_size) 1541 phi_cor(1:npts,iln)=phitmp(1:npts) 1542 if (npts<radmesh%mesh_size) phi_cor(npts+1:radmesh%mesh_size,iln)=zero 1543 end if 1544 LIBPAW_DEALLOCATE(phitmp) 1545 end do 1546 LIBPAW_DEALLOCATE(meshsz) 1547 LIBPAW_DEALLOCATE(meshtp) 1548 LIBPAW_DEALLOCATE(radstp) 1549 LIBPAW_DEALLOCATE(logstp) 1550 close(unt) 1551 end if 1552 1553 !Core WF file is in old (proprietary) format 1554 if ((oldformat).and.(.not.usexml)) then 1555 unt = libpaw_get_free_unit() 1556 open(unt,file=trim(filename_),form='formatted',action="read") 1557 do while (dum/='atompaw ') 1558 read(unt,'(a8)') dum 1559 end do 1560 read(unt,'(2i4)') npts,nphicor 1561 LIBPAW_ALLOCATE(ncor,(nphicor)) 1562 LIBPAW_ALLOCATE(lcor,(nphicor)) 1563 LIBPAW_ALLOCATE(energy_cor,(nphicor)) 1564 LIBPAW_ALLOCATE(phi_cor,(npts,nphicor)) 1565 LIBPAW_ALLOCATE(rad,(npts)) 1566 do iln=1,nphicor 1567 read(unt,'(a4,i4,a3,i4,a6,f15.7,a8,f15.7)') & 1568 & dum1,ncor(iln),dum2,lcor(iln),dum3,noccor,dum4,energy_cor(iln) 1569 energy_cor(iln)=energy_cor(iln)*half ! For consistency reasons (in the legacy coreWF format, ernegies are in Ry) 1570 1571 do jln=1,npts 1572 read(unt,*) rad(jln),phi_cor(jln,iln) 1573 end do 1574 read(unt,*) 1575 end do 1576 LIBPAW_DEALLOCATE(rad) 1577 close(unt) 1578 end if 1579 1580 !Set an array 'a la' indlmn 1581 lmncmax=0 1582 do ib=1,nphicor 1583 il=lcor(ib) 1584 lmncmax=lmncmax+2*il+1 1585 end do 1586 LIBPAW_ALLOCATE(indlmn_core,(6,lmncmax)) 1587 indlmn_core=0;ilmn=0;iln=0 1588 do ib=1,nphicor 1589 il=lcor(ib) 1590 iln=iln+1 1591 do ilm=1,2*il+1 1592 indlmn_core(1,ilmn+ilm)=il 1593 indlmn_core(2,ilmn+ilm)=ilm-(il+1) 1594 indlmn_core(3,ilmn+ilm)=1 1595 indlmn_core(4,ilmn+ilm)=il*il+ilm 1596 indlmn_core(5,ilmn+ilm)=iln 1597 indlmn_core(6,ilmn+ilm)=1 1598 end do 1599 ilmn=ilmn+2*il+1 1600 end do 1601 1602 end subroutine pawpsp_read_corewf
m_pawpsp/pawpsp_read_header [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_read_header
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
4091 subroutine pawpsp_read_header(funit,lloc,lmax,mmax,pspcod,pspxc,r2well,zion,znucl) 4092 4093 4094 !This section has been created automatically by the script Abilint (TD). 4095 !Do not modify the following lines by hand. 4096 #undef ABI_FUNC 4097 #define ABI_FUNC 'pawpsp_read_header' 4098 !End of the abilint section 4099 4100 implicit none 4101 !Arguments ------------------------------------ 4102 !scalars 4103 integer,intent(in):: funit 4104 integer,intent(out):: lloc,lmax,mmax,pspcod,pspxc 4105 real(dp),intent(out):: r2well,zion,znucl 4106 !Local variables------------------------------- 4107 integer:: pspdat 4108 character(len=fnlen):: title 4109 character(len=500) :: msg 4110 4111 ! ************************************************************************* 4112 4113 !Read and write some description of file from first line (character data) 4114 read (funit,'(a)') title 4115 write(msg, '(a,a)' ) '- ',trim(title) 4116 call wrtout(ab_out,msg,'COLL') 4117 call wrtout(std_out, msg,'COLL') 4118 4119 !Read and write more data describing psp parameters 4120 read (funit,*) znucl,zion,pspdat 4121 write(msg, '(a,f9.5,f10.5,2x,i8,t47,a)' ) & 4122 & '-',znucl,zion,pspdat,'znucl, zion, pspdat' 4123 call wrtout(ab_out,msg,'COLL') 4124 call wrtout(std_out, msg,'COLL') 4125 4126 read (funit,*) pspcod,pspxc,lmax,lloc,mmax,r2well 4127 if(pspxc<0) then 4128 write(msg, '(i5,i8,2i5,i10,f10.5,t47,a)' ) & 4129 & pspcod,pspxc,lmax,lloc,mmax,r2well,& 4130 & 'pspcod,pspxc,lmax,lloc,mmax,r2well' 4131 else 4132 write(msg, '(4i5,i10,f10.5,t47,a)' ) & 4133 & pspcod,pspxc,lmax,lloc,mmax,r2well,& 4134 & 'pspcod,pspxc,lmax,lloc,mmax,r2well' 4135 end if 4136 call wrtout(ab_out,msg,'COLL') 4137 call wrtout(std_out, msg,'COLL') 4138 4139 end subroutine pawpsp_read_header
m_pawpsp/pawpsp_read_header_2 [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_read_header_2
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
Reads pspversion, basis_size and lmn_size
PARENTS
m_pawpsp
CHILDREN
SOURCE
4168 subroutine pawpsp_read_header_2(funit,pspversion,basis_size,lmn_size) 4169 4170 4171 !This section has been created automatically by the script Abilint (TD). 4172 !Do not modify the following lines by hand. 4173 #undef ABI_FUNC 4174 #define ABI_FUNC 'pawpsp_read_header_2' 4175 !End of the abilint section 4176 4177 implicit none 4178 4179 !Arguments ------------------------------------ 4180 !scalars 4181 integer,intent(in):: funit 4182 integer,intent(out) :: pspversion,basis_size,lmn_size 4183 4184 !Local variables------------------------------- 4185 integer :: creatorid 4186 character(len=80) :: pspline 4187 character(len=500) :: msg 4188 4189 ! ************************************************************************* 4190 4191 !Read psp version in line 4 of the header 4192 pspversion=1 4193 read (funit,'(a80)') pspline;pspline=adjustl(pspline) 4194 if (pspline(1:3)=="paw".or.pspline(1:3)=="PAW") & 4195 & read(unit=pspline(4:80),fmt=*) pspversion 4196 if (pspversion<1.or.pspversion>5) then 4197 write(msg, '(a,i2,a,a,a)' )& 4198 & 'This version of PAW psp file (',pspversion,') is not compatible with',ch10,& 4199 & 'current version of Abinit.' 4200 MSG_ERROR(msg) 4201 end if 4202 4203 if (pspversion==1) then 4204 read (unit=pspline,fmt=*) basis_size,lmn_size 4205 else 4206 ! Here psp file for Abinit 4.3+ 4207 read (unit=pspline(5:80),fmt=*) creatorid 4208 read (funit,*) basis_size,lmn_size 4209 end if 4210 4211 end subroutine pawpsp_read_header_2
m_pawpsp/pawpsp_read_header_xml [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_read_header_xml
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
This is done instead of: call pawpsxml2ab( psxml, pspheads,1) since pspheads does not exist in PAW library. should we include it to avoid the following code replica? check pspheads commented out in pawpsp_17in, and routine pawpsp_read_xml_2
PARENTS
m_pawpsp,pawpsxml2ab,pspatm
CHILDREN
SOURCE
4364 subroutine pawpsp_read_header_xml(lloc,lmax,pspcod,pspxc,& 4365 & psxml,r2well,zion,znucl) 4366 4367 4368 !This section has been created automatically by the script Abilint (TD). 4369 !Do not modify the following lines by hand. 4370 #undef ABI_FUNC 4371 #define ABI_FUNC 'pawpsp_read_header_xml' 4372 !End of the abilint section 4373 4374 implicit none 4375 !Arguments ------------------------------------ 4376 !scalars 4377 type(paw_setup_t),intent(in) :: psxml 4378 integer,intent(out):: lloc,lmax,pspcod,pspxc 4379 real(dp),intent(out):: r2well,zion,znucl 4380 !Local variables------------------------------- 4381 integer :: il 4382 #if defined LIBPAW_HAVE_LIBXC 4383 integer :: ii 4384 #endif 4385 character(len=100) :: xclibxc 4386 character(len=500) :: msg 4387 !arrays 4388 4389 ! ************************************************************************* 4390 4391 lloc = 0 4392 r2well = 0 4393 pspcod=17 4394 znucl=psxml%atom%znucl 4395 zion =psxml%atom%zval 4396 4397 !lmax: 4398 lmax = 0 4399 do il=1,psxml%valence_states%nval 4400 if(psxml%valence_states%state(il)%ll>lmax) lmax=psxml%valence_states%state(il)%ll 4401 end do 4402 !pspxc 4403 select case(trim(psxml%xc_functional%name)) 4404 case('PZ') 4405 pspxc = 2 4406 #if defined LIBPAW_HAVE_LIBXC 4407 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4408 & +libxc_functionals_getid('XC_LDA_C_PZ')) 4409 #endif 4410 case('W') 4411 pspxc = 4 4412 #if defined LIBPAW_HAVE_LIBXC 4413 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4414 & +libxc_functionals_getid('XC_LDA_C_WIGNER')) 4415 #endif 4416 case('HL') 4417 pspxc = 5 4418 #if defined LIBPAW_HAVE_LIBXC 4419 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4420 & +libxc_functionals_getid('XC_LDA_C_HL')) 4421 #endif 4422 case('GL') 4423 #if defined LIBPAW_HAVE_LIBXC 4424 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4425 & +libxc_functionals_getid('XC_LDA_C_GL')) 4426 #else 4427 write(msg, '(7a)' )& 4428 & 'The exchange and correlation functional by Gunnarson-Lundqvist', ch10,& 4429 & 'is not implemented in Abinit.',ch10,& 4430 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4431 & ' generation or compile ABINIT with the libXC library.' 4432 MSG_ERROR(msg) 4433 #endif 4434 case('VWN') 4435 #if defined LIBPAW_HAVE_LIBXC 4436 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4437 & +libxc_functionals_getid('XC_LDA_C_VWN')) 4438 #else 4439 write(msg, '(7a)' )& 4440 & 'The exchange and correlation functional by Vosko,Wilk and Nusair', ch10,& 4441 & 'is not implemented in Abinit.',ch10,& 4442 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4443 & ' generation or compile ABINIT with the libXC library.' 4444 MSG_ERROR(msg) 4445 #endif 4446 case('PW') 4447 pspxc = 7 4448 #if defined LIBPAW_HAVE_LIBXC 4449 pspxc = -(libxc_functionals_getid('XC_LDA_X')*1000 & 4450 & +libxc_functionals_getid('XC_LDA_C_PW')) 4451 #endif 4452 case('PBE') 4453 pspxc = 11 4454 #if defined LIBPAW_HAVE_LIBXC 4455 pspxc = -(libxc_functionals_getid('XC_GGA_X_PBE')*1000 & 4456 & +libxc_functionals_getid('XC_GGA_C_PBE')) 4457 #endif 4458 case('revPBE') 4459 pspxc = 14 4460 #if defined LIBPAW_HAVE_LIBXC 4461 pspxc = -(libxc_functionals_getid('XC_GGA_X_PBE_R')*1000 & 4462 & +libxc_functionals_getid('XC_GGA_C_PBE')) 4463 #endif 4464 case('RPBE') 4465 pspxc = 15 4466 #if defined LIBPAW_HAVE_LIBXC 4467 pspxc = -(libxc_functionals_getid('XC_GGA_X_RPBE')*1000 & 4468 & +libxc_functionals_getid('XC_GGA_C_PBE')) 4469 #endif 4470 case('PW91') 4471 #if defined LIBPAW_HAVE_LIBXC 4472 pspxc = -(libxc_functionals_getid('XC_GGA_X_PW91')*1000 & 4473 & +libxc_functionals_getid('XC_GGA_C_PW91')) 4474 #else 4475 write(msg, '(7a)' )& 4476 & 'The exchange and correlation functional by Perdew and Wang 91', ch10,& 4477 & 'is not implemented in Abinit.',ch10,& 4478 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4479 & ' generation or compile ABINIT with the libXC library.' 4480 MSG_ERROR(msg) 4481 #endif 4482 case('BLYP') 4483 #if defined LIBPAW_HAVE_LIBXC 4484 pspxc = -(libxc_functionals_getid('XC_GGA_X_B88')*1000 & 4485 & +libxc_functionals_getid('XC_GGA_C_LYP')) 4486 #else 4487 write(msg, '(7a)' )& 4488 & 'The exchange and correlation functional BLYP', ch10,& 4489 & 'is not implemented in Abinit.',ch10,& 4490 & 'Action : choose another XC functional in the pseudopotential',ch10, & 4491 & ' generation or compile ABINIT with the libXC library.' 4492 MSG_ERROR(msg) 4493 #endif 4494 case DEFAULT 4495 xclibxc=trim(psxml%xc_functional%name) 4496 if (xclibxc(1:3)=='XC_' .or.xclibxc(1:3)=='xc_' .or. & 4497 & xclibxc(1:5)=='LDA_X'.or.xclibxc(1:5)=='LDA_C'.or. & 4498 & xclibxc(1:5)=='lda_x'.or.xclibxc(1:5)=='lda_c'.or. & 4499 & xclibxc(1:5)=='GGA_X'.or.xclibxc(1:5)=='GGA_C'.or. & 4500 & xclibxc(1:5)=='gga_x'.or.xclibxc(1:5)=='gga_c') then 4501 #if defined LIBPAW_HAVE_LIBXC 4502 ii=index(xclibxc,'+') 4503 if (ii>0) then 4504 pspxc=-(libxc_functionals_getid(xclibxc(1:ii-1))*1000 & 4505 & +libxc_functionals_getid(xclibxc(ii+1:))) 4506 else 4507 pspxc=-libxc_functionals_getid(xclibxc) 4508 end if 4509 #else 4510 msg='Cannot use LibXC functional because ABINIT is not compiled with LibXC !' 4511 MSG_ERROR(msg) 4512 #endif 4513 ! To be eliminated later (temporary) 4514 else if(trim(psxml%xc_functional%functionaltype)=='LIBXC')then 4515 #if defined LIBPAW_HAVE_LIBXC 4516 xclibxc=trim(psxml%xc_functional%name) 4517 read(unit=xclibxc,fmt=*) pspxc 4518 pspxc=-pspxc 4519 #else 4520 msg='Cannot use LibXC functional because ABINIT is not compiled with LibXC !' 4521 MSG_ERROR(msg) 4522 #endif 4523 else 4524 write(msg, '(3a)') 'Unknown XC functional in psp file: ',trim(xclibxc),' !' 4525 MSG_ERROR(msg) 4526 end if 4527 end select 4528 4529 end subroutine pawpsp_read_header_xml
m_pawpsp/pawpsp_read_pawheader [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_read_pawheader
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp,pawpsxml2ab,pspatm
CHILDREN
SOURCE
4555 subroutine pawpsp_read_pawheader(basis_size,lmax,lmn_size,& 4556 & l_size,mesh_size,pspversion,psxml,rpaw,rshp,shape_type) 4557 4558 4559 !This section has been created automatically by the script Abilint (TD). 4560 !Do not modify the following lines by hand. 4561 #undef ABI_FUNC 4562 #define ABI_FUNC 'pawpsp_read_pawheader' 4563 !End of the abilint section 4564 4565 implicit none 4566 !Arguments ------------------------------------ 4567 !scalars 4568 integer,intent(in):: lmax 4569 integer,intent(out):: basis_size,mesh_size,lmn_size,l_size 4570 integer,intent(out):: pspversion,shape_type 4571 real(dp),intent(out)::rpaw,rshp 4572 type(paw_setup_t),intent(in) :: psxml 4573 !Local variables------------------------------- 4574 integer::il 4575 4576 ! ************************************************************************* 4577 4578 !All of this was moved from pawpsxml2ab, 4579 !basis_size 4580 basis_size=psxml%valence_states%nval 4581 !mesh_size 4582 do il=1,psxml%ngrid 4583 if(psxml%radial_grid(il)%id==psxml%idgrid) & 4584 & mesh_size=psxml%radial_grid(il)%iend-psxml%radial_grid(il)%istart+1 4585 end do 4586 !lmn_size: 4587 lmn_size=0 4588 do il=1,psxml%valence_states%nval 4589 lmn_size=lmn_size+2*psxml%valence_states%state(il)%ll+1 4590 end do 4591 !lsize 4592 l_size=2*lmax+1 4593 !pspversion 4594 pspversion=10 4595 !rpaw: 4596 rpaw=0.d0 4597 if (psxml%rpaw<0.d0) then 4598 do il=1,psxml%valence_states%nval 4599 if(psxml%valence_states%state(il)%rc>rpaw) rpaw=psxml%valence_states%state(il)%rc 4600 end do 4601 else 4602 rpaw=psxml%rpaw 4603 end if 4604 !shape_type, rshp: 4605 select case(trim(psxml%shape_function%gtype)) 4606 case('gauss') 4607 shape_type=1 4608 rshp=rpaw 4609 case('bessel') 4610 shape_type=3 4611 rshp=psxml%shape_function%rc 4612 case('sinc') 4613 shape_type=2 4614 rshp=psxml%shape_function%rc 4615 case('exp') 4616 shape_type=1 4617 rshp=rpaw 4618 case('num') 4619 shape_type=-1 4620 rshp=rpaw 4621 end select 4622 4623 end subroutine pawpsp_read_pawheader
m_pawpsp/pawpsp_rw_atompaw [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_rw_atompaw
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
1629 subroutine pawpsp_rw_atompaw(basis_size,filpsp,wvl) 1630 1631 1632 !This section has been created automatically by the script Abilint (TD). 1633 !Do not modify the following lines by hand. 1634 #undef ABI_FUNC 1635 #define ABI_FUNC 'pawpsp_rw_atompaw' 1636 !End of the abilint section 1637 1638 implicit none 1639 1640 !Arguments ------------------------------------ 1641 integer,intent(in):: basis_size 1642 type(wvlpaw_type),intent(in)::wvl 1643 character(len=fnlen),intent(in)::filpsp 1644 !arrays 1645 character(strlen) :: pspline 1646 character(len=fnlen)::fname 1647 1648 !Local variables------------------------------- 1649 integer :: ib,ii,ios,jj,step,iunt,ount 1650 !arrays 1651 1652 ! ************************************************************************* 1653 iunt = libpaw_get_free_unit() 1654 ount = libpaw_get_free_unit() 1655 1656 step=0 1657 ! Open psp file for reading 1658 open(unit=iunt,file=trim(filpsp),form='formatted',status='old',action="read") 1659 ! Open the file for writing 1660 write(fname,'(2a)') libpaw_basename(trim(filpsp)),".wvl" 1661 open(unit=ount,file=fname,form='formatted',status='unknown',action="write") 1662 1663 read_loop: do 1664 if(step==0) then 1665 read(iunt,'(a)',IOSTAT=ios) pspline 1666 if ( ios /= 0 ) exit read_loop 1667 if(index(trim(pspline),'CORE_DENSITY')/=0 .and. & 1668 & index(trim(pspline),'PSEUDO_CORE_DENSITY')==0 .and. & 1669 & index(trim(pspline),'TCORE_DENSITY')==0 ) then 1670 step=1 1671 else 1672 write(ount,'(a)') trim(pspline) 1673 end if 1674 elseif(step==1) then 1675 !Write Gaussian projectors: 1676 jj=0 1677 do ib=1,basis_size 1678 write(ount,'(a,i1,a)') "===== GAUSSIAN_TPROJECTOR ",ib,& 1679 & " ===== " 1680 write(ount,'(i5,1x,i5,1x,a)')wvl%pngau(ib),wvl%ptotgau, ":ngauss, total ngauss" 1681 write(ount,'(3(1x,es23.16))')(wvl%parg(:,ii),& 1682 & ii=jj+1,jj+wvl%pngau(ib)) 1683 write(ount,'(3(1x,es23.16))')(wvl%pfac(:,ii),& 1684 & ii=jj+1,jj+wvl%pngau(ib)) 1685 jj=jj+wvl%pngau(ib) 1686 end do 1687 write(ount,'(a)') trim(pspline) 1688 step=0 1689 end if 1690 end do read_loop 1691 1692 close(iunt) 1693 close(ount) 1694 1695 end subroutine pawpsp_rw_atompaw
m_pawpsp/pawpsp_vhar2rho [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_vhar2rho
FUNCTION
gets rho(r) from v(r), solving the Poisson equation \lap v(r) = 4 \pi rho(r)
INPUTS
radmesh = radial grid (datastructure) vv(:)= potential
OUTPUT
rho(:)= density
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
2716 subroutine pawpsp_vhar2rho(radmesh,rho,vv) 2717 2718 2719 !This section has been created automatically by the script Abilint (TD). 2720 !Do not modify the following lines by hand. 2721 #undef ABI_FUNC 2722 #define ABI_FUNC 'pawpsp_vhar2rho' 2723 !End of the abilint section 2724 2725 implicit none 2726 2727 !Arguments ------------------------------------ 2728 type(pawrad_type),intent(in) :: radmesh 2729 real(dp), intent(in) :: vv(:) 2730 real(dp), intent(out):: rho(:) 2731 2732 !Local variables------------------------------- 2733 integer :: nr 2734 real(dp) :: dfdr(radmesh%mesh_size),d2fdr(radmesh%mesh_size) 2735 2736 ! ************************************************************************* 2737 2738 nr=size(vv) 2739 if (nr/=size(rho)) then 2740 MSG_BUG('wrong sizes!') 2741 end if 2742 2743 !Laplacian = 2744 !\frac{\partial^2}{\partial r^2} + 2/r \frac{\partial}{\partial r} 2745 2746 !Calculate derivatives 2747 call nderiv_gen(dfdr(1:nr),vv,radmesh,der2=d2fdr(1:nr)) 2748 2749 rho(2:nr)=d2fdr(2:nr) + 2._dp*dfdr(2:nr)/radmesh%rad(2:nr) 2750 call pawrad_deducer0(rho,nr,radmesh) 2751 2752 rho(1:nr)=-rho(1:nr)/(4._dp*pi) 2753 2754 end subroutine pawpsp_vhar2rho
m_pawpsp/pawpsp_wvl [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_wvl
FUNCTION
WVL+PAW related operations
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp,pspatm
CHILDREN
SOURCE
4239 subroutine pawpsp_wvl(filpsp,pawrad, pawtab,usewvl, wvl_ngauss, comm_mpi) 4240 4241 4242 !This section has been created automatically by the script Abilint (TD). 4243 !Do not modify the following lines by hand. 4244 #undef ABI_FUNC 4245 #define ABI_FUNC 'pawpsp_wvl' 4246 !End of the abilint section 4247 4248 implicit none 4249 4250 !Arguments------------------------------------ 4251 !scalars 4252 integer, optional,intent(in):: comm_mpi 4253 integer, intent(in):: usewvl, wvl_ngauss(2) 4254 character(len=fnlen),intent(in)::filpsp 4255 type(pawrad_type),intent(in) :: pawrad 4256 type(pawtab_type),intent(inout):: pawtab 4257 !arrays 4258 4259 !Local variables------------------------------- 4260 !scalars 4261 integer:: ii, me, mparam, nterm_bounds(2) 4262 type(pawrad_type)::tproj_mesh 4263 character(len=500) :: msg 4264 !arrays 4265 integer,allocatable:: ngauss_param(:) 4266 real(dp),allocatable:: gauss_param(:,:) 4267 4268 ! ************************************************************************* 4269 4270 me=0; if (present(comm_mpi))me=xmpi_comm_rank(comm_mpi) 4271 4272 !If usewvl flag is on, we must have the pawtab%wvl pointer allocated 4273 if (usewvl==1.and.pawtab%has_wvl==0) then 4274 call wvlpaw_allocate(pawtab%wvl) 4275 pawtab%has_wvl=1 4276 end if 4277 4278 !Fit projectors to a sum of Gaussians: 4279 if (usewvl ==1 .and. pawtab%wvl%ptotgau==0 ) then 4280 4281 if (pawtab%has_tproj==0) then 4282 msg='pawtab%tproj must be allocated' 4283 MSG_BUG(msg) 4284 end if 4285 4286 ! 1) fit projectors to gaussians 4287 write(msg,'(a,a)')ch10,'Fitting tproj to Gaussians' 4288 call wrtout(std_out,msg,'COLL') 4289 4290 ! See fit_gen (option==4): 4291 do ii=1,2 4292 nterm_bounds(ii)=ceiling(wvl_ngauss(ii)/2.0) 4293 end do 4294 mparam=nterm_bounds(2)*4 4295 LIBPAW_ALLOCATE(gauss_param,(mparam,pawtab%basis_size)) 4296 LIBPAW_ALLOCATE(ngauss_param,(pawtab%basis_size)) 4297 ! compute tproj_mesh 4298 call pawrad_init(tproj_mesh,mesh_size=size(pawtab%tproj,1),& 4299 & mesh_type=pawrad%mesh_type,rstep=pawrad%rstep, lstep=pawrad%lstep) 4300 4301 if(present(comm_mpi)) then 4302 call gaussfit_projector(pawtab%basis_size,mparam,& 4303 & ngauss_param,nterm_bounds,pawtab%orbitals,& 4304 & gauss_param,tproj_mesh,& 4305 & pawtab%rpaw,pawtab%tproj,comm_mpi) 4306 else 4307 call gaussfit_projector(pawtab%basis_size,mparam,& 4308 & ngauss_param,nterm_bounds,pawtab%orbitals,& 4309 & gauss_param,tproj_mesh,& 4310 & pawtab%rpaw,pawtab%tproj) 4311 end if 4312 ! tproj is now as a sum of sin+cos functions, 4313 ! convert it to a sum of complex gaussians and fill %wvl object: 4314 call pawpsp_wvl_sin2gauss(pawtab%basis_size,mparam,& 4315 & ngauss_param,gauss_param,pawtab%wvl) 4316 LIBPAW_DEALLOCATE(gauss_param) 4317 LIBPAW_DEALLOCATE(ngauss_param) 4318 4319 if(me==0) then 4320 call pawpsp_rw_atompaw(pawtab%basis_size,filpsp,pawtab%wvl) 4321 end if 4322 4323 pawtab%has_wvl=2 4324 4325 end if 4326 4327 !Projectors in real space are no more needed 4328 call pawrad_free(tproj_mesh) 4329 if(allocated(pawtab%tproj)) then 4330 LIBPAW_DEALLOCATE(pawtab%tproj) 4331 pawtab%has_tproj=0 4332 end if 4333 4334 end subroutine pawpsp_wvl
m_pawpsp/pawpsp_wvl_calc [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_wvl_calc
FUNCTION
Performs tests and compute data related to pspcod=7 or 17 ("PAW pseudopotentials")
INPUTS
tnvale(vale_mesh%mesh_size)= pseudo valence density (+ nhat in output) usewvl= flag for wavelets method vale_mesh<type(pawrad_type)>= radial mesh for the valence density vloc_mesh<type(pawrad_type)>= radial mesh for the local potential vlocr(vloc_mesh%mesh_size)= local potential according to vlocopt.
OUTPUT
Sets pawtab%rholoc
SIDE EFFECTS
pawtab <type(pawtab_type)>= objects are modified
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
2788 subroutine pawpsp_wvl_calc(pawtab,tnvale,usewvl,vale_mesh,vloc_mesh,vlocr) 2789 2790 2791 !This section has been created automatically by the script Abilint (TD). 2792 !Do not modify the following lines by hand. 2793 #undef ABI_FUNC 2794 #define ABI_FUNC 'pawpsp_wvl_calc' 2795 !End of the abilint section 2796 2797 implicit none 2798 2799 !Arguments ------------------------------------ 2800 !scalars 2801 integer,intent(in)::usewvl 2802 type(pawrad_type),intent(in) :: vale_mesh 2803 type(pawtab_type),intent(inout) :: pawtab 2804 type(pawrad_type),intent(in) ::vloc_mesh 2805 2806 !arrays 2807 real(dp),intent(in) :: tnvale(vale_mesh%mesh_size*pawtab%has_tvale) 2808 real(dp),intent(in) :: vlocr(vloc_mesh%mesh_size) 2809 2810 2811 !Local variables ------------------------------ 2812 !scalars 2813 integer :: msz 2814 character(len=500) :: msg 2815 !arrays 2816 2817 ! ************************************************************************* 2818 2819 !If usewvl flag is on, we must have the pawtab%wvl pointer allocated 2820 if (pawtab%has_wvl==0) then 2821 msg='pawtab%has_wvl flag should be on o entry' 2822 MSG_BUG(msg) 2823 end if 2824 call wvlpaw_allocate(pawtab%wvl) 2825 2826 !========================================================== 2827 !Change mesh_size of tvalespl 2828 !Compute second derivative from tNvale(r) 2829 2830 if (pawtab%has_tvale/=0) then 2831 if(usewvl==1) then 2832 if(allocated(pawtab%tvalespl)) then 2833 LIBPAW_DEALLOCATE(pawtab%tvalespl) 2834 end if 2835 LIBPAW_ALLOCATE(pawtab%tvalespl,(vale_mesh%mesh_size,2)) 2836 pawtab%tnvale_mesh_size=vale_mesh%mesh_size 2837 pawtab%tvalespl(:,1)=tnvale 2838 ! Compute second derivative of tvalespl(r) 2839 call paw_spline(vale_mesh%rad,pawtab%tvalespl(:,1),vale_mesh%mesh_size,zero,zero,pawtab%tvalespl(:,2)) 2840 end if 2841 else 2842 pawtab%dnvdq0=zero 2843 pawtab%tnvale_mesh_size=0 2844 end if 2845 2846 !========================================================== 2847 !Save rholoc: 2848 !Get local density from local potential 2849 !use the poisson eq. 2850 msz=vloc_mesh%mesh_size 2851 call wvlpaw_rholoc_free(pawtab%wvl%rholoc) 2852 LIBPAW_ALLOCATE(pawtab%wvl%rholoc%d,(msz,4)) 2853 LIBPAW_ALLOCATE(pawtab%wvl%rholoc%rad,(msz)) 2854 pawtab%wvl%rholoc%msz=msz 2855 pawtab%wvl%rholoc%rad(1:msz)=vloc_mesh%rad(1:msz) 2856 2857 !get rho from v: 2858 call pawpsp_vhar2rho(vloc_mesh,pawtab%wvl%rholoc%d(:,1),vlocr) 2859 ! 2860 !get second derivative, and store it 2861 call paw_spline(pawtab%wvl%rholoc%rad,pawtab%wvl%rholoc%d(:,1),pawtab%wvl%rholoc%msz,& 2862 & zero,zero,pawtab%wvl%rholoc%d(:,2)) 2863 2864 !save also vlocr: 2865 pawtab%wvl%rholoc%d(:,3)=vlocr 2866 2867 !get second derivative, and store it 2868 call paw_spline(pawtab%wvl%rholoc%rad,vlocr,pawtab%wvl%rholoc%msz,& 2869 & zero,zero,pawtab%wvl%rholoc%d(:,4)) 2870 2871 !Test 2872 !do ii=1,pawtab%wvl%rholoc%msz 2873 !write(503,'(3(f16.10,x))')pawtab%wvl%rholoc%rad(ii),pawtab%wvl%rholoc%d(ii,1),pawtab%wvl%rholoc%d(ii,3) 2874 !end do 2875 ! 2876 !Do splint 2877 ! 2878 !nmesh=4000 2879 !rread1= (9.9979999d0/real(nmesh-1,dp)) ! 0.0025001d0 !step 2880 !allocate(raux1(nmesh),raux2(nmesh)) 2881 !do ii=1,nmesh 2882 !raux1(ii)=rread1*real(ii-1,dp) !mesh 2883 !end do 2884 !call splint(pawtab%wvl%rholoc%msz,pawtab%wvl%rholoc%rad,pawtab%wvl%rholoc%d(:,1),pawtab%wvl%rholoc%d(:,2),& 2885 !& nmesh,raux1,raux2,ierr) 2886 !do ii=1,nmesh 2887 !write(401,'(10(f20.7,x))')raux1(ii),raux2(ii),raux2(ii)*raux1(ii)**2 2888 !end do 2889 !deallocate(raux1,raux2) 2890 2891 end subroutine pawpsp_wvl_calc
m_pawpsp/pawpsp_wvl_sin2gauss [ Functions ]
[ Top ] [ m_pawpsp ] [ Functions ]
NAME
pawpsp_wvl_sin2gauss
FUNCTION
Converts a f(x)=sum_i^N_i a_i sin(b_i x)+ c_i cos( d_i x) to
f(x)=sum_j e_j exp(f_j x), where e and f are complex numbers.
INPUTS
basis_size = size of the lmn basis mparam = number of terms in the summatory (N_i, see the expression above) nparam = Array containing the parameters (a_i, b_i,c_i,d_i) wvl = wavelets data type
OUTPUT
SIDE EFFECTS
On output wvl%pfac and wvl%parg are filled with complex parameters (e_i, f_i)
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
3930 subroutine pawpsp_wvl_sin2gauss(basis_size,mparam,nparam,& 3931 & param,wvl) 3932 3933 3934 !This section has been created automatically by the script Abilint (TD). 3935 !Do not modify the following lines by hand. 3936 #undef ABI_FUNC 3937 #define ABI_FUNC 'pawpsp_wvl_sin2gauss' 3938 !End of the abilint section 3939 3940 implicit none 3941 ! 3942 !Arguments ------------------------------------ 3943 integer,intent(in) :: mparam,basis_size 3944 integer,intent(in) :: nparam(basis_size) 3945 real(dp),intent(in) :: param(mparam,basis_size) 3946 type(wvlpaw_type),intent(inout):: wvl 3947 3948 !Local variables ------------------------------ 3949 integer :: i,ii,ib,ngauss,nterm 3950 real(dp) :: sep 3951 real(dp) :: a1(mparam),a2(mparam),a3(mparam),a4(mparam),a5(mparam) 3952 real(dp) :: b1r(mparam),b2r(mparam),b1i(mparam),b2i(mparam) 3953 character(len=500) :: message 3954 ! 3955 !extra variables, use to debug 3956 ! 3957 !integer::igau,nr,unitp 3958 !real(dp)::step,rmax 3959 !real(dp),allocatable::r(:), y(:) 3960 !complex::fac,arg 3961 !complex(dp),allocatable::f(:) 3962 !************************************************************************ 3963 3964 ! Convert from \sum(sin+cos) expressions to sums of complex gaussians 3965 ! (only works for option=4, see fit_gen) 3966 3967 ! get number of coefficients: 3968 ii=0 3969 do ib=1,basis_size 3970 nterm=nparam(ib)/4 !option=4, there are 4 parameters for each term 3971 ii=ii+nterm*2 !two gaussians for each term 3972 end do 3973 ! 3974 ! Allocate objects 3975 ! 3976 ngauss=ii 3977 wvl%ptotgau=ngauss !total number of complex gaussians 3978 LIBPAW_ALLOCATE(wvl%pfac,(2,ngauss)) 3979 LIBPAW_ALLOCATE(wvl%parg,(2,ngauss)) 3980 LIBPAW_ALLOCATE(wvl%pngau,(basis_size)) 3981 wvl%pngau(1:basis_size)=nparam(1:basis_size)/2 !option=4 3982 ! 3983 ii=0 3984 do ib=1,basis_size 3985 ! 3986 ! Get parameters in sin+cos expansion: 3987 ! Option4: \sum a1 exp(-a2 x^2) ( a3 sin(k x^2) + a4 cos(k x^2)) 3988 ! 3989 nterm=nparam(ib)/4 !option=4 3990 ! 3991 a1(1:nterm)=param(1:nterm,ib) 3992 a2(1:nterm)=param(nterm+1:nterm*2,ib) 3993 a3(1:nterm)=param(nterm*2+1:nterm*3,ib) 3994 a4(1:nterm)=param(nterm*3+1:nterm*4,ib) 3995 sep=1.1d0 3996 do i=1,nterm 3997 a5(i)=sep**(i) 3998 end do 3999 4000 ! First check that "a2" is a positive number (it is multiplied by -1, so 4001 ! that gaussians decay to zero: 4002 if( any(a2(1:nterm) < tol12) ) then 4003 message = 'Real part of Gaussians should be a negative number (they should go to zero at infty)' 4004 MSG_ERROR(message) 4005 end if 4006 4007 ! 4008 ! Now translate them to a sum of complex gaussians: 4009 ! pngau(ib)=nterm*2 4010 ! Two gaussians by term: 4011 ! 4012 ! First gaussian 4013 b1r(1:nterm)= a1(1:nterm)*a4(1:nterm)/2.d0 !coefficient, real 4014 b1i(1:nterm)=-a1(1:nterm)*a3(1:nterm)/2.d0 !coefficient, imag 4015 b2r(1:nterm)=-a2(1:nterm) !exponential, real 4016 b2i(1:nterm)= a5(1:nterm) !exponential, imag 4017 ! 4018 wvl%pfac(1,ii+1:ii+nterm)=b1r(1:nterm) 4019 wvl%pfac(2,ii+1:ii+nterm)=b1i(1:nterm) 4020 wvl%parg(1,ii+1:ii+nterm)=b2r(1:nterm) 4021 wvl%parg(2,ii+1:ii+nterm)=b2i(1:nterm) 4022 ! Second gaussian 4023 wvl%pfac(1,ii+nterm+1:ii+nterm*2)= b1r(1:nterm) 4024 wvl%pfac(2,ii+nterm+1:ii+nterm*2)=-b1i(1:nterm) 4025 wvl%parg(1,ii+nterm+1:ii+nterm*2)= b2r(1:nterm) 4026 wvl%parg(2,ii+nterm+1:ii+nterm*2)=-b2i(1:nterm) 4027 ! 4028 ii=ii+nterm*2 4029 end do 4030 4031 ! begin debug 4032 ! write(*,*)'pawpsp_wvl_sin2gauss, comment me' 4033 ! nr=3000 4034 ! rmax=10.d0 4035 ! LIBPAW_ALLOCATE(r,(nr)) 4036 ! LIBPAW_ALLOCATE(f,(nr)) 4037 ! LIBPAW_ALLOCATE(y,(nr)) 4038 ! step=rmax/real(nr-1,dp) 4039 ! do ir=1,nr 4040 ! r(ir)=real(ir-1,dp)*step 4041 ! end do 4042 ! ! 4043 ! ii=0 4044 ! do ib=1,basis_size 4045 ! unitp=500+ib 4046 ! f(:)=czero 4047 ! ! 4048 ! do igau=1,wvl%pngau(ib) 4049 ! ii=ii+1 4050 ! arg=cmplx(wvl%parg(1,ii),wvl%parg(2,ii)) 4051 ! fac=cmplx(wvl%pfac(1,ii),wvl%pfac(2,ii)) 4052 ! f(:)=f(:)+fac*exp(arg*r(:)**2) 4053 ! end do 4054 ! do ir=1,nr 4055 ! write(unitp,'(3f16.7)')r(ir),real(f(ir))!,y(ir) 4056 ! end do 4057 ! end do 4058 ! LIBPAW_DEALLOCATE(r) 4059 ! LIBPAW_DEALLOCATE(f) 4060 ! LIBPAW_DEALLOCATE(y) 4061 ! end debug 4062 4063 end subroutine pawpsp_wvl_sin2gauss
pawpsp_main/pawpsp_check_xml_upf [ Functions ]
[ Top ] [ pawpsp_main ] [ Functions ]
NAME
pawpsp_main_checks
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
4965 subroutine pawpsp_check_xml_upf(filpsp) 4966 4967 4968 !This section has been created automatically by the script Abilint (TD). 4969 !Do not modify the following lines by hand. 4970 #undef ABI_FUNC 4971 #define ABI_FUNC 'pawpsp_check_xml_upf' 4972 !End of the abilint section 4973 4974 implicit none 4975 4976 !Arguments ------------------------------------ 4977 !scalars 4978 character(len=fnlen),intent(in):: filpsp ! name of the psp file 4979 4980 !Local variables------------------------------- 4981 integer :: unt 4982 character(len=70):: testxml 4983 4984 ! ************************************************************************* 4985 4986 ! Check if the file pseudopotential file is written in XML 4987 usexml = 0 4988 unt = libpaw_get_free_unit() 4989 open (unit=unt,file=filpsp,form='formatted',status='old',action="read") 4990 rewind (unit=unt) 4991 read(unt,*) testxml 4992 if(testxml(1:5)=='<?xml')then 4993 usexml = 1 4994 read(unt,*) testxml 4995 if(testxml(1:4)/='<paw')then 4996 msg='Reading a NC pseudopotential for a PAW calculation?' 4997 MSG_BUG(msg) 4998 end if 4999 else 5000 usexml = 0 5001 end if 5002 close (unit=unt) 5003 5004 ! Check if pseudopotential file is a Q-espresso UPF file 5005 unt = libpaw_get_free_unit() 5006 open (unit=unt,file=filpsp,form='formatted',status='old',action="read") 5007 rewind (unit=unt) 5008 read(unt,*) testxml ! just a string, no relation to xml. 5009 if(testxml(1:9)=='<PP_INFO>')then 5010 msg='UPF format not allowed with PAW (USPP part not read yet)!' 5011 MSG_ERROR(msg) 5012 end if 5013 close (unit=unt) 5014 5015 end subroutine pawpsp_check_xml_upf
pawpsp_main/pawpsp_consistency [ Functions ]
[ Top ] [ pawpsp_main ] [ Functions ]
NAME
pawpsp_consistency
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
m_pawpsp
CHILDREN
SOURCE
5042 subroutine pawpsp_consistency() 5043 5044 5045 !This section has been created automatically by the script Abilint (TD). 5046 !Do not modify the following lines by hand. 5047 #undef ABI_FUNC 5048 #define ABI_FUNC 'pawpsp_consistency' 5049 !End of the abilint section 5050 5051 implicit none 5052 5053 !Local variables------------------------------- 5054 5055 ! ************************************************************************* 5056 5057 !Check pspcod=7 or 17 5058 if(pspcod/=7 .and. pspcod/=17)then 5059 write(msg, '(a,i2,a,a)' )& 5060 & 'In reading atomic psp file, finds pspcod=',pspcod,ch10,& 5061 & 'This is not an allowed value within PAW.' 5062 MSG_BUG(msg) 5063 end if 5064 5065 !Does nuclear charge znuclpsp agree with psp input znucl 5066 if (abs(znuclpsp-znucl)>tol8) then 5067 write(msg, '(a,f10.5,2a,f10.5,5a)' )& 5068 & 'Pseudopotential file znucl=',znucl,ch10,& 5069 & 'does not equal input znuclpsp=',znuclpsp,' better than 1e-08 .',ch10,& 5070 & 'znucl is read from the psp file in pspatm_abinit, while',ch10,& 5071 & 'znuclpsp is read in iofn2.' 5072 MSG_BUG(msg) 5073 end if 5074 5075 !Does nuclear charge zionpsp agree with psp input zion 5076 if (abs(zionpsp-zion)>tol8) then 5077 write(msg, '(a,f10.5,2a,f10.5,5a)' )& 5078 & 'Pseudopotential file zion=',zion,ch10,& 5079 & 'does not equal input zionpsp=',zionpsp,' better than 1e-08 .',ch10,& 5080 & 'zion is read from the psp file in pawpsp_main, while',ch10,& 5081 & 'zionpsp is read in iofn2.' 5082 MSG_BUG(msg) 5083 end if 5084 5085 !Check several choices for ixc against pspxc 5086 !ixc is from ABINIT code; pspxc is from atomic psp file 5087 if (ixc==0) then 5088 msg='Note that input ixc=0 => no xc is being used.' 5089 MSG_WARNING(msg) 5090 else if(ixc/=pspxc) then 5091 write(msg, '(a,i8,a,a,a,i8,a,a,a,a,a,a,a,a,a,a)' ) & 5092 & 'Pseudopotential file pspxc=',pspxc,',',ch10,& 5093 & 'not equal to input ixc=',ixc,'.',ch10,& 5094 & 'These parameters must agree to get the same xc ',ch10,& 5095 & 'in ABINIT code as in psp construction.',ch10,& 5096 & 'Action : check psp design or input file.',ch10,& 5097 & 'Assume experienced user. Execution will continue.',ch10 5098 MSG_WARNING(msg) 5099 end if 5100 5101 if (lloc>lmax ) then 5102 write(msg, '(a,2i12,a,a,a,a)' )& 5103 & 'lloc,lmax=',lloc,lmax,ch10,& 5104 & 'chosen l of local psp exceeds range from input data.',ch10,& 5105 & 'Action : check pseudopotential input file.' 5106 MSG_ERROR(msg) 5107 end if 5108 5109 end subroutine pawpsp_consistency