TABLE OF CONTENTS
ABINIT/m_upf2abinit [ Modules ]
NAME
m_upf2abinit
FUNCTION
Procedures to read NC pseudos in UPF1/UPF2 format and convert data into the internal ABINIT representation.
COPYRIGHT
Copyright (C) 2009-2022 ABINIT group (MJV, MG, 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 .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 module m_upf2abinit 24 25 use defs_basis 26 use m_abicore 27 use m_errors 28 use m_atomdata 29 use m_splines 30 31 use defs_datatypes, only : pseudopotential_type 32 use m_io_tools, only : open_file 33 use m_numeric_tools, only : smooth, nderiv, ctrap 34 use m_copy, only : alloc_copy 35 use m_paw_numeric, only : jbessel => paw_jbessel 36 use m_pawpsp, only : pawpsp_nl 37 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free, simp_gen 38 use m_psptk, only : cc_derivatives, psp8lo, psp8nl 39 40 implicit none 41 42 private
ABINIT/psp11lo [ Functions ]
NAME
psp11lo
FUNCTION
Compute sine transform to transform from V(r) to q^2 V(q). Computes integrals on logarithmic grid using related uniform grid in exponent and corrected trapezoidal integration. Generalized from psp5lo for non-log grids using dr/di weights.
INPUTS
drdi=derivative of radial grid wrt index mmax=number of radial r grid points mqgrid=number of grid points in q from 0 to qmax. qgrid(mqgrid)=q grid values (bohr**-1). rad(mmax)=r grid values (bohr). vloc(mmax)=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=derivative of q^2 V(q) wrt q at q=0 and q=qmax (needed for spline fitter).
SOURCE
861 subroutine psp11lo(drdi,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion) 862 863 !Arguments---------------------------------------------------------- 864 !scalars 865 integer,intent(in) :: mmax,mqgrid 866 real(dp),intent(in) :: zion 867 real(dp),intent(out) :: epsatm,yp1,ypn 868 !arrays 869 real(dp),intent(in) :: drdi(mmax) 870 real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax) 871 real(dp),intent(out) :: q2vq(mqgrid) 872 873 !Local variables------------------------------- 874 !scalars 875 integer :: iq,ir 876 real(dp),parameter :: scale=10.0d0 877 real(dp) :: arg,result_ctrap,test,ztor1 878 !arrays 879 real(dp),allocatable :: work(:) 880 881 ! ************************************************************************* 882 883 ABI_MALLOC(work,(mmax)) 884 885 ! Do q=0 separately (compute epsatm) 886 ! Do integral from 0 to r1 887 ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2 888 889 ! Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$ 890 ! with extra factor of drdi to convert to uniform grid 891 do ir = 1, mmax 892 ! First handle tail region 893 test=vloc(ir)+zion/rad(ir) 894 ! Ignore small contributions, or impose a cut-off in the case 895 ! the pseudopotential data are in single precision. 896 ! (it is indeed expected that vloc is very close to zero beyond 20, 897 ! so a value larger than 2.0d-8 is considered anomalous) 898 if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then 899 work(ir)=zero 900 else 901 work(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion) 902 end if 903 work(ir)=work(ir)*drdi(ir) 904 end do 905 906 ! Do integral from r(1) to r(max) 907 call ctrap(mmax,work,one,result_ctrap) 908 epsatm=4.d0*pi*(result_ctrap+ztor1) 909 910 q2vq(1)=-zion/pi 911 912 ! Loop over q values 913 do iq=2,mqgrid 914 arg=2.d0*pi*qgrid(iq) 915 ! ztor1=$ -Zv/\pi + 2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$ 916 ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion) * cos(arg*rad(1)) )/pi 917 918 ! set up integrand 919 do ir=1,mmax 920 !test=vloc(ir)+zion/rad(ir) 921 !Ignore contributions within decade of machine precision (suppressed ...) 922 !if ((scale+abs(test)).eq.scale) then 923 !work(ir)=zero 924 !else 925 work(ir)=sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion) 926 !end if 927 work(ir)=work(ir)*drdi(ir) 928 end do 929 ! do integral from r(1) to r(mmax) 930 call ctrap(mmax,work,one,result_ctrap) 931 932 ! store q^2 v(q) 933 ! FIXME: I only see one factor q, not q^2, but the same is done in other pspXlo.F90 934 q2vq(iq)=ztor1+2.d0*qgrid(iq)*result_ctrap 935 936 end do 937 938 ! Compute derivatives of q^2 v(q) at ends of interval 939 yp1=0.0d0 940 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$ 941 !integral from 0 to r1 942 arg=2.0d0*pi*qgrid(mqgrid) 943 ztor1=zion*rad(1)*sin(arg*rad(1)) 944 ztor1=ztor1+ 3.d0*rad(1)*vloc(1)*cos(arg*rad(1))/arg + (rad(1)**2-1.0d0/arg**2)*vloc(1)*sin(arg*rad(1)) 945 !integral from r(mmax) to infinity is overkill; ignore 946 !set up integrand 947 do ir=1,mmax 948 !test=vloc(ir)+zion/rad(ir) 949 !Ignore contributions within decade of machine precision (supressed ...) 950 !if ((scale+abs(test)).eq.scale) then 951 !work(ir)=0.0d0 952 !else 953 work(ir)=(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * (rad(ir)*vloc(ir)+zion) 954 !end if 955 work(ir)=work(ir)*drdi(ir) 956 end do 957 958 call ctrap(mmax,work,one,result_ctrap) 959 ypn=2.0d0 * (ztor1 + result_ctrap) 960 ABI_FREE(work) 961 962 end subroutine psp11lo
ABINIT/upf1_to_abinit [ Functions ]
NAME
upf1_to_abinit
FUNCTION
This routine wraps a call to a PWSCF module, which reads in a UPF1 (PWSCF / Espresso) format pseudopotential, then transfers data to abinit internal variables. "UPF1 PWSCF format" (pspcod=11)
INPUTS
filpsp = name of file with UPF data psps = sturcture with global dimension data for pseudopotentials, header info ... used contents: psps%lmnmax psps%mqgrid_ff psps%mqgrid_vl psps%dimekb psps%n1xccc psps%qgrid_ff psps%qgrid_vl
OUTPUT
pspxc = index of xc functional for this pseudo lmax_ = maximal angular momentum lloc = local component chosen for pseudopotential mmax = maximum number of points in real space radial grid znucl = charge of species nucleus zion = valence charge epsatm = integral of local potential - coulomb potential of zion xcccrc = radius for non linear core correction ekb(dimekb)= Kleinman Bylander energies, see pspatm.F90 epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree) indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector; if any, spin-orbit components begin at l=mpsang+1 nproj= number of projectors for each channel xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives, from psp file (used in NC only)
SOURCE
97 subroutine upf1_to_abinit(filpsp, znucl, zion, pspxc, lmax_, lloc, mmax, & 98 psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj_l, vlspl, xccc1d) 99 100 101 use pseudo_pwscf ! pwscf module with all data explicit! 102 use m_read_upf_pwscf, only : read_pseudo 103 use m_pspheads, only : upfxc2abi 104 105 !Arguments ------------------------------- 106 character(len=fnlen), intent(in) :: filpsp 107 type(pseudopotential_type),intent(in) :: psps 108 integer, intent(out) :: pspxc, lmax_, lloc, mmax 109 real(dp), intent(out) :: znucl, zion 110 real(dp), intent(out) :: epsatm, xcccrc 111 !arrays 112 integer, intent(out) :: indlmn(6,psps%lmnmax) 113 integer, intent(out) :: nproj_l(psps%mpssoang) 114 real(dp), intent(inout) :: ekb(psps%dimekb) 115 real(dp), intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) 116 real(dp), intent(out) :: vlspl(psps%mqgrid_vl,2) 117 real(dp), intent(inout) :: xccc1d(psps%n1xccc,6) 118 119 !Local variables ------------------------- 120 integer :: ir, iproj, ll, iunit 121 real(dp) :: yp1, ypn 122 character(len=500) :: msg 123 type(atomdata_t) :: atom 124 logical, allocatable :: found_l(:) 125 real(dp), allocatable :: work_space(:),work_spl(:) 126 real(dp), allocatable :: ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:) 127 128 ! ********************************************************************* 129 130 ! ######### in module pseudo: ############ 131 ! 132 ! only npsx = 1 is used here 133 ! grids are allocated for much larger fixed length (ndm=2000) 134 ! number of species (6) and projectors (8) as well... 135 ! 136 ! psd(npsx) = specied string 137 ! pseudotype = uspp / nc string 138 ! dft(npsx) = exchange correlation string (20 chars) 139 ! lmax(npsx) = maximum l channel 140 ! mesh(npsx) = number of points for local pot 141 ! nbeta(npsx) = number of projectors (beta functions for uspp) 142 ! nlcc(npsx) = flag for presence of NL core correction 143 ! zp(npsx) = valence ionic charge 144 ! r(ndm,npsx) = radial mesh 145 ! rab(ndm,npsx) = dr / di for radial mesh 146 ! rho_atc(ndm,npsx) = NLCC pseudocharge density 147 ! vloc0(ndm,npsx) = local pseudopotential 148 ! betar(ndm, nbrx, npsx) = projector functions in real space mesh 149 ! lll(nbrx,npsx) = angular momentum channel for each projector 150 ! ikk2(nbrx,npsx) = maximum index for each projector function 151 ! dion(nbrx,nbrx,npsx) = dij or Kleinman Bylander energies 152 ! 153 ! ######## end description of pseudo module contents ########## 154 155 if (open_file (filpsp,msg,newunit=iunit,status='old',form='formatted') /= 0) then 156 ABI_ERROR(msg) 157 end if 158 159 ! read in psp data to static data in pseudo module, for ipsx == 1 160 call read_pseudo(1,iunit) 161 close (iunit) 162 163 ! convert from Rydberg to Ha units 164 vloc0 = half * vloc0 165 dion = half * dion 166 167 ! if upf file is a USPP one, stop 168 if (pseudotype == 'US') then 169 ABI_ERROR('upf1_to_abinit: USPP UPF files not supported') 170 end if 171 172 ! copy over to abinit internal arrays and vars 173 ! FIXME: The API is broken. It does not recognize PBEsol 174 ! should use upfdft_to_ixc 175 call upfxc2abi(dft(1), pspxc) 176 lmax_ = lmax(1) 177 178 ! Check if the local component is one of the angular momentum channels 179 ! effectively if one of the ll is absent from the NL projectors 180 ABI_MALLOC(found_l, (0:lmax_)) 181 found_l = .true. 182 do ll = 0, lmax_ 183 if (any(lll(1:nbeta(1),1) == ll)) found_l(ll) = .false. 184 end do 185 186 if (count(found_l) /= 1) then 187 lloc = -1 188 else 189 do ll = 0, lmax_ 190 if (found_l(ll)) then 191 lloc = ll 192 exit 193 end if 194 end do 195 end if 196 197 ABI_FREE(found_l) 198 !FIXME: do something about lloc == -1 199 200 call atomdata_from_symbol(atom,psd(1)) 201 znucl = atom%znucl 202 zion = zp(1) 203 mmax = mesh(1) 204 205 call psp11lo(rab(1:mmax,1), epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl, & 206 vlspl(:,1), r(1:mmax,1), vloc0(1:mmax,1), yp1, ypn, zion) 207 208 209 ! Fit spline to q^2 V(q) (Numerical Recipes subroutine) 210 ABI_MALLOC(work_space, (psps%mqgrid_vl)) 211 ABI_MALLOC(work_spl, (psps%mqgrid_vl)) 212 213 call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl) 214 215 vlspl(:,2) = work_spl(:) 216 ABI_FREE(work_space) 217 ABI_FREE(work_spl) 218 219 ! this has to do the FT of the projectors to reciprocal space 220 ! allocate proj to avoid temporary copy. 221 ABI_MALLOC(proj, (mmax,1:nbeta(1))) 222 proj = betar(1:mmax,1:nbeta(1), 1) 223 224 call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, & 225 nbeta(1), proj, lll(1:nbeta(1),1), ikk2(1:nbeta(1),1), & 226 psps%qgrid_ff, r(1:mmax,1), rab(1:mmax,1), psps%useylm) 227 228 ABI_FREE(proj) 229 230 nproj_l = 0 231 do iproj = 1, nbeta(1) 232 ll = lll(iproj,1) 233 nproj_l(ll+1) = nproj_l(ll+1) + 1 234 end do 235 236 ! shape = dimekb vs. shape = n_proj 237 do ll = 1, nbeta(1) 238 ekb(ll) = dion(ll,ll,1) 239 end do 240 241 xcccrc = zero; xccc1d = zero 242 ! if we find a core density, do something about it 243 ! rho_atc contains the nlcc density 244 ! rho_at contains the total density 245 if (nlcc(1)) then 246 ABI_MALLOC(ff, (mmax)) 247 ABI_MALLOC(ff1, (mmax)) 248 ABI_MALLOC(ff2, (mmax)) 249 ff(1:mmax) = rho_atc(1:mmax,1) ! model core charge without derivative factor 250 251 ff1 = zero 252 call nderiv(one,ff,ff1,mmax,1) ! first derivative 253 ff1(1:mmax) = ff1(1:mmax) / rab(1:mmax,1) 254 call smooth(ff1, mmax, 15) ! run 15 iterations of smoothing 255 256 ff2 = zero 257 call nderiv(one, ff1, ff2, mmax, 1) ! second derivative 258 ff2(1:mmax) = ff2(1:mmax) / rab(1:mmax,1) 259 call smooth(ff2, mmax, 15) ! run 15 iterations of smoothing? 260 261 ! determine a good rchrg = xcccrc 262 do ir = mmax, 1, -1 263 if (abs(ff(ir)) > 1.e-6) then 264 xcccrc = r(ir,1) 265 exit 266 end if 267 end do 268 ABI_MALLOC(rad_cc, (mmax)) 269 rad_cc = r(1:mmax,1) 270 rad_cc(1) = zero ! force this so that the core charge covers whole spline interval. 271 272 call cc_derivatives(rad_cc,ff,ff1,ff2,mmax,psps%n1xccc,xcccrc,xccc1d) 273 274 ABI_FREE(rad_cc) 275 ABI_FREE(ff) 276 ABI_FREE(ff1) 277 ABI_FREE(ff2) 278 end if ! nlcc present 279 280 end subroutine upf1_to_abinit
ABINIT/upf2_to_abinit [ Functions ]
NAME
upf2_to_abinit
FUNCTION
This routine wraps a call to a PWSCF module, which reads in a UPF2 (PWSCF / Espresso) format pseudopotential, then transfers data to abinit internal variables. "UPF2 PWSCF format" (pspcod=12)
INPUTS
filpsp = name of file with UPF2 data psps = sturcture with global dimension data for pseudopotentials, header info ... used contents: psps%lmnmax psps%mqgrid_ff psps%mqgrid_vl psps%dimekb psps%n1xccc psps%qgrid_ff psps%qgrid_vl
OUTPUT
pspxc = index of xc functional for this pseudo lmax_ = maximal angular momentum lloc = local component chosen for pseudopotential mmax = maximum number of points in real space radial grid znucl = charge of species nucleus zion = valence charge epsatm = integral of local potential - coulomb potential of zion xcccrc = radius for non linear core correction ekb(dimekb)= Kleinman Bylander energies, see pspatm.F90 epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$(hartree) indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit ffspl(mqgrid_ff,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector; if any, spin-orbit components begin at l=mpsang+1 nproj= number of projectors for each channel xccc1d(n1xccc*(1-usepaw),6)=1D core charge function and five derivatives, from psp file (used in NC only) nctab<nctab_t>=NC tables %has_tvale=True if the pseudo contains the pseudo valence charge %tvalespl(mqgrid_vl,2)=the pseudo valence density and 2nd derivative in reciprocal space on a regular grid
SOURCE
331 subroutine upf2_to_abinit(ipsp, filpsp, znucl, zion, pspxc, lmax, lloc, mmax, & 332 psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj_l, vlspl, xccc1d, nctab, maxrad) 333 334 use pseudo_types, only : pseudo_upf, deallocate_pseudo_upf !, pseudo_config 335 use read_upf_new_module, only : read_upf_new 336 use defs_datatypes, only : nctab_t 337 use m_psps, only : nctab_eval_tvalespl 338 use m_pspheads, only : upfdft_to_ixc, upf2_jl2srso 339 340 !Arguments ------------------------------- 341 integer,intent(in) :: ipsp 342 character(len=fnlen), intent(in) :: filpsp 343 type(pseudopotential_type),intent(in) :: psps 344 type(nctab_t),intent(inout) :: nctab 345 integer, intent(out) :: pspxc, lmax, lloc, mmax 346 real(dp), intent(out) :: znucl, zion 347 real(dp), intent(out) :: epsatm, xcccrc, maxrad 348 !arrays 349 integer, intent(out) :: indlmn(6,psps%lmnmax), nproj_l(psps%mpssoang) 350 real(dp), intent(inout) :: ekb(psps%dimekb) 351 real(dp), intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) 352 real(dp), intent(out) :: vlspl(psps%mqgrid_vl,2) 353 real(dp), intent(inout) :: xccc1d(psps%n1xccc,6) 354 355 !Local variables ------------------------- 356 integer :: ierr, ir, irad, iprj, il, ll, smooth_niter, nso, nn, iln, kk, mm, pspindex, iwfc !, iq 357 integer :: atmwfc_lmax 358 real(dp) :: yp1, ypn, amesh, damesh, intg 359 character(len=500) :: msg 360 logical :: linear_mesh 361 type(pseudo_upf) :: upf 362 type(atomdata_t) :: atom 363 type(pawrad_type) :: mesh 364 integer :: my_nproj_l(0:3), my_nprojso_l(1:3) 365 !integer :: nproj_tmp(psps%mpssoang) 366 integer,allocatable :: awfc_indlmn(:,:) 367 logical,allocatable :: found_l(:) 368 real(dp),allocatable :: work_space(:), work_spl(:) 369 real(dp),allocatable :: ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:) 370 real(dp),allocatable :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:) 371 372 ! ********************************************************************* 373 374 ! See also https://github.com/QEF/qeschemas/blob/master/UPF/qe_pp-0.99.xsd 375 ! and https://github.com/QEF/qeschemas/files/9497267/pp.md 376 call read_upf_new(filpsp, upf, ierr) 377 ABI_CHECK(ierr == 0, sjoin("read_upf_new returned ierr:", itoa(ierr))) 378 379 call atomdata_from_symbol(atom, upf%psd) 380 znucl = atom%znucl 381 zion = upf%zp 382 mmax = upf%mesh 383 ! FIXME Temporary hack 384 !mmax = 300 ! H 385 !mmax = 600 ! Si 386 maxrad = upf%rmax 387 !maxrad = upf%r(mmax) 388 389 ABI_CHECK(upfdft_to_ixc(upf%dft, pspxc, msg) == 0, msg) 390 lmax = upf%lmax 391 392 ! Write some description of file 393 write(msg, '(3(a,1x))' ) '-',trim(upf%psd), trim(upf%generated) 394 call wrtout([std_out, ab_out], msg) 395 write(msg,'(a,f9.5,f10.5,2x,a,t47,a)')'-',znucl,zion,trim(upf%date),'znucl, zion, pspdat' 396 call wrtout([std_out, ab_out], msg) 397 write(msg, '(5(i0,1x),t47,a)' ) 12, pspxc, lmax, upf%lloc, mmax,'pspcod,pspxc,lmax,lloc,mmax' 398 call wrtout([std_out, ab_out], msg) 399 400 ! Check that rad grid is linear starting at zero 401 linear_mesh = .True. 402 amesh = upf%r(2) - upf%r(1); damesh = zero 403 do irad=2,mmax-1 404 damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1))) 405 end do 406 linear_mesh = damesh < tol8 407 408 if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then 409 write(msg,'(3a)')& 410 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,& 411 'Action: check your pseudopotential file.' 412 ABI_ERROR(msg) 413 end if 414 415 ! Check if the local component is one of the angular momentum channels 416 ! effectively if one of the ll is absent from the NL projectors 417 ABI_MALLOC(found_l, (0:lmax)) 418 found_l = .true. 419 do ll=0,lmax 420 if (any(upf%lll(1:upf%nbeta) == ll)) found_l(ll) = .false. 421 end do 422 423 if (count(found_l) /= 1) then 424 lloc = -1 425 else 426 do ll=0,lmax 427 if (found_l(ll)) then 428 lloc = ll 429 exit 430 end if 431 end do 432 end if 433 ABI_FREE(found_l) 434 !FIXME: do something about lloc == -1 435 436 ! convert vloc from Rydberg to Ha 437 upf%vloc = half * upf%vloc 438 439 if (linear_mesh) then 440 call psp8lo(amesh, epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl, & 441 vlspl(:,1), upf%r(1:mmax), upf%vloc(1:mmax), yp1, ypn, zion) 442 else 443 call psp11lo(upf%rab(1:mmax), epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl,& 444 vlspl(:,1), upf%r(1:mmax), upf%vloc(1:mmax), yp1, ypn, zion) 445 end if 446 447 ! Fit spline to q^2 V(q) (Numerical Recipes subroutine) 448 ABI_MALLOC(work_space, (psps%mqgrid_vl)) 449 ABI_MALLOC(work_spl, (psps%mqgrid_vl)) 450 451 call spline(psps%qgrid_vl, vlspl(:,1), psps%mqgrid_vl, yp1, ypn, work_spl) 452 453 vlspl(:,2) = work_spl(:) 454 ABI_FREE(work_space) 455 ABI_FREE(work_spl) 456 457 nproj_l = 0 458 459 if (.not. upf%has_so) then 460 do iprj=1,upf%nbeta 461 ll = upf%lll(iprj) 462 nproj_l(ll+1) = nproj_l(ll+1) + 1 463 end do 464 write(msg, '(a,*(i6))' ) ' nproj',nproj_l 465 call wrtout([std_out, ab_out], msg) 466 467 ! shape = dimekb vs. shape = n_proj 468 ! convert from Rydberg to Ha 469 do ll=1,upf%nbeta 470 ekb(ll) = upf%dion(ll,ll) * half 471 end do 472 473 ! this has to do the FT of the projectors to reciprocal space 474 ! allocate proj to avoid temporary copy. 475 ABI_MALLOC(proj, (mmax,1:upf%nbeta)) 476 proj = upf%beta(1:mmax,1:upf%nbeta) 477 478 call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, & 479 upf%nbeta, proj, upf%lll(1:upf%nbeta), upf%kbeta(1:upf%nbeta), & 480 psps%qgrid_ff, upf%r(1:mmax), upf%rab(1:mmax), psps%useylm) 481 482 ! This to reproduce psp8in version with linear meshes. 483 ! Compute Vanderbilt-KB form factors and fit splines 484 call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, & 485 psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj) 486 487 ABI_FREE(proj) 488 489 else 490 call upf2_jl2srso(upf, my_nproj_l, my_nprojso_l, vsr, esr, vso, eso) 491 492 write(msg, '(a,*(i6))' ) ' nproj',my_nproj_l 493 call wrtout([std_out, ab_out], msg) 494 write(msg, '(5x,a)' ) "spin-orbit psp" 495 call wrtout([std_out, ab_out], msg) 496 write(msg, '(5x,a,*(i6))' ) ' nprojso',my_nprojso_l 497 call wrtout([std_out, ab_out], msg) 498 499 ABI_CALLOC(proj, (mmax, psps%lnmax)) 500 pspindex = 0; iln=0; indlmn(:,:)=0 501 nso = 2 502 503 if (psps%pspso(ipsp) == 0) then 504 write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,& 505 ' but internal variable pspso is 0' 506 ABI_COMMENT(msg) 507 nso = 1 508 end if 509 510 do nn=1,nso 511 !do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 512 ! ll = ipsang-(nn-1)*lmax-1 513 do il=1,lmax+1 514 ll = il - 1 515 if (nn == 1) then 516 ! SR part 517 do iprj=1,my_nproj_l(ll) 518 iln = iln + 1 519 ekb(iln) = esr(iprj, il) 520 proj(:,iln) = vsr(:, iprj, il) 521 nproj_l(il) = my_nproj_l(ll) 522 kk = iprj 523 do mm=1,2*ll*psps%useylm+1 524 pspindex = pspindex + 1 525 indlmn(1,pspindex) = ll 526 indlmn(2,pspindex) = mm-ll*psps%useylm-1 527 indlmn(3,pspindex) = kk 528 indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm 529 indlmn(5,pspindex) = iln 530 indlmn(6,pspindex) = nn 531 !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex 532 end do 533 end do 534 535 else 536 ! SOC part 537 if (ll == 0) cycle 538 do iprj=1,my_nprojso_l(ll) 539 iln = iln + 1 540 ekb(iln) = eso(iprj, il) 541 proj(:,iln) = vso(:,iprj,il) 542 ! Note ll in nproj_l i.e. the s channel in the SOC part is not included in nproj. 543 nproj_l(ll + psps%mpsang) = my_nprojso_l(ll) 544 kk = iprj !+ my_nproj_l(ll) 545 do mm=1,2*ll*psps%useylm+1 546 pspindex = pspindex + 1 547 indlmn(1,pspindex) = ll 548 indlmn(2,pspindex) = mm-ll*psps%useylm-1 549 indlmn(3,pspindex) = kk 550 indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm 551 indlmn(5,pspindex) = iln 552 indlmn(6,pspindex) = nn 553 !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex 554 end do 555 end do 556 end if 557 558 end do ! il 559 end do ! nn 560 561 ! This to reproduce psp8in version with linear meshes. 562 ! Compute Vanderbilt-KB form factors and fit splines 563 call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, & 564 psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj) 565 566 ABI_FREE(proj) 567 ABI_FREE(vsr) 568 ABI_FREE(esr) 569 ABI_FREE(vso) 570 ABI_FREE(eso) 571 !ABI_WARNING("upf2_to_abinit: UPF2 with SOC") 572 end if 573 574 ! if we find a core density, do something about it 575 ! rho_atc contains the nlcc density 576 ! rho_at contains the total density 577 xcccrc = zero; xccc1d = zero 578 579 if (upf%nlcc) then 580 ABI_MALLOC(ff, (mmax)) 581 ABI_MALLOC(ff1, (mmax)) 582 ABI_MALLOC(ff2, (mmax)) 583 ff(1:mmax) = upf%rho_atc(1:mmax) ! model core charge without derivative factor 584 !smooth_niter = 15 ! run 15 iterations of smoothing? 585 smooth_niter = 0 ! Don't smooth core charges to be consistent with the treatment done in psp8in 586 587 ff1 = zero 588 call nderiv(one, ff, ff1, mmax, 1) ! first derivative 589 ff1(1:mmax) = ff1(1:mmax) / upf%rab(1:mmax) 590 call smooth(ff1, mmax, smooth_niter) 591 592 ff2 = zero 593 call nderiv(one, ff1, ff2, mmax, 1) ! second derivative 594 ff2(1:mmax) = ff2(1:mmax) / upf%rab(1:mmax) 595 call smooth(ff2, mmax, smooth_niter) 596 597 ! determine a good rchrg = xcccrc 598 do ir = mmax, 1, -1 599 !if (abs(ff(ir)) > 1.e-6) then 600 if (abs(ff(ir)) > tol20) then 601 xcccrc = upf%r(ir) 602 exit 603 end if 604 end do 605 !xcccrc = upf%r(mmax) 606 607 ABI_MALLOC(rad_cc, (mmax)) 608 rad_cc = upf%r(1:mmax) 609 rad_cc(1) = zero ! force this so that the core charge covers whole spline interval. 610 611 call cc_derivatives(rad_cc, ff, ff1, ff2, mmax, psps%n1xccc, xcccrc, xccc1d) 612 !call psp8cc(mmax, psps%n1xccc, xcccrc, xccc1d) 613 614 ABI_FREE(rad_cc) 615 ABI_FREE(ff) 616 ABI_FREE(ff1) 617 ABI_FREE(ff2) 618 end if ! nlcc present 619 620 ! Read pseudo valence charge in real space on the linear mesh 621 ! and transform it to reciprocal space on a regular grid 622 ! TODO: Spline input data on linear mesh if not linear 623 ABI_MALLOC(ff, (mmax)) 624 ff = upf%rho_at / four_pi 625 where (abs(upf%r) > tol16) 626 ff = ff / upf%r ** 2 627 else where 628 ff = zero 629 end where 630 631 ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space. 632 call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh) 633 call nctab_eval_tvalespl(nctab, zion, mesh, ff, psps%mqgrid_vl, psps%qgrid_vl) 634 635 nctab%num_tphi = upf%nwfc 636 if (upf%nwfc > 0) then 637 ! Store atomic wavefunctions and metadata in nctab, then compute form factors for spline. 638 ! NB: lchi is the radial part of the KS equation, multiplied by r. 639 ABI_MALLOC(awfc_indlmn, (6, upf%nwfc)) 640 awfc_indlmn = huge(1) 641 atmwfc_lmax = -1 642 643 do iwfc=1, upf%nwfc 644 !print *, "label: ", upf%els(iwfc) 645 !print *, "n", upf%nchi(iwfc) 646 !print *, "nn", upf%nn(iwfc) 647 !print *, "j", upf%jchi(iwfc) 648 !print *, "l", upf%lchi(iwfc) 649 !print *, "occ:", upf%oc(iwfc) 650 atmwfc_lmax = max(atmwfc_lmax, upf%lchi(iwfc)) 651 ff = upf%chi(:, iwfc) ** 2; call simp_gen(intg, ff, mesh) 652 write(std_out, *)" wavefunction (before rescaling) integrates to: ",intg 653 upf%chi(:, iwfc) = upf%chi(:, iwfc) / sqrt(intg) 654 ff = upf%chi(:, iwfc) ** 2; call simp_gen(intg, ff, mesh) 655 write(std_out, *)" wavefunction (after rescaling) integrates to: ",intg 656 657 ! NB: we only need ll (1), and iln (5) in psp8nl 658 awfc_indlmn(1, iwfc) = upf%lchi(iwfc) 659 awfc_indlmn(5, iwfc) = iwfc 660 end do 661 662 ! All this sfree/remalloc stuff is for handling memory in multi dataset mode! 663 ABI_REMALLOC(nctab%tphi_qspl, (psps%mqgrid_ff, 2, upf%nwfc)) 664 ABI_SFREE(nctab%tphi_n) 665 ABI_SFREE(nctab%tphi_l) 666 ABI_SFREE(nctab%tphi_occ) 667 ABI_SFREE(nctab%tphi_jtot) 668 call alloc_copy(upf%nchi, nctab%tphi_n) 669 call alloc_copy(upf%lchi, nctab%tphi_l) 670 call alloc_copy(upf%oc, nctab%tphi_occ) 671 nctab%has_jtot = upf%has_so 672 if (upf%has_so) call alloc_copy(upf%jchi, nctab%tphi_jtot) 673 674 !call psp8nl(amesh, nctab%tphi_qspl, awfc_indlmn, atmwfc_lmax, upf%nwfc, upf%nwfc, mmax, & 675 ! psps%mqgrid_ff, psps%qgrid_ff, upf%r, upf%chi) 676 677 !do iwfc=1, upf%nwfc 678 ! do iq=1,psps%mqgrid_ff 679 ! write(555, *) nctab%tphi_qspl(iq,:,iwfc) 680 ! end do 681 !end do 682 683 !do iwfc=1, upf%nwfc 684 ! where (abs(upf%r) > tol16) 685 ! upf%chi(:,iwfc) = upf%chi(:,iwfc) / upf%r 686 ! else where 687 ! upf%chi(:,iwfc) = zero 688 ! end where 689 !end do 690 691 call pawpsp_nl(nctab%tphi_qspl, awfc_indlmn, upf%nwfc, upf%nwfc, psps%mqgrid_ff, psps%qgrid_ff, mesh, upf%chi) 692 693 !do iwfc=1, upf%nwfc 694 ! do iq=1,psps%mqgrid_ff 695 ! write(666, *) nctab%tphi_qspl(iq,:,iwfc) 696 ! end do 697 !end do 698 699 ABI_FREE(awfc_indlmn) 700 !stop "nwfc" 701 end if ! upf%nwfc > 0 702 703 ABI_FREE(ff) 704 call pawrad_free(mesh) 705 call deallocate_pseudo_upf(upf) 706 707 end subroutine upf2_to_abinit
m_upf2abinit/psp11nl [ Functions ]
[ Top ] [ m_upf2abinit ] [ Functions ]
NAME
psp11nl
FUNCTION
Fourier transform the real space UPF projector functions to reciprocal space
INPUTS
lmax=maximum ang momentum for which nonlocal form factor is desired. Usually lmax=1, sometimes = 0 (e.g. for oxygen); lmax <= 2 allowed. mmax=number of radial grid points for atomic grid lnmax= maximum index for all l channel projectors, dimension of ffspl lmnmax= maximum index for all projectors, dimension of indlmn mqgrid=number of grid points for q grid n_proj = total number of NL projectors read in proj = projector functions times r, on a real space grid proj_l = angular momentum channel for each projector proj_nr = max number of r-points used for each projector qgrid(mqgrid)=q-values at which form factors are returned r(mmax)=radial grid values drdi=derivative of grid point wrt index useylm = input to use m dependency of NL part, or only Legendre polynomials
OUTPUT
ffspl(mqgrid,2,mpsang)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum indlmn = indexing of each projector, for n, l, m, s, ln, lmn (see pspatm.F90)
SOURCE
740 subroutine psp11nl(ffspl,indlmn, mmax, lnmax, lmnmax, mqgrid, n_proj, & 741 proj, proj_l, proj_nr, qgrid, r, drdi, useylm) 742 743 !Arguments ------------------------------------ 744 !scalars 745 integer,intent(in) :: mmax, lnmax, lmnmax, mqgrid, useylm, n_proj 746 !arrays 747 integer, intent(in) :: proj_l(n_proj) 748 integer, intent(in) :: proj_nr(n_proj) 749 integer, intent(out) :: indlmn(6,lmnmax) 750 real(dp),intent(in) :: r(mmax) 751 real(dp),intent(in) :: drdi(mmax) 752 real(dp),intent(in) :: proj(mmax,n_proj) 753 real(dp),intent(in) :: qgrid(mqgrid) 754 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i 755 756 !Local variables------------------------------- 757 !scalars 758 integer,parameter :: bessorder = 0 ! never calculate derivatives of bessel functions 759 integer :: iproj, nr, ll, llold, ipsang, i_indlmn 760 integer :: iproj_1l, ir, iq, mm 761 real(dp) :: res, arg, besfact, dummy, dummy2 762 real(dp), allocatable :: work(:) 763 character(len=500) :: msg 764 765 !************************************************************************* 766 767 ffspl = zero; indlmn = 0; i_indlmn = 0 768 llold = -1; iproj_1l = 1 769 770 ! loop over all projectors 771 do iproj=1,n_proj 772 if (iproj > lmnmax) then 773 write(msg,'(a,2i0)') ' Too many projectors found. n_proj, lmnmax = ',n_proj, lmnmax 774 ABI_ERROR(msg) 775 end if 776 777 nr = proj_nr(iproj) 778 ABI_MALLOC(work, (nr)) 779 ll = proj_l(iproj) 780 if (ll < llold) then 781 ABI_ERROR('UPF projectors are not in order of increasing ll') 782 else if (ll == llold) then 783 iproj_1l = iproj_1l + 1 784 else 785 iproj_1l = 1 786 llold = ll 787 end if 788 789 ! determine indlmn for this projector (keep in UPF order and enforce that they are in increasing ll) 790 ! indlmn(6,lmnmax,ntypat) 791 ! For each type of psp, 792 ! array giving l,m,n,lm,ln,spin for i=ln (if useylm=0) 793 ! or i=lmn (if useylm=1) 794 ! NB: spin is used for NC pseudos with SOC term: 1 if scalar term (spin diagonal), 2 if SOC term. 795 do mm=1,2*ll*useylm+1 796 i_indlmn = i_indlmn + 1 797 indlmn(1, i_indlmn) = ll 798 indlmn(2, i_indlmn) = mm-ll*useylm-1 799 indlmn(3, i_indlmn) = iproj_1l 800 indlmn(4, i_indlmn) = ll*ll+(1-useylm)*ll+mm 801 indlmn(5, i_indlmn) = iproj 802 indlmn(6, i_indlmn) = 1 !spin? FIXME: to get j for relativistic cases 803 end do 804 805 ! FT projectors to reciprocal space q 806 do iq=1,mqgrid 807 arg = two_pi * qgrid(iq) 808 809 ! FIXME: add semianalytic form for integral from 0 to first point 810 do ir=1,nr 811 call jbessel(besfact, dummy, dummy2, ll, bessorder, arg*r(ir)) 812 work(ir) = drdi(ir) * besfact * proj(ir, iproj) * r(ir) !* r(ir) 813 end do 814 call ctrap (nr, work, one, res) 815 ffspl(iq, 1, iproj) = res 816 end do 817 ABI_FREE(work) 818 end do ! iproj 819 820 ! add derivative of ffspl(:,1,:) for spline interpolation later 821 ABI_MALLOC(work, (mqgrid)) 822 do ipsang=1,lnmax 823 call spline(qgrid,ffspl(:,1,ipsang),mqgrid,zero,zero,ffspl(:,2,ipsang)) 824 end do 825 ABI_FREE(work) 826 827 end subroutine psp11nl