TABLE OF CONTENTS
ABINIT/m_psp9 [ Modules ]
NAME
m_psp9
FUNCTION
Initialize pspcod=9 (pseudopotentials from the PSML XML format):
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (JJ, MVer, YP) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_psp9 23 24 use defs_basis 25 use m_splines 26 use m_errors 27 use m_abicore 28 #if defined HAVE_LIBPSML 29 use m_psml 30 #endif 31 32 use defs_datatypes, only : nctab_t 33 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 34 use m_psps, only : nctab_eval_tvalespl 35 use m_psptk, only : psp8lo, psp8nl 36 37 implicit none 38 39 private
ABINIT/psp9cc [ Functions ]
NAME
psp9cc
FUNCTION
Compute the core charge density, for use in the XC core correction, following the function definition valid for format 9 of the pseudopotentials (PSML).
INPUTS
mmax=maximum number of points in real space grid in the psp file n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
OUTPUT
rchrg=cut-off radius for the core density xccc1d(n1xccc,6)= 1D core charge function and its four first derivatives
NOTES
This routine will be built only if PSML support is enabled.
SOURCE
680 #if defined HAVE_LIBPSML 681 682 subroutine psp9cc(psxml,mmax,n1xccc,rad,rchrg,xccc1d) 683 684 !Arguments ------------------------------------ 685 !scalars 686 integer,intent(in) :: mmax,n1xccc 687 real(dp),intent(out) :: rchrg 688 type(ps_t),intent(in) :: psxml 689 !arrays 690 real(dp),intent(in) :: rad(mmax) 691 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 692 693 !Local variables------------------------------- 694 !scalars 695 integer :: i1xccc,idum,irad,jj 696 real(dp) :: amesh,c1,c2,c3,c4,damesh,dri,pi4i,tff,xp,xpm1,xpm2,xpp1,xx,twelvth 697 character(len=500) :: message 698 !arrays 699 integer :: iwork(8) 700 real(dp) :: rscale(5),dpoly(6,6),vpoly(6) 701 real(dp),allocatable :: ff(:,:) 702 703 !********************************************************************** 704 705 !Check that rad grid is linear starting at zero 706 amesh=rad(2)-rad(1) 707 damesh=zero 708 do irad=2,mmax-1 709 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 710 end do 711 712 if(damesh>tol8 .or. rad(1)/=zero) then 713 write(message, '(5a)' )& 714 & 'Pseudopotential input file requires linear radial mesh',ch10,& 715 & 'starting at zero.',ch10,& 716 & 'Action: check your pseudopotential input file.' 717 ABI_ERROR(message) 718 end if 719 720 ABI_MALLOC(ff,(mmax,5)) 721 722 dri = one / amesh 723 pi4i = quarter / pi 724 twelvth = one / 12.0_dp 725 726 !Read from pp file the model core charge and calculate its first 4 derivatives 727 !assumed to be on a linear grid starting at zero. 728 !The input functions contain the 4pi factor, and must be rescaled. 729 730 !Store the value of the pseudo-core charge. 731 ff(:,:) = zero 732 do jj=1,mmax 733 ff(jj,1) = ps_CoreCharge_Value(psxml,rad(jj)) 734 end do 735 736 !Calculate 4 first derivatives with 5-point stencil, except borders 737 do irad=3,mmax-2 738 ff(irad,2) = (-ff(irad+2,1) + 8.0d0*ff(irad+1,1) - & 739 & 8.0d0*ff(irad-1,1) + ff(irad-2,1)) * twelvth * dri 740 ff(irad,3) = (-ff(irad+2,1) + 16.0d0*ff(irad+1,1) - 30.0d0*ff(irad,1) + & 741 & 16.0d0*ff(irad-1,1) - ff(irad-2,1)) * twelvth * dri * dri 742 ff(irad,4) = (ff(irad+2,1) - 2.0d0*ff(irad+1,1) + & 743 & 2.0d0*ff(irad-1,1) - ff(irad-2,1)) * half * dri * dri * dri 744 ff(irad,5) = (ff(irad+2,1) - 4.0d0*ff(irad+1,1) + 6.0d0*ff(irad,1) - & 745 & 4.0d0*ff(irad-1,1) + ff(irad-2,1)) * dri * dri * dri * dri 746 end do 747 748 !Add border near zero using polynomial fit 749 dpoly(:,:) = zero 750 dpoly(:,1) = one 751 vpoly(:) = zero 752 vpoly(1) = ff(1,1) 753 do irad=2,6 754 do jj=1,6 755 dpoly(irad,jj) = rad(irad)**(jj-1) 756 end do 757 vpoly(irad) = ff(irad,1) 758 end do 759 call dgesv(6,1,dpoly,6,iwork,vpoly,6,idum) 760 761 do irad=1,2 762 ff(irad,2) = & 763 & vpoly(2) + 2.0d0*vpoly(3)*rad(irad) + & 764 & 3.0d0*vpoly(4)*rad(irad)*rad(irad) + & 765 & 4.0d0*vpoly(5)*rad(irad)*rad(irad)*rad(irad) + & 766 & 5.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad)*rad(irad) 767 ff(irad,3) = & 768 & 2.0d0*vpoly(3)*rad(irad) + & 769 & 6.0d0*vpoly(4)*rad(irad) + & 770 & 12.0d0*vpoly(5)*rad(irad)*rad(irad) + & 771 & 20.0d0*vpoly(6)*rad(irad)*rad(irad)*rad(irad) 772 ff(irad,4) = & 773 & 6.0d0*vpoly(4) + & 774 & 24.0d0*vpoly(5)*rad(irad) + & 775 & 60.0d0*vpoly(6)*rad(irad)*rad(irad) 776 ff(irad,5) = 24.0d0*vpoly(5) + & 777 & 120.0d0*vpoly(6)*rad(irad) 778 end do 779 780 !Make linear approximation for the tail near mmax 781 do irad=1,2 782 ff(mmax-2+irad,2) = ff(mmax-2,2) + irad * (ff(mmax-2,2) - ff(mmax-3,2)) 783 ff(mmax-2+irad,3) = ff(mmax-2,3) + irad * (ff(mmax-2,3) - ff(mmax-3,3)) 784 ff(mmax-2+irad,4) = ff(mmax-2,4) + irad * (ff(mmax-2,4) - ff(mmax-3,4)) 785 ff(mmax-2+irad,5) = ff(mmax-2,5) + irad * (ff(mmax-2,5) - ff(mmax-3,5)) 786 end do 787 788 !Renormalize core charge 789 ! ff(:,:) = ff(:,:) * pi4i 790 791 !determine xcccrc where the pseudocore becomes 0 792 !This is a difference with respect the Hamann's treatment of the core 793 !charge when reading PSP8. 794 !In Hamann's case (PSP8), xcccrc = rchrg, and this value is 795 !introduced in the pseudopotential input file. 796 !rchrg is not included in the PSML format 797 rchrg = zero 798 do jj=mmax,1,-1 799 if (ff(jj,1) > tol13) then 800 rchrg=rad(jj) 801 exit 802 end if 803 end do 804 805 !Check that input rchrg is consistent with last grid point 806 if(rchrg>rad(mmax)) then 807 write(message, '(5a)' )& 808 & 'Pseudopotential input file core charge mesh',ch10,& 809 & 'is inconsistent with rchrg in header.',ch10,& 810 & 'Action: check your pseudopotential input file.' 811 ABI_ERROR(message) 812 end if 813 814 !Factors for unit range scaling 815 do jj = 1, 5 816 rscale(jj)=rchrg**(jj-1) 817 end do 818 819 !Generate uniform mesh xx in the box cut by rchrg 820 !and interpolate the core charge and derivatives 821 !Cubic polynomial interpolation is used which is consistent 822 !with the original interpolation of these functions from 823 !a log grid to the input linear grid. 824 825 dri=1.d0/amesh 826 do i1xccc=1,n1xccc 827 xx=(i1xccc-1)* rchrg/dble(n1xccc-1) 828 829 ! index to find bracketing input mesh points 830 irad = int(dri * xx) + 1 831 irad = max(irad,2) 832 irad = min(irad,mmax-2) 833 ! interpolation coefficients 834 xp = dri * (xx - rad(irad)) 835 xpp1 = xp + one 836 xpm1 = xp - one 837 xpm2 = xp - two 838 c1 = -xp * xpm1 * xpm2 * sixth 839 c2 = xpp1 * xpm1 * xpm2 * half 840 c3 = - xp * xpp1 * xpm2 * half 841 c4 = xp * xpp1 * xpm1 * sixth 842 ! Now do the interpolation on all derivatives for this grid point 843 ! Include 1/4pi normalization and unit range scaling 844 do jj=1,5 845 tff = c1 * ff(irad - 1, jj) & 846 & + c2 * ff(irad , jj) & 847 & + c3 * ff(irad + 1, jj) & 848 & + c4 * ff(irad + 2, jj) 849 xccc1d(i1xccc,jj)=pi4i*rscale(jj)*tff 850 end do 851 end do 852 853 !5th derivative is apparently not in use, so set to zero 854 xccc1d(:,6)=zero 855 856 ABI_FREE(ff) 857 858 end subroutine psp9cc
ABINIT/psp9in [ Functions ]
NAME
psp9in
FUNCTION
Initialize pspcod=9 (pseudopotentials from the PSML XML format): continue to read the corresponding file, then compute the local and non-local potentials.
INPUTS
filpsp=filename of the PSML pseudopotential
lloc=angular momentum choice of local pseudopotential
lmax=value of lmax mentioned at the second line of the psp file
lmnmax=if useylm=1, max number of (l,m,n) comp. over all type of psps
=if useylm=0, max number of (l,n) comp. over all type of psps
lnmax=max. number of (l,n) components over all type of psps
mmax=maximum number of points in real space grid in the psp file
angular momentum of nonlocal pseudopotential
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mpssoang= 2*maximum angular momentum for nonlocal pseudopotentials - 1
mqgrid=dimension of q (or G) grid for arrays.
mqgrid_vl=dimension of q (or G) grid for valence charge (array qgrid_vl)
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
qgrid_vl(psps%mqgrid_vl)=values of q on grid from 0 to qmax (bohr^-1) for valence charge
pspso=spin-orbit characteristics, govern the content of ffspl and ekb
if =0 : this input requires NO spin-orbit characteristics of the psp
if =2 : this input requires HGH or psp8 characteristics of the psp
if =3 : this input requires HFN characteristics of the psp
useylm=governs the way the nonlocal operator is to be applied:
1=using Ylm, 0=using Legendre polynomials
zion=nominal valence of atom as specified in psp file
znucl=nuclear number of atom as specified in psp file
OUTPUT
ekb(lnmax)=Kleinman-Bylander energy, {{\ \begin{equation} \frac{\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r))^2 dr]} {\int_0^\infty [Rl(r)^2 (Vl(r)-Vloc(r)) dr]} \end{equation} }} for each (l,n) epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree) ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and second derivative from spline fit for each angular momentum and each projector indlmn(6,i)= array giving l,m,n,lm,ln,s for i=ln (if useylm=0) or i=lmn (if useylm=1) nproj(mpssoang)=number of projection functions for each angular momentum qchrg is not used, and could be suppressed later vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit xcccrc=XC core correction cutoff radius (bohr) xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file 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
106 subroutine psp9in(filpsp,ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 107 & mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,& 108 & useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad) 109 110 !Arguments ------------------------------------ 111 !scalars 112 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mpsang,mpssoang,mqgrid,mqgrid_vl 113 integer,intent(in) :: pspso,n1xccc,useylm 114 integer,intent(out) :: mmax 115 real(dp),intent(in) :: zion,znucl 116 real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad 117 type(nctab_t),intent(inout) :: nctab 118 character(len=fnlen),intent(in) :: filpsp 119 !arrays 120 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang) 121 real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl) 122 real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2) 123 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 124 125 !Local variables------------------------------- 126 !scalars 127 #if defined HAVE_LIBPSML 128 integer :: iln,pspindex,ipsang,irad,kk,ll 129 integer :: mm,nn,nso,ii,ir,il 130 integer :: nshells 131 integer :: iproj,irelt,nders 132 integer :: np_dn, np_lj, np_nr, np_so, np_sr, np_up, val_l, val_n 133 real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn,zval 134 real(dp) :: rmax,rmatch,z,chgvps 135 real(dp) :: val_occ 136 logical :: has_nlcc,has_spin 137 logical :: has_tvale,oncvpsp 138 character(len=500) :: message 139 character(len=30) :: creator 140 character(len=7), parameter :: oncvpsp_name = "ONCVPSP" 141 type(pawrad_type) :: mesh 142 #endif 143 !arrays 144 #if defined HAVE_LIBPSML 145 integer, allocatable :: idx_so(:),idx_sr(:) 146 real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_spl(:) 147 type(ps_t) :: psxml 148 #endif 149 150 ! *************************************************************************** 151 152 #if defined HAVE_LIBPSML 153 154 call ps_destroy(psxml) 155 call psml_reader(filpsp,psxml,debug=.true.) 156 157 !Identify the atomic code that generated the pseudopotential 158 call ps_Provenance_Get(psxml, 1, creator=creator) 159 !Check whether the pseudopotential has been created with ONCVPSP, 160 !Don Hamann's code 161 oncvpsp = (trim(creator(1:7)) .eq. trim(oncvpsp_name)) 162 !DEBUG 163 !write(std_out,*)' psp9in : creator : ', creator 164 !write(std_out,*)' psp9in : oncvpsp : ', oncvpsp 165 !ENDDEBUG 166 167 ! SIESTA's ATOM uses spherical harmonics, while ONCVPSP uses Legendre 168 ! polynomials, which means we have to check the consistency of input variables 169 ! wrt the pseudos 170 ! 171 ! Note: commented because NC pseudos do not have non-diagonal terms 172 ! 173 ! if ( oncvpsp ) then 174 ! if ( useylm /= 0 ) then 175 ! write(message,'(3a)') "ONCVPSP pseudos use Legendre polynomials but we use spherical harmonics", & 176 !& ch10, "ACTION: set useylm to 0 in your input file" 177 ! ABI_ERROR(message) 178 ! endif 179 ! else 180 ! if ( useylm == 0 ) then 181 ! write(message,'(3a)') "ATOM pseudos use spherical harmonics but we use Legendre polynomials", & 182 !& ch10, "ACTION: set useylm to 1 in your input file" 183 ! ABI_ERROR(message) 184 ! endif 185 ! endif 186 187 ! The atomic number is a real number instead of a simple integer 188 ! z (in Abinit), atomic-number in the header of the PSML file. 189 ! z = ps_AtomicNumber(psxml) 190 ! 191 ! The difference between the number of protons in the nucleus and the 192 ! sum of the populations of the core shells is the effective atomic number 193 ! of the pseudo-atom, Zval (in Abinit), z-pseudo in the header of the 194 ! PSML file. 195 ! zval = ps_Zpseudo(psxml) 196 197 call ps_PseudoAtomSpec_Get(psxml, & 198 & atomic_number=z, z_pseudo=zval, & 199 & spin_dft=has_spin, core_corrections=has_nlcc) 200 201 !--- 202 203 !Feb 2015: shifted to Hamann grid for convenience - libpsml interpolates anyway 204 ! 205 ! The following lines are taken from the oncvpsp.f90 subroutine of the oncvpsp 206 ! code implemented by D. Hamann 207 ! The atomic number of the element is read from the header of the XML file 208 ! Logarithmic grid defined by Hamann in oncvpsp code 209 ! z = psxml%header%z 210 ! amesh = 1.012d0 211 ! al = dlog(amesh) 212 ! rr1 = .0005d0/z 213 ! mmax = dlog(45.0d0 /rr1)/al 214 ! 215 ! ABI_MALLOC( rad,(mmax) ) 216 ! 217 ! do ir = 1, mmax 218 ! rad(ir) = rr1 * dexp(al*(ir-1)) 219 ! end do 220 221 !Determine the maximum number of points in the grid --- 222 rmax = 6.0_dp 223 amesh = 0.01_dp 224 mmax = int(rmax/amesh) 225 ! if(mod(mmax,2) .eq. 0) mmax = mmax + 1 226 227 !Print core charge info, for compatibility with psp8 228 rchrg = zero 229 fchrg = zero 230 if (has_nlcc) then 231 rchrg = amesh * (mmax - 2) 232 ! PSML does not store fchrg for now but we know we have core corrections, 233 ! then let's set it arbitrarily to 1.0 234 fchrg = one 235 else 236 write(message, '(a)' ) '- psp9in: No XC core correction.' 237 call wrtout(std_out,message,'COLL') 238 end if 239 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,zero,'rchrg,fchrg,qchrg' 240 call wrtout(ab_out,message,'COLL') 241 call wrtout(std_out,message,'COLL') 242 243 !Do we have a valence charge? 244 call ps_ValenceConfiguration_Get(psxml, nshells=nshells) 245 has_tvale = (nshells > 0) 246 247 ! Compute the valence charge of the reference configuration used to 248 ! generate the pseudopotential 249 chgvps = 0.0_dp 250 ! Loop on all the shells included in the valence 251 do il = 1, nshells 252 ! Sum the corresponding occupation of each shell 253 ! FIXME: What if there is spin? 254 call ps_ValenceShell_Get(psxml, il, n=val_n, l=val_l, occupation=val_occ) 255 chgvps = chgvps + val_occ 256 write(std_out,*)' psp9in : n, l, occupation = ', & 257 & val_n, val_l, val_occ 258 end do 259 260 !DEBUG 261 !write(std_out,*)' psp9in : atomic number' 262 !write(std_out,*)' psp9in : z = ', z 263 !write(std_out,*)' psp9in : valence charge of the reference configuration' 264 !write(std_out,*)' psp9in : chgvps = ', chgvps 265 !write(std_out,*)' psp9in : nominal valence charge' 266 !write(std_out,*)' psp9in : zval = ', zval 267 !write(std_out,*)' psp9in : mqgrid_vl = ', mqgrid_vl 268 !write(std_out,*)' psp9in : parameters to define the points of the grid' 269 !write(std_out,*)' psp9in : amesh = ', amesh 270 !write(std_out,*)' psp9in : rmax = ', rmax 271 !write(std_out,*)' psp9in : mmax = ', mmax 272 !ENDDEBUG 273 274 ! TODO: should be simple to average these and get difference for SREL+SOC, 275 ! but also the Ekb etc... 276 call ps_NonlocalProjectors_Filter(psxml, set=SET_DOWN, number=np_dn) 277 call ps_NonlocalProjectors_Filter(psxml, set=SET_LJ, number=np_lj) 278 call ps_NonlocalProjectors_Filter(psxml, set=SET_NONREL, number=np_nr) 279 call ps_NonlocalProjectors_Filter(psxml, set=SET_SO, number=np_so) 280 call ps_NonlocalProjectors_Filter(psxml, set=SET_SREL, number=np_sr) 281 call ps_NonlocalProjectors_Filter(psxml, set=SET_UP, number=np_up) 282 if (np_lj > 0) then 283 message = 'For the moment LJ format projectors are not supported; SREL + SO is the internal abinit format' 284 ABI_BUG(message) 285 end if 286 287 if (np_up > 0 .or. np_dn > 0) then 288 write (message,'(3a)') 'For the moment separate spin up and down format projectors are not supported;',ch10,& 289 & ' spin average is the internal abinit format' 290 ABI_BUG(message) 291 end if 292 293 !-------------------------------------------------------------------- 294 295 !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn 296 if(pspso==2) then 297 nso=2 298 else 299 nso=1 300 end if 301 302 !Find the number of projectors per angular momentum shell 303 nproj(:)=0 304 if (np_nr > 0) then 305 call ps_NonlocalProjectors_Filter(psxml, set=SET_NONREL, indexes=idx_sr) 306 do iproj = 1, np_nr 307 call ps_Projector_Get(psxml, idx_sr(iproj), l=il) 308 nproj(il+1) = nproj(il+1) + 1 309 end do 310 else 311 if (np_sr > 0) then 312 call ps_NonlocalProjectors_Filter(psxml, set=SET_SREL, indexes=idx_sr) 313 do iproj = 1, np_sr 314 call ps_Projector_Get(psxml, idx_sr(iproj), l=il) 315 nproj(il+1) = nproj(il+1) + 1 316 end do 317 else ! this should not happen 318 ABI_BUG('Your psml potential should have either scalar- or non- relativistic projectors') 319 end if 320 end if 321 322 write(message, '(a,5i6)' ) ' nproj',nproj(1:lmax+1) 323 call wrtout(ab_out,message,'COLL') 324 call wrtout(std_out, message,'COLL') 325 326 irelt = 0 327 if (nso == 2) then 328 call ps_NonlocalProjectors_Filter(psxml, set=SET_SO, indexes=idx_so) 329 do iproj = 1, np_so 330 call ps_Projector_Get(psxml, idx_so(iproj), l=il) 331 nproj(il+lmax+2) = nproj(il+lmax+2) + 1 332 irelt = 1 333 end do 334 end if 335 336 pspindex=0;iln=0;indlmn(:,:)=0 337 do nn=1,nso 338 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 339 ll=ipsang-(nn-1)*lmax-1 340 if (nproj(ipsang)>0) then 341 do kk=1,nproj(ipsang) 342 iln=iln+1 343 do mm=1,2*ll*useylm+1 344 pspindex=pspindex+1 345 indlmn(1,pspindex)=ll ! l angular momentum channel 346 indlmn(2,pspindex)=mm-ll*useylm-1 ! hash of position in m 347 indlmn(3,pspindex)=kk ! index of projector 348 indlmn(4,pspindex)=ll*ll+(1-useylm)*ll+mm ! hash of position in l(l+1) array 349 indlmn(5,pspindex)=iln ! absolute index of l, n disregarding m values 350 indlmn(6,pspindex)=nn ! spin orbit index!!! NOT the n shell index 351 end do 352 end do 353 end if 354 end do 355 end do 356 357 ! Determine whether the atomic calculation to generate the pseudopotential 358 ! is relativistic or not 359 360 !DEBUG 361 !write(std_out,*)' psp9in : pseudopotential generation relativity ', ps_Relativity(psxml) 362 !write(std_out,*)' psp9in : SOC pseudopotential? (1=yes, 0 =no) ' 363 !write(std_out,*)' psp9in : irelt = ', irelt 364 !write(ab_out,*)' psp9in : irelt = ', irelt 365 !ENDDEBUG 366 367 !Can now allocate grids, potentials and projectors 368 ABI_MALLOC(rad,(mmax)) 369 ABI_MALLOC(vloc,(mmax)) 370 ABI_MALLOC(vpspll,(mmax,lnmax)) 371 372 !Feb 2015: shifted to Hamann grid for convenience - libpsml interpolates anyway 373 do ir=1,mmax 374 rad(ir) = amesh * (ir - 1) 375 end do 376 !! DEBUG 377 ! do ir = 2, mmax 378 ! write(std_out,'(i5,f20.12)')ir, rad(ir) 379 ! end do 380 !! ENDDEBUG 381 !--- 382 write(message, '(a,i5,es16.6,es16.6)')' psp9in : mmax, amesh, rad(mmax) = ', mmax, amesh, rad(mmax) 383 call wrtout(ab_out,message,'COLL') 384 call wrtout(std_out,message,'COLL') 385 386 !Check that rad grid is linear starting at zero 387 amesh=rad(2)-rad(1) 388 damesh=zero 389 do irad=2,mmax-1 390 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 391 end do 392 if(damesh>tol8 .or. rad(1)/=zero) then 393 write(message, '(5a)' )& 394 & 'Pseudopotential input file requires linear radial mesh',ch10,& 395 & 'starting at zero.',ch10,& 396 & 'Action: check your pseudopotential input file.' 397 ABI_ERROR(message) 398 end if 399 400 !Take care of the non-linear core corrections 401 !---------------------------------------------------------------------------- 402 ! xcccrc : XC core correction cutoff radius (bohr) 403 ! It is defined as the radius where the pseudo-core 404 ! charge density becomes zero 405 ! (here we have set up a tolerance of 1.d-12). 406 407 rmatch = zero 408 nders = 0 409 if (has_nlcc) then 410 411 ! In Abinit, at least for the Troullier-Martins pseudopotential, 412 ! the pseudocore charge density and its derivatives (xccc1d) 413 ! are introduced in a linear grid. 414 ! This grid is normalized, so the radial coordinates run between 415 ! from 0 and 1 (from 0 to xcccrc, where xcccrc is the radius 416 ! where the pseudo-core becomes zero). 417 418 call ps_CoreCharge_get(psxml, rc=rmatch, nderivs=nders) 419 write (message,'(1X,A,A,5X,A,1X,F8.3,A,5X,A,I8,A)') & 420 & "Reading pseudocore charge",ch10, & 421 & "- matching radius:",rmatch,ch10, & 422 & "- number of continuous derivatives",nders,ch10 423 call wrtout(std_out,message,'COLL') 424 425 !Get core charge function and derivatives, if needed 426 if(fchrg>1.0d-15)then 427 call psp9cc(psxml,mmax,n1xccc,rad,rchrg,xccc1d) 428 ! The core charge function for pspcod=9 429 ! becomes zero beyond rchrg. Thus xcccrc must be set 430 ! equal to rchrg. 431 xcccrc=rchrg 432 else 433 xccc1d(:,:) = zero 434 xcccrc = zero 435 fchrg = zero 436 qchrg = zero 437 end if 438 439 maxrad = rad(mmax) 440 441 end if ! has_nlcc 442 443 !! DEBUG 444 ! write(std_out,*)' xcccrc = ', xcccrc, rchrg 445 ! write(std_out,*) 446 ! write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc 447 ! do ii = 1, n1xccc 448 ! write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),& 449 ! & xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6) 450 ! enddo 451 ! write(std_out,*) 452 ! stop 453 !! ENDDEBUG 454 455 456 !-------------------------------------------------------------------- 457 !Carry out calculations for local (lloc) pseudopotential. 458 !Obtain Fourier transform (1-d sine transform) 459 !to get q^2 V(q). 460 461 !Read and process vlocal: 462 !The local potential is given by a <radfunc> element under the <local-potential> 463 !element. 464 !After reading, this is a copy of the treatment to the 465 !local part carry out in psp8 466 !i.e. (as in Hamann pseudopotential) 467 ! 468 !Read the local component of the pseudopotential 469 vloc = zero 470 do ir = 1, mmax 471 vloc(ir) = ps_LocalPotential_Value(psxml, rad(ir)) 472 end do 473 474 call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,& 475 & vlspl(:,1),rad,vloc,yp1,ypn,zion) 476 477 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 478 ABI_MALLOC(work_spl,(mqgrid)) 479 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 480 vlspl(:,2)=work_spl(:) 481 ABI_FREE(work_spl) 482 483 !! DEBUG 484 ! write(std_out,*)'# Vlocal = ' 485 ! write(std_out,*)' amesh = ', amesh 486 ! write(std_out,*)' epsatm = ', epsatm 487 ! write(std_out,*)' mmax = ', mmax 488 ! write(std_out,*)' mqgrid = ', mqgrid 489 ! do ir = 1, mqgrid 490 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 491 ! enddo 492 ! do ir = 1, mqgrid 493 ! write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2) 494 ! enddo 495 ! write(std_out,*) 496 ! do ir = 1, mmax 497 ! write(std_out,*)' rad = ', rad(ir), vloc(ir) 498 ! enddo 499 ! write(std_out,*) 500 ! write(std_out,*)' yp1 = ', yp1 501 ! write(std_out,*)' ypn = ', ypn 502 ! write(std_out,*)' zion = ', zion 503 ! stop 504 !! ENDDEBUG 505 506 507 !-------------------------------------------------------------------- 508 !Take care of non-local part 509 510 !Zero out all Kleinman-Bylander energies to initialize 511 do ii = 1, lmnmax ! loop over all possible projectors 512 if (indlmn(6,ii) == 1) then 513 call ps_Projector_Get(psxml, idx_sr(indlmn(5,ii)), ekb=ekb(indlmn(5,ii))) 514 else if (indlmn(6,ii) == 2) then 515 call ps_Projector_Get(psxml, idx_so(indlmn(5,ii)), ekb=ekb(indlmn(5,ii))) 516 end if 517 end do 518 519 !Read the KB projectors from the PSML file 520 !Note than in the PSML file the radial part of the projector is stored, 521 !while Abinit expects the radial part times the radii. 522 !We have to multiply by r after reading it. 523 !Note than in Hamann's format (psp8), Abinit directly reads r * radial_part_KB 524 vpspll = zero 525 do ii = 1, lmnmax 526 if (indlmn(6,ii) == 1) then 527 do ir = 1, mmax 528 vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_sr(indlmn(5,ii)), rad(ir)) 529 vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii)) 530 end do 531 else if (indlmn(6,ii) == 2) then 532 do ir = 1, mmax 533 vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_so(indlmn(5,ii)), rad(ir)) 534 vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii)) 535 end do 536 end if 537 end do 538 539 !Allow for option of no nonlocal corrections (lloc=lmax=0) 540 if (lloc==0.and.lmax==0) then 541 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 542 call wrtout(ab_out,message,'COLL') 543 call wrtout(std_out,message,'COLL') 544 else 545 546 ! ---------------------------------------------------------------------- 547 ! Compute Vanderbilt-KB form factors and fit splines 548 549 call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,& 550 & mqgrid,qgrid,rad,vpspll) 551 552 end if 553 554 !! DEBUG 555 ! write(std_out,*)'# KB Projectors = ' 556 ! write(std_out,*)' amesh = ', amesh 557 ! do ir = 1, mqgrid 558 ! do il = 1, lnmax 559 ! write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il) 560 ! enddo 561 ! enddo 562 ! do il = 1, lmnmax 563 ! write(std_out,*)' indlmn = ', il, indlmn(:,il) 564 ! enddo 565 ! write(std_out,*)' lmax = ', lmax 566 ! write(std_out,*)' lmnmax = ', lmnmax 567 ! write(std_out,*)' lnmax = ', lnmax 568 ! write(std_out,*)' mmax = ', mmax 569 ! write(std_out,*)' mqgrid = ', mqgrid 570 ! do ir = 1, mqgrid 571 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 572 ! enddo 573 ! do il = 1, lnmax 574 ! write(std_out,*) 575 ! write(std_out,*)'# il = ', il 576 ! do ir = 1, mmax 577 ! write(std_out,*)' rad = ', rad(ir), vpspll(ir,il) 578 ! enddo 579 ! enddo 580 ! stop 581 !! ENDDEBUG 582 583 ! Read pseudo valence charge in real space on the linear mesh 584 ! and transform it to reciprocal space on a regular grid. Use vloc as workspace. 585 vloc(:) = zero 586 if (has_tvale) then 587 do irad=1,mmax 588 vloc(irad) = ps_ValenceCharge_Value(psxml,rad(irad)) 589 vloc(irad) = vloc(irad) / four_pi 590 end do 591 592 !! DEBUG 593 ! do irad = 1, mmax 594 ! write(std_out,*)' Valence Charge = ', rad(irad), vloc(irad) 595 ! enddo 596 ! stop 597 !! ENDDEBUG 598 599 600 ! Check that rad grid is linear starting at zero 601 amesh=rad(2)-rad(1) 602 damesh=zero 603 do irad=2,mmax-1 604 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 605 end do 606 if(damesh>tol8 .or. rad(1)/=zero) then 607 write(message, '(5a)' )& 608 & 'Pseudopotential input file requires linear radial mesh',ch10,& 609 & 'starting at zero.',ch10,& 610 & 'Action: check your pseudopotential input file.' 611 ABI_ERROR(message) 612 end if 613 614 ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space. 615 call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh) 616 call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl) 617 call pawrad_free(mesh) 618 end if 619 620 ABI_FREE(vpspll) 621 ABI_FREE(vloc) 622 ABI_FREE(rad) 623 if (allocated(idx_sr)) then 624 ABI_FREE_NOCOUNT(idx_sr) 625 end if 626 if (allocated(idx_so)) then 627 ABI_FREE_NOCOUNT(idx_so) 628 end if 629 630 call ps_destroy(psxml) 631 632 !-------------------------------------------------------------------- 633 634 #else 635 ABI_UNUSED(mpsang) 636 ABI_UNUSED(pspso) 637 ABI_UNUSED(qgrid_vl) 638 ABI_UNUSED(nctab%mqgrid_vl) 639 !Initialize some arguments, for portability at compile time 640 indlmn=0 ; mmax=0 ; nproj=0 641 ekb=zero ; epsatm=zero ; ffspl=zero ; qchrg=zero ; vlspl=zero ; xcccrc=zero ; xccc1d=zero 642 643 if(.false.)write(std_out,*)filpsp ! Just to keep filpsp when HAVE_LIBPSML is false 644 if(.false.)write(std_out,*)lloc ! Just to keep lloc when HAVE_LIBPSML is false 645 if(.false.)write(std_out,*)lmax ! Just to keep lmax when HAVE_LIBPSML is false 646 if(.false.)write(std_out,*)mpsang ! Just to keep mpsang when HAVE_LIBPSML is false 647 if(.false.)write(std_out,*)pspso ! Just to keep pspso when HAVE_LIBPSML is false 648 if(.false.)write(std_out,*)qgrid ! Just to keep qgrid when HAVE_LIBPSML is false 649 if(.false.)write(std_out,*)qgrid_vl ! Just to keep qgrid_vl when HAVE_LIBPSML is false 650 if(.false.)write(std_out,*)useylm ! Just to keep useylm when HAVE_LIBPSML is false 651 if(.false.)write(std_out,*)zion ! Just to keep zion when HAVE_LIBPSML is false 652 if(.false.)write(std_out,*)znucl ! Just to keep znucl when HAVE_LIBPSML is false 653 #endif 654 655 end subroutine psp9in