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-2024 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
912 subroutine psp11lo(drdi,epsatm,mmax,mqgrid,qgrid,q2vq,rad,vloc,yp1,ypn,zion) 913 914 !Arguments---------------------------------------------------------- 915 !scalars 916 integer,intent(in) :: mmax,mqgrid 917 real(dp),intent(in) :: zion 918 real(dp),intent(out) :: epsatm,yp1,ypn 919 !arrays 920 real(dp),intent(in) :: drdi(mmax) 921 real(dp),intent(in) :: qgrid(mqgrid),rad(mmax),vloc(mmax) 922 real(dp),intent(out) :: q2vq(mqgrid) 923 924 !Local variables------------------------------- 925 !scalars 926 integer :: iq,ir 927 real(dp),parameter :: scale=10.0d0 928 real(dp) :: arg,result_ctrap,test,ztor1 929 !arrays 930 real(dp),allocatable :: work(:) 931 932 ! ************************************************************************* 933 934 ABI_MALLOC(work,(mmax)) 935 936 ! Do q=0 separately (compute epsatm) 937 ! Do integral from 0 to r1 938 ztor1=(zion/2.0d0+rad(1)*vloc(1)/3.d0)*rad(1)**2 939 940 ! Set up integrand for q=0: $ \int[r^2 (V(r)+\frac{Zv}{r}) dr]$ 941 ! with extra factor of drdi to convert to uniform grid 942 do ir = 1, mmax 943 ! First handle tail region 944 test=vloc(ir)+zion/rad(ir) 945 ! Ignore small contributions, or impose a cut-off in the case 946 ! the pseudopotential data are in single precision. 947 ! (it is indeed expected that vloc is very close to zero beyond 20, 948 ! so a value larger than 2.0d-8 is considered anomalous) 949 if (abs(test)<1.0d-20 .or. (rad(ir)>20.0d0 .and. abs(test)>2.0d-8) ) then 950 work(ir)=zero 951 else 952 work(ir)=rad(ir)*(rad(ir)*vloc(ir)+zion) 953 end if 954 work(ir)=work(ir)*drdi(ir) 955 end do 956 957 ! Do integral from r(1) to r(max) 958 call ctrap(mmax,work,one,result_ctrap) 959 epsatm=4.d0*pi*(result_ctrap+ztor1) 960 961 q2vq(1)=-zion/pi 962 963 ! Loop over q values 964 do iq=2,mqgrid 965 arg=2.d0*pi*qgrid(iq) 966 ! ztor1=$ -Zv/\pi + 2q \int_0^{r1}[\sin(2\pi q r)(rV(r)+Zv) dr]$ 967 ztor1=(vloc(1)*sin(arg*rad(1))/arg-(rad(1)*vloc(1)+zion) * cos(arg*rad(1)) )/pi 968 969 ! set up integrand 970 do ir=1,mmax 971 !test=vloc(ir)+zion/rad(ir) 972 !Ignore contributions within decade of machine precision (suppressed ...) 973 !if ((scale+abs(test)).eq.scale) then 974 !work(ir)=zero 975 !else 976 work(ir)=sin(arg*rad(ir))*(rad(ir)*vloc(ir)+zion) 977 !end if 978 work(ir)=work(ir)*drdi(ir) 979 end do 980 ! do integral from r(1) to r(mmax) 981 call ctrap(mmax,work,one,result_ctrap) 982 983 ! store q^2 v(q) 984 ! FIXME: I only see one factor q, not q^2, but the same is done in other pspXlo.F90 985 q2vq(iq)=ztor1+2.d0*qgrid(iq)*result_ctrap 986 987 end do 988 989 ! Compute derivatives of q^2 v(q) at ends of interval 990 yp1=0.0d0 991 !ypn=$ 2\int_0^\infty[(\sin(2\pi qmax r)+(2\pi qmax r)*\cos(2\pi qmax r)(r V(r)+Z) dr]$ 992 !integral from 0 to r1 993 arg=2.0d0*pi*qgrid(mqgrid) 994 ztor1=zion*rad(1)*sin(arg*rad(1)) 995 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)) 996 !integral from r(mmax) to infinity is overkill; ignore 997 !set up integrand 998 do ir=1,mmax 999 !test=vloc(ir)+zion/rad(ir) 1000 !Ignore contributions within decade of machine precision (supressed ...) 1001 !if ((scale+abs(test)).eq.scale) then 1002 !work(ir)=0.0d0 1003 !else 1004 work(ir)=(sin(arg*rad(ir))+arg*rad(ir)*cos(arg*rad(ir))) * (rad(ir)*vloc(ir)+zion) 1005 !end if 1006 work(ir)=work(ir)*drdi(ir) 1007 end do 1008 1009 call ctrap(mmax,work,one,result_ctrap) 1010 ypn=2.0d0 * (ztor1 + result_ctrap) 1011 ABI_FREE(work) 1012 1013 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_spl(:), ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:) 126 127 ! ********************************************************************* 128 129 ! ######### in module pseudo: ############ 130 ! 131 ! only npsx = 1 is used here 132 ! grids are allocated for much larger fixed length (ndm=2000) 133 ! number of species (6) and projectors (8) as well... 134 ! 135 ! psd(npsx) = specied string 136 ! pseudotype = uspp / nc string 137 ! dft(npsx) = exchange correlation string (20 chars) 138 ! lmax(npsx) = maximum l channel 139 ! mesh(npsx) = number of points for local pot 140 ! nbeta(npsx) = number of projectors (beta functions for uspp) 141 ! nlcc(npsx) = flag for presence of NL core correction 142 ! zp(npsx) = valence ionic charge 143 ! r(ndm,npsx) = radial mesh 144 ! rab(ndm,npsx) = dr / di for radial mesh 145 ! rho_atc(ndm,npsx) = NLCC pseudocharge density 146 ! vloc0(ndm,npsx) = local pseudopotential 147 ! betar(ndm, nbrx, npsx) = projector functions in real space mesh 148 ! lll(nbrx,npsx) = angular momentum channel for each projector 149 ! ikk2(nbrx,npsx) = maximum index for each projector function 150 ! dion(nbrx,nbrx,npsx) = dij or Kleinman Bylander energies 151 ! 152 ! ######## end description of pseudo module contents ########## 153 154 if (open_file (filpsp,msg,newunit=iunit,status='old',form='formatted') /= 0) then 155 ABI_ERROR(msg) 156 end if 157 158 ! read in psp data to static data in pseudo module, for ipsx == 1 159 call read_pseudo(1,iunit) 160 close (iunit) 161 162 ! convert from Rydberg to Ha units 163 vloc0 = half * vloc0 164 dion = half * dion 165 166 ! if upf file is a USPP one, stop 167 if (pseudotype == 'US') then 168 ABI_ERROR('upf1_to_abinit: USPP UPF files not supported') 169 end if 170 171 ! copy over to abinit internal arrays and vars 172 ! FIXME: The API is broken. It does not recognize PBEsol 173 ! should use upfdft_to_ixc 174 call upfxc2abi(dft(1), pspxc) 175 lmax_ = lmax(1) 176 177 ! Check if the local component is one of the angular momentum channels 178 ! effectively if one of the ll is absent from the NL projectors 179 ABI_MALLOC(found_l, (0:lmax_)) 180 found_l = .true. 181 do ll = 0, lmax_ 182 if (any(lll(1:nbeta(1),1) == ll)) found_l(ll) = .false. 183 end do 184 185 if (count(found_l) /= 1) then 186 lloc = -1 187 else 188 do ll = 0, lmax_ 189 if (found_l(ll)) then 190 lloc = ll 191 exit 192 end if 193 end do 194 end if 195 196 ABI_FREE(found_l) 197 !FIXME: do something about lloc == -1 198 199 call atomdata_from_symbol(atom,psd(1)) 200 znucl = atom%znucl 201 zion = zp(1) 202 mmax = mesh(1) 203 204 call psp11lo(rab(1:mmax,1), epsatm, mmax, psps%mqgrid_vl, psps%qgrid_vl, & 205 vlspl(:,1), r(1:mmax,1), vloc0(1:mmax,1), yp1, ypn, zion) 206 207 208 ! Fit spline to q^2 V(q) (Numerical Recipes subroutine) 209 ABI_MALLOC(work_spl, (psps%mqgrid_vl)) 210 211 call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl) 212 213 vlspl(:,2) = work_spl(:) 214 ABI_FREE(work_spl) 215 216 ! this has to do the FT of the projectors to reciprocal space 217 ! allocate proj to avoid temporary copy. 218 ABI_MALLOC(proj, (mmax,1:nbeta(1))) 219 proj = betar(1:mmax,1:nbeta(1), 1) 220 221 call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, & 222 nbeta(1), proj, lll(1:nbeta(1),1), ikk2(1:nbeta(1),1), & 223 psps%qgrid_ff, r(1:mmax,1), rab(1:mmax,1), psps%useylm) 224 225 ABI_FREE(proj) 226 227 nproj_l = 0 228 do iproj = 1, nbeta(1) 229 ll = lll(iproj,1) 230 nproj_l(ll+1) = nproj_l(ll+1) + 1 231 end do 232 233 ! shape = dimekb vs. shape = n_proj 234 do ll = 1, nbeta(1) 235 ekb(ll) = dion(ll,ll,1) 236 end do 237 238 xcccrc = zero; xccc1d = zero 239 ! if we find a core density, do something about it 240 ! rho_atc contains the nlcc density 241 ! rho_at contains the total density 242 if (nlcc(1)) then 243 ABI_MALLOC(ff, (mmax)) 244 ABI_MALLOC(ff1, (mmax)) 245 ABI_MALLOC(ff2, (mmax)) 246 ff(1:mmax) = rho_atc(1:mmax,1) ! model core charge without derivative factor 247 248 ff1 = zero 249 call nderiv(one,ff,ff1,mmax,1) ! first derivative 250 ff1(1:mmax) = ff1(1:mmax) / rab(1:mmax,1) 251 call smooth(ff1, mmax, 15) ! run 15 iterations of smoothing 252 253 ff2 = zero 254 call nderiv(one, ff1, ff2, mmax, 1) ! second derivative 255 ff2(1:mmax) = ff2(1:mmax) / rab(1:mmax,1) 256 call smooth(ff2, mmax, 15) ! run 15 iterations of smoothing? 257 258 ! determine a good rchrg = xcccrc 259 do ir = mmax, 1, -1 260 if (abs(ff(ir)) > 1.e-6) then 261 xcccrc = r(ir,1) 262 exit 263 end if 264 end do 265 ABI_MALLOC(rad_cc, (mmax)) 266 rad_cc = r(1:mmax,1) 267 rad_cc(1) = zero ! force this so that the core charge covers whole spline interval. 268 269 call cc_derivatives(rad_cc,ff,ff1,ff2,mmax,psps%n1xccc,xcccrc,xccc1d) 270 271 ABI_FREE(rad_cc) 272 ABI_FREE(ff) 273 ABI_FREE(ff1) 274 ABI_FREE(ff2) 275 end if ! nlcc present 276 277 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 vloc_rcut= Real-space cutoff for local part 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
329 subroutine upf2_to_abinit(ipsp, filpsp, vloc_rcut, znucl, zion, pspxc, lmax, lloc, mmax, & 330 psps, epsatm, xcccrc, indlmn, ekb, ffspl, nproj_l, vlspl, xccc1d, nctab, maxrad) 331 332 use pseudo_types, only : pseudo_upf, deallocate_pseudo_upf !, pseudo_config 333 use read_upf_new_module, only : read_upf_new 334 use defs_datatypes, only : nctab_t 335 use m_psps, only : nctab_eval_tvalespl 336 use m_pspheads, only : upfdft_to_ixc, upf2_jl2srso 337 338 !Arguments ------------------------------- 339 integer,intent(in) :: ipsp 340 character(len=fnlen), intent(in) :: filpsp 341 real(dp),intent(in) :: vloc_rcut 342 type(pseudopotential_type),intent(in) :: psps 343 type(nctab_t),intent(inout) :: nctab 344 integer, intent(out) :: pspxc, lmax, lloc, mmax 345 real(dp), intent(out) :: znucl, zion, epsatm, xcccrc, maxrad 346 !arrays 347 integer, intent(out) :: indlmn(6,psps%lmnmax), nproj_l(psps%mpssoang) 348 real(dp), intent(inout) :: ekb(psps%dimekb) 349 real(dp), intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) 350 real(dp), intent(out) :: vlspl(psps%mqgrid_vl,2) 351 real(dp), intent(inout) :: xccc1d(psps%n1xccc,6) 352 353 !Local variables ------------------------- 354 integer :: ierr, ir, irad, iprj, il, ll, smooth_niter, nso, nn, iln, kk, mm, pspindex, iwfc !, iq 355 integer :: atmwfc_lmax, mmax_cut 356 !real(dp),parameter :: vloc_rcut = 10.0_dp ! QE Value 357 !real(dp),parameter :: vloc_rcut = 6.0_dp ! PseudoDojo value used to generate psp8 files 358 real(dp) :: yp1, ypn, amesh, damesh, intg 359 character(len=500) :: msg 360 logical :: linear_mesh, debug 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), units(2) 365 !integer :: nproj_tmp(psps%mpssoang) 366 integer,allocatable :: awfc_indlmn(:,:) 367 logical,allocatable :: found_l(:) 368 real(dp),allocatable :: work_spl(:), ff(:), ff1(:), ff2(:), rad_cc(:), proj(:,:), chi_tmp(:) 369 real(dp),allocatable :: vsr(:,:,:), esr(:,:), vso(:,:,:), eso(:,:) 370 371 ! ********************************************************************* 372 373 units = [std_out, ab_out] 374 375 ! See also https://github.com/QEF/qeschemas/blob/master/UPF/qe_pp-0.99.xsd 376 ! and https://github.com/QEF/qeschemas/files/9497267/pp.md 377 call read_upf_new(filpsp, upf, ierr) 378 ABI_CHECK(ierr == 0, sjoin("read_upf_new returned ierr:", itoa(ierr))) 379 380 call atomdata_from_symbol(atom, upf%psd) 381 znucl = atom%znucl 382 zion = upf%zp 383 mmax = upf%mesh 384 maxrad = upf%rmax 385 386 ABI_CHECK(upfdft_to_ixc(upf%dft, pspxc, msg) == 0, msg) 387 lmax = upf%lmax 388 389 ! Write some description of file 390 write(msg, '(3(a,1x))' ) '-',trim(upf%psd), trim(upf%generated) 391 call wrtout(units, msg) 392 write(msg,'(a,f9.5,f10.5,2x,a,t47,a)')'-',znucl,zion,trim(upf%date),'znucl, zion, pspdat' 393 call wrtout(units, msg) 394 write(msg, '(5(i0,1x),t47,a)' ) 12, pspxc, lmax, upf%lloc, mmax,'pspcod,pspxc,lmax,lloc,mmax' 395 call wrtout(units, msg) 396 397 ! Check that rad grid is linear starting at zero 398 linear_mesh = .True. 399 amesh = upf%r(2) - upf%r(1); damesh = zero 400 do irad=2,mmax-1 401 damesh = max(damesh, abs(upf%r(irad)+amesh-upf%r(irad+1))) 402 end do 403 linear_mesh = damesh < tol8 404 405 if (.not. linear_mesh .or. abs(upf%r(1)) > tol16) then 406 write(msg,'(3a)')& 407 'Assuming pseudized valence charge given on linear radial mesh starting at zero.',ch10,& 408 'Action: check your pseudopotential file.' 409 ABI_ERROR(msg) 410 end if 411 412 ! Check if the local component is one of the angular momentum channels 413 ! effectively if one of the ll is absent from the NL projectors 414 ABI_MALLOC(found_l, (0:lmax)) 415 found_l = .true. 416 do ll=0,lmax 417 if (any(upf%lll(1:upf%nbeta) == ll)) found_l(ll) = .false. 418 end do 419 420 if (count(found_l) /= 1) then 421 lloc = -1 422 else 423 do ll=0,lmax 424 if (found_l(ll)) then 425 lloc = ll 426 exit 427 end if 428 end do 429 end if 430 ABI_FREE(found_l) 431 !FIXME: do something about lloc == -1 432 433 ! convert vloc from Rydberg to Ha 434 upf%vloc = half * upf%vloc 435 436 ! ================================================= 437 ! This comment is from q-e/Modules/read_pseudo.F90 438 ! ================================================= 439 440 ! the radial grid is defined up to r(mesh) but we introduce 441 ! an auxiliary variable msh to limit the grid up to rcut=10 a.u. 442 ! This is used to cut off the numerical noise arising from the 443 ! large-r tail in cases like the integration of V_loc-Z/r 444 ! 445 ! NB: In QE, the default value for vloc_rcut is 10 Bohr. 446 call wrtout(std_out, sjoin(" Cutting radial-mesh for vloc using vloc_rcut:", ftoa(vloc_rcut), "(Bohr)")) 447 mmax_cut = mmax 448 do ir=1,upf%mesh 449 if (upf%r(ir) > vloc_rcut) then 450 mmax_cut = ir 451 ! msh is forced to be odd for simpson integration (maybe obsolete?) 452 mmax_cut = 2 * ( (mmax_cut + 1) / 2) - 1 453 exit 454 end if 455 end do 456 457 !write(std_out,*)" UPF file with mmax: ", mmax, " with r_max:", upf%r(mmax) 458 !write(std_out,*)" Using mmax_cut: ", mmax_cut, " with r_cut:", upf%r(mmax_cut) 459 460 ! Note mmax_cut here. 461 if (linear_mesh) then 462 call psp8lo(amesh, epsatm, mmax_cut, psps%mqgrid_vl, psps%qgrid_vl, & 463 vlspl(:,1), upf%r, upf%vloc, yp1, ypn, zion) 464 else 465 call psp11lo(upf%rab, epsatm, mmax_cut, psps%mqgrid_vl, psps%qgrid_vl,& 466 vlspl(:,1), upf%r, upf%vloc, yp1, ypn, zion) 467 end if 468 469 ! Fit spline to q^2 V(q) (Numerical Recipes subroutine) 470 ABI_MALLOC(work_spl, (psps%mqgrid_vl)) 471 472 call spline(psps%qgrid_vl, vlspl(:,1), psps%mqgrid_vl, yp1, ypn, work_spl) 473 474 vlspl(:,2) = work_spl(:) 475 ABI_FREE(work_spl) 476 477 debug = .False.!; debug = .True. 478 if (debug) then 479 write(std_out,*)'# Vlocal upf = ' 480 write(std_out,*)' amesh = ', amesh 481 write(std_out,*)' epsatm = ', epsatm 482 write(std_out,*)' mmax = ', mmax 483 write(std_out,*)' mqgrid = ', psps%mqgrid_vl 484 do ir = 1, psps%mqgrid_vl 485 write(std_out,*)' qgrid = ', ir, psps%qgrid_vl(ir) 486 enddo 487 do ir = 1, psps%mqgrid_vl 488 write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2) 489 enddo 490 write(std_out,*) 491 do ir = 1, mmax 492 write(std_out,*)' rad = ', upf%r(ir), upf%vloc(ir) 493 enddo 494 write(std_out,*) 495 write(std_out,*)' yp1 = ', yp1 496 write(std_out,*)' ypn = ', ypn 497 write(std_out,*)' zion = ', zion 498 !stop 499 end if 500 501 nproj_l = 0 502 503 if (.not. upf%has_so) then 504 do iprj=1,upf%nbeta 505 ll = upf%lll(iprj) 506 nproj_l(ll+1) = nproj_l(ll+1) + 1 507 end do 508 write(msg, '(a,*(i6))' ) ' nproj',nproj_l 509 call wrtout(units, msg) 510 511 ! shape = dimekb vs. shape = n_proj 512 ! convert from Rydberg to Ha 513 do ll=1,upf%nbeta 514 ekb(ll) = upf%dion(ll,ll) * half 515 end do 516 517 ! this has to do the FT of the projectors to reciprocal space 518 ! allocate proj to avoid temporary copy. 519 ABI_MALLOC(proj, (mmax,1:upf%nbeta)) 520 proj = upf%beta(1:mmax,1:upf%nbeta) 521 522 call psp11nl(ffspl, indlmn, mmax, psps%lnmax, psps%lmnmax, psps%mqgrid_ff, & 523 upf%nbeta, proj, upf%lll(1:upf%nbeta), upf%kbeta(1:upf%nbeta), & 524 psps%qgrid_ff, upf%r(1:mmax), upf%rab(1:mmax), psps%useylm) 525 526 ! This to reproduce psp8in version with linear meshes. 527 ! Compute Vanderbilt-KB form factors and fit splines 528 call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, & 529 psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj) 530 531 ABI_FREE(proj) 532 533 else 534 call upf2_jl2srso(upf, my_nproj_l, my_nprojso_l, vsr, esr, vso, eso) 535 536 write(msg, '(a,*(i6))' ) ' nproj',my_nproj_l 537 call wrtout(units, msg) 538 write(msg, '(5x,a)' ) "spin-orbit psp" 539 call wrtout(units, msg) 540 write(msg, '(5x,a,*(i6))' ) ' nprojso',my_nprojso_l 541 call wrtout(units, msg) 542 543 ABI_CALLOC(proj, (mmax, psps%lnmax)) 544 pspindex = 0; iln=0; indlmn(:,:)=0 545 nso = 2 546 547 if (psps%pspso(ipsp) == 0) then 548 write (msg, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,& 549 ' but internal variable pspso is 0' 550 ABI_COMMENT(msg) 551 nso = 1 552 end if 553 554 do nn=1,nso 555 !do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 556 ! ll = ipsang-(nn-1)*lmax-1 557 do il=1,lmax+1 558 ll = il - 1 559 if (nn == 1) then 560 ! SR part 561 do iprj=1,my_nproj_l(ll) 562 iln = iln + 1 563 ekb(iln) = esr(iprj, il) 564 proj(:,iln) = vsr(:, iprj, il) 565 nproj_l(il) = my_nproj_l(ll) 566 kk = iprj 567 do mm=1,2*ll*psps%useylm+1 568 pspindex = pspindex + 1 569 indlmn(1,pspindex) = ll 570 indlmn(2,pspindex) = mm-ll*psps%useylm-1 571 indlmn(3,pspindex) = kk 572 indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm 573 indlmn(5,pspindex) = iln 574 indlmn(6,pspindex) = nn 575 !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex 576 end do 577 end do 578 579 else 580 ! SOC part 581 if (ll == 0) cycle 582 do iprj=1,my_nprojso_l(ll) 583 iln = iln + 1 584 ekb(iln) = eso(iprj, il) 585 proj(:,iln) = vso(:,iprj,il) 586 ! Note ll in nproj_l i.e. the s channel in the SOC part is not included in nproj. 587 nproj_l(ll + psps%mpsang) = my_nprojso_l(ll) 588 kk = iprj !+ my_nproj_l(ll) 589 do mm=1,2*ll*psps%useylm+1 590 pspindex = pspindex + 1 591 indlmn(1,pspindex) = ll 592 indlmn(2,pspindex) = mm-ll*psps%useylm-1 593 indlmn(3,pspindex) = kk 594 indlmn(4,pspindex) = ll*ll+(1-psps%useylm)*ll+mm 595 indlmn(5,pspindex) = iln 596 indlmn(6,pspindex) = nn 597 !print *, "indlmn:", indlmn(:,pspindex), "pspindex:", pspindex 598 end do 599 end do 600 end if 601 602 end do ! il 603 end do ! nn 604 605 ! This to reproduce psp8in version with linear meshes. 606 ! Compute Vanderbilt-KB form factors and fit splines 607 call psp8nl(amesh, ffspl, indlmn, lmax, psps%lmnmax, psps%lnmax, mmax, & 608 psps%mqgrid_ff, psps%qgrid_ff, upf%r, proj) 609 610 ABI_FREE(proj) 611 ABI_FREE(vsr) 612 ABI_FREE(esr) 613 ABI_FREE(vso) 614 ABI_FREE(eso) 615 !ABI_WARNING("upf2_to_abinit: UPF2 with SOC") 616 end if 617 618 ! if we find a core density, do something about it 619 ! rho_atc contains the nlcc density 620 ! rho_at contains the total density 621 622 ! In Abinit, at least for the Troullier-Martins pseudopotential, 623 ! the pseudocore charge density and its derivatives (xccc1d) 624 ! are introduced in a linear grid. 625 ! This grid is normalized, so the radial coordinates run between 626 ! from 0 and 1 (from 0 to xcccrc, where xcccrc is the radius 627 ! where the pseudo-core becomes zero). 628 629 xcccrc = zero; xccc1d = zero 630 631 if (upf%nlcc) then 632 ABI_MALLOC(ff, (mmax)) 633 ABI_MALLOC(ff1, (mmax)) 634 ABI_MALLOC(ff2, (mmax)) 635 ! model core charge without derivative factor 636 ff(1:mmax) = upf%rho_atc(1:mmax) 637 !smooth_niter = 15 ! run 15 iterations of smoothing? 638 smooth_niter = 0 ! Don't smooth core charges to be consistent with the treatment done in psp8in 639 640 ff1 = zero 641 call nderiv(one, ff, ff1, mmax, 1) ! first derivative 642 ff1(1:mmax) = ff1(1:mmax) / upf%rab(1:mmax) 643 call smooth(ff1, mmax, smooth_niter) 644 645 ff2 = zero 646 call nderiv(one, ff1, ff2, mmax, 1) ! second derivative 647 ff2(1:mmax) = ff2(1:mmax) / upf%rab(1:mmax) 648 call smooth(ff2, mmax, smooth_niter) 649 650 ! determine a good rchrg = xcccrc 651 do ir = mmax, 1, -1 652 !if (abs(ff(ir)) > tol6) then 653 if (abs(ff(ir)) > tol20) then 654 xcccrc = upf%r(ir); exit 655 end if 656 end do 657 !xcccrc = upf%r(mmax) 658 659 ABI_MALLOC(rad_cc, (mmax)) 660 rad_cc = upf%r(1:mmax) 661 rad_cc(1) = zero ! force this so that the core charge covers whole spline interval. 662 663 call cc_derivatives(rad_cc, ff, ff1, ff2, mmax, psps%n1xccc, xcccrc, xccc1d) 664 !call psp8cc(mmax, psps%n1xccc, xcccrc, xccc1d) 665 666 ABI_FREE(rad_cc) 667 ABI_FREE(ff) 668 ABI_FREE(ff1) 669 ABI_FREE(ff2) 670 end if ! nlcc present 671 672 ! Read pseudo valence charge in real space on the linear mesh 673 ! and transform it to reciprocal space on a regular grid 674 ! TODO: Spline input data on linear mesh if not linear 675 ABI_MALLOC(ff, (mmax)) 676 ff = upf%rho_at(1:mmax) / four_pi 677 where (abs(upf%r) > tol16) 678 ff = ff / upf%r(1:mmax) ** 2 679 else where 680 ff = zero 681 end where 682 683 ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space. 684 call pawrad_init(mesh, mesh_size=mmax, mesh_type=1, rstep=amesh) 685 call nctab_eval_tvalespl(nctab, zion, mesh, ff, psps%mqgrid_vl, psps%qgrid_vl) 686 687 nctab%num_tphi = upf%nwfc 688 if (upf%nwfc > 0) then 689 ! Store atomic wavefunctions and metadata in nctab, then compute form factors for spline. 690 ! NB: lchi is the radial part of the KS equation, multiplied by r. 691 ABI_MALLOC(awfc_indlmn, (6, upf%nwfc)) 692 awfc_indlmn = huge(1) 693 atmwfc_lmax = -1 694 ABI_MALLOC(chi_tmp, (upf%mesh)) 695 696 do iwfc=1, upf%nwfc 697 !print *, "label: ", upf%els(iwfc); print *, "n", upf%nchi(iwfc) 698 !print *, "nn", upf%nn(iwfc); print *, "j", upf%jchi(iwfc) 699 !print *, "l", upf%lchi(iwfc); print *, "occ:", upf%oc(iwfc) 700 atmwfc_lmax = max(atmwfc_lmax, upf%lchi(iwfc)) 701 chi_tmp = upf%chi(:, iwfc) ** 2; call simp_gen(intg, chi_tmp, mesh) 702 !write(std_out, *)" wavefunction (before rescaling) integrates to: ",intg 703 upf%chi(:, iwfc) = upf%chi(:, iwfc) / sqrt(intg) 704 chi_tmp = upf%chi(:, iwfc) ** 2; call simp_gen(intg, chi_tmp, mesh) 705 !write(std_out, *)" wavefunction (after rescaling) integrates to: ",intg 706 707 ! NB: we only need ll (1), and iln (5) in psp8nl 708 awfc_indlmn(1, iwfc) = upf%lchi(iwfc) 709 awfc_indlmn(5, iwfc) = iwfc 710 end do 711 ABI_FREE(chi_tmp) 712 713 ! All this sfree/remalloc stuff is for handling memory in multi dataset mode! 714 ABI_REMALLOC(nctab%tphi_qspl, (psps%mqgrid_ff, 2, upf%nwfc)) 715 ABI_SFREE(nctab%tphi_n) 716 ABI_SFREE(nctab%tphi_l) 717 ABI_SFREE(nctab%tphi_occ) 718 ABI_SFREE(nctab%tphi_jtot) 719 call alloc_copy(upf%nchi, nctab%tphi_n) 720 call alloc_copy(upf%lchi, nctab%tphi_l) 721 call alloc_copy(upf%oc, nctab%tphi_occ) 722 nctab%has_jtot = upf%has_so 723 if (upf%has_so) call alloc_copy(upf%jchi, nctab%tphi_jtot) 724 725 !call psp8nl(amesh, nctab%tphi_qspl, awfc_indlmn, atmwfc_lmax, upf%nwfc, upf%nwfc, mmax, & 726 ! psps%mqgrid_ff, psps%qgrid_ff, upf%r, upf%chi) 727 728 !do iwfc=1, upf%nwfc 729 ! do iq=1,psps%mqgrid_ff 730 ! write(555, *) nctab%tphi_qspl(iq,:,iwfc) 731 ! end do 732 !end do 733 734 !do iwfc=1, upf%nwfc 735 ! where (abs(upf%r) > tol16) 736 ! upf%chi(:,iwfc) = upf%chi(:,iwfc) / upf%r 737 ! else where 738 ! upf%chi(:,iwfc) = zero 739 ! end where 740 !end do 741 742 call pawpsp_nl(nctab%tphi_qspl, awfc_indlmn, upf%nwfc, upf%nwfc, psps%mqgrid_ff, psps%qgrid_ff, mesh, upf%chi) 743 744 !do iwfc=1, upf%nwfc 745 ! do iq=1,psps%mqgrid_ff 746 ! write(666, *) nctab%tphi_qspl(iq,:,iwfc) 747 ! end do 748 !end do 749 750 ABI_FREE(awfc_indlmn) 751 !stop "nwfc" 752 end if ! upf%nwfc > 0 753 754 ABI_FREE(ff) 755 call pawrad_free(mesh) 756 call deallocate_pseudo_upf(upf) 757 758 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
791 subroutine psp11nl(ffspl,indlmn, mmax, lnmax, lmnmax, mqgrid, n_proj, & 792 proj, proj_l, proj_nr, qgrid, r, drdi, useylm) 793 794 !Arguments ------------------------------------ 795 !scalars 796 integer,intent(in) :: mmax, lnmax, lmnmax, mqgrid, useylm, n_proj 797 !arrays 798 integer, intent(in) :: proj_l(n_proj) 799 integer, intent(in) :: proj_nr(n_proj) 800 integer, intent(out) :: indlmn(6,lmnmax) 801 real(dp),intent(in) :: r(mmax) 802 real(dp),intent(in) :: drdi(mmax) 803 real(dp),intent(in) :: proj(mmax,n_proj) 804 real(dp),intent(in) :: qgrid(mqgrid) 805 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i 806 807 !Local variables------------------------------- 808 !scalars 809 integer,parameter :: bessorder = 0 ! never calculate derivatives of bessel functions 810 integer :: iproj, nr, ll, llold, ipsang, i_indlmn 811 integer :: iproj_1l, ir, iq, mm 812 real(dp) :: res, arg, besfact, dummy, dummy2 813 real(dp), allocatable :: work(:) 814 character(len=500) :: msg 815 816 !************************************************************************* 817 818 ffspl = zero; indlmn = 0; i_indlmn = 0 819 llold = -1; iproj_1l = 1 820 821 ! loop over all projectors 822 do iproj=1,n_proj 823 if (iproj > lmnmax) then 824 write(msg,'(a,2i0)') ' Too many projectors found. n_proj, lmnmax = ',n_proj, lmnmax 825 ABI_ERROR(msg) 826 end if 827 828 nr = proj_nr(iproj) 829 ABI_MALLOC(work, (nr)) 830 ll = proj_l(iproj) 831 if (ll < llold) then 832 ABI_ERROR('UPF projectors are not in order of increasing ll') 833 else if (ll == llold) then 834 iproj_1l = iproj_1l + 1 835 else 836 iproj_1l = 1 837 llold = ll 838 end if 839 840 ! determine indlmn for this projector (keep in UPF order and enforce that they are in increasing ll) 841 ! indlmn(6,lmnmax,ntypat) 842 ! For each type of psp, 843 ! array giving l,m,n,lm,ln,spin for i=ln (if useylm=0) 844 ! or i=lmn (if useylm=1) 845 ! NB: spin is used for NC pseudos with SOC term: 1 if scalar term (spin diagonal), 2 if SOC term. 846 do mm=1,2*ll*useylm+1 847 i_indlmn = i_indlmn + 1 848 indlmn(1, i_indlmn) = ll 849 indlmn(2, i_indlmn) = mm-ll*useylm-1 850 indlmn(3, i_indlmn) = iproj_1l 851 indlmn(4, i_indlmn) = ll*ll+(1-useylm)*ll+mm 852 indlmn(5, i_indlmn) = iproj 853 indlmn(6, i_indlmn) = 1 !spin? FIXME: to get j for relativistic cases 854 end do 855 856 ! FT projectors to reciprocal space q 857 do iq=1,mqgrid 858 arg = two_pi * qgrid(iq) 859 860 ! FIXME: add semianalytic form for integral from 0 to first point 861 do ir=1,nr 862 call jbessel(besfact, dummy, dummy2, ll, bessorder, arg*r(ir)) 863 work(ir) = drdi(ir) * besfact * proj(ir, iproj) * r(ir) !* r(ir) 864 end do 865 call ctrap (nr, work, one, res) 866 ffspl(iq, 1, iproj) = res 867 end do 868 ABI_FREE(work) 869 end do ! iproj 870 871 ! add derivative of ffspl(:,1,:) for spline interpolation later 872 ABI_MALLOC(work, (mqgrid)) 873 do ipsang=1,lnmax 874 call spline(qgrid,ffspl(:,1,ipsang),mqgrid,zero,zero,ffspl(:,2,ipsang)) 875 end do 876 ABI_FREE(work) 877 878 end subroutine psp11nl