TABLE OF CONTENTS
ABINIT/psp8in [ Functions ]
NAME
psp8in
FUNCTION
Initialize pspcod=8 (pseudopotentials in the format generated by DRH): continue to read the corresponding file, then compute the local and non-local potentials.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (DRH, XG, AF) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
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, as read from input file 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
PARENTS
pspatm
CHILDREN
nctab_eval_tvalespl,pawrad_free,pawrad_init,psp8cc,psp8lo,psp8nl,spline wrtout
SOURCE
68 #if defined HAVE_CONFIG_H 69 #include "config.h" 70 #endif 71 72 #include "abi_common.h" 73 74 subroutine psp8in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 75 & mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,& 76 & useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad) 77 78 use defs_basis 79 use m_splines 80 use m_errors 81 use m_profiling_abi 82 83 use defs_datatypes, only : nctab_t 84 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 85 use m_psps, only : nctab_eval_tvalespl 86 87 !This section has been created automatically by the script Abilint (TD). 88 !Do not modify the following lines by hand. 89 #undef ABI_FUNC 90 #define ABI_FUNC 'psp8in' 91 use interfaces_14_hidewrite 92 use interfaces_64_psp, except_this_one => psp8in 93 !End of the abilint section 94 95 implicit none 96 97 !Arguments ------------------------------------ 98 !scalars 99 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid,mqgrid_vl 100 integer,intent(in) :: pspso,n1xccc,useylm 101 real(dp),intent(in) :: zion,znucl 102 real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad 103 type(nctab_t),intent(inout) :: nctab 104 !arrays 105 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang) 106 real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl) 107 real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2) 108 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 109 110 !Local variables------------------------------- 111 !scalars 112 integer :: extension_switch,iln,iln0,pspindex,ipsang,irad,jj,kk,ll,ll_err,llin 113 integer :: mm,nn,nso 114 real(dp) :: amesh,damesh,fchrg,rchrg,yp1,ypn 115 logical :: has_tvale 116 character(len=500) :: message,errmsg 117 type(pawrad_type) :: mesh 118 !arrays 119 integer, allocatable :: nproj_tmp(:) 120 real(dp),allocatable :: rad(:),vloc(:),vpspll(:,:),work_spl(:) 121 122 ! *************************************************************************** 123 124 !DEBUG 125 !write(std_out,*)' psp8in : enter and stop ' 126 !stop 127 !ENDDEBUG 128 129 130 !File format of formatted drh psp input, as adapted for use 131 !by the ABINIT code (the 3 first lines 132 !have already been read in calling -pspatm- routine) : 133 134 !(1) title (character) line 135 !(2) znucl,zion,pspdat 136 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well (r2well not used) 137 !(4) rchrg,fchrg,qchrg (fchrg /=0 if core charge, qchrg not used) 138 !(5) nproj(0:lmax) (several projectors allowed for each l) 139 !(6) extension_switch(2) (spin-orbit parameters) 140 !Then, for ll=0,lmax : 141 !if(nproj(ll)>0) 142 !1/<u1|vbkb1>, 1/<u2|vbkb2>, ... 143 !for irad=1,mmax : irad, r(irad), vbkb1(irad,ll), vbkb2(irad,ll), ... 144 !elseif ll=lloc 145 !for irad=1,mmax : irad, r(irad), vloc(irad) 146 !end if 147 ! 148 !If(lloc>lmax) 149 !for irad=1,mmax : irad, r(irad), vloc(irad) 150 !end if 151 ! 152 !vbkb are Bloechl-Kleinman-Bylander projectors,(vpsp(r,ll)-vloc(r))*u(r,ll), 153 !unnormalized 154 !Note that an arbitrary local potential is allowed. Set lloc>lmax, and 155 !provide projectors for all ll<=lmax 156 ! 157 !Finally, if fchrg>0, 158 !for irad=1,mmax : irad, r(irad), xccc(irad), 159 !xccc'(irac), xccc''(irad), xccc'''(irad), xccc''''(irad) 160 ! 161 !Model core charge for nonlinear core xc correction, and 4 derivatives 162 163 read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg 164 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 165 call wrtout(ab_out,message,'COLL') 166 call wrtout(std_out, message,'COLL') 167 168 ABI_ALLOCATE(nproj_tmp,(mpssoang)) 169 170 nproj_tmp(:)=0 171 read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(1:lmax+1) 172 write(message, '(a,5i6)' ) ' nproj',nproj_tmp(1:lmax+1) 173 call wrtout(ab_out,message,'COLL') 174 call wrtout(std_out,message,'COLL') 175 176 !place holder for future implementation of additional optional header 177 !lines without invalidating existing psp files 178 !Now (12/2014) extended to include spin-orbit projectors 179 180 ! The integer labeled "extension switch" on line 6 181 ! of the *.psp8 file will be set to 1 (non- or scalar-relativistic) 182 ! or 3 (relativistic) to signal this to Abinit that the file contains the pseudo valence charge. 183 184 has_tvale = .False. 185 read (tmp_unit,*, err=10, iomsg=errmsg) extension_switch 186 if (any(extension_switch==[2, 3])) then 187 read (tmp_unit,*, err=10, iomsg=errmsg) nproj_tmp(lmax+2:2*lmax+1) 188 write(message, '(5x,a,i6)' ) 'spin-orbit psp, extension_switch',extension_switch 189 call wrtout(ab_out,message,'COLL') 190 call wrtout(std_out, message,'COLL') 191 write(message, '(5x,a,5i6)' ) ' nprojso',nproj_tmp(lmax+2:2*lmax+1) 192 call wrtout(ab_out,message,'COLL') 193 call wrtout(std_out, message,'COLL') 194 has_tvale = (extension_switch == 3) 195 else if (any(extension_switch==[0,1])) then 196 write(message, '(5x,a,i6)' ) 'extension_switch',extension_switch 197 call wrtout(ab_out,message,'COLL') 198 call wrtout(std_out, message,'COLL') 199 has_tvale = (extension_switch == 1) 200 else 201 write(message, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,& 202 & 'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit' 203 MSG_ERROR(message) 204 end if 205 206 if(lloc<4) then 207 if (nproj_tmp(lloc+1)/=0) then 208 write(message, '(a,i4,a,a,i4,5a)' )& 209 & 'Pseudopotential input file has nproj=',nproj_tmp(lloc+1),ch10,& 210 & 'for angular momentum',lloc,' which is the local potential.',ch10,& 211 & 'Should be 0 for the local potential',ch10,& 212 & 'Action: check your pseudopotential input file.' 213 MSG_ERROR(message) 214 end if 215 end if 216 217 !-------------------------------------------------------------------- 218 219 !Initialize array indlmn giving l,m,n,lm,ln,s for i=lmn 220 ! if(pspso==2) then 221 if (any(extension_switch == [0,1])) then 222 nso=1 223 else if (any(extension_switch == [2,3])) then 224 nso=2 225 if (pspso==0) then 226 write (message, '(3a)') 'You are reading a pseudopotential file with spin orbit projectors',ch10,& 227 & ' but internal variable pspso is 0 - this is not possible for oncvpsp. Set nspinor to 2 and so_psp 1' 228 MSG_COMMENT(message) 229 end if 230 else 231 write(message, '(a,i0,2a)' ) 'invalid extension_switch: ',extension_switch,ch10,& 232 & 'Should be [0,1] for scalar-relativistic psp or [2,3] to include spin-orbit' 233 MSG_ERROR(message) 234 end if 235 236 pspindex=0;iln=0;indlmn(:,:)=0 237 do nn=1,nso 238 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 239 ll=ipsang-(nn-1)*lmax-1 240 if (nproj_tmp(ipsang)>0) then 241 do kk=1,nproj_tmp(ipsang) 242 iln=iln+1 243 do mm=1,2*ll*useylm+1 244 pspindex=pspindex+1 245 indlmn(1,pspindex)=ll 246 indlmn(2,pspindex)=mm-ll*useylm-1 247 indlmn(3,pspindex)=kk 248 indlmn(4,pspindex)=ll*ll+(1-useylm)*ll+mm 249 indlmn(5,pspindex)=iln 250 indlmn(6,pspindex)=nn 251 end do 252 end do 253 end if 254 end do 255 end do 256 257 ! repackage nproj_tmp for proper use by pspatm 258 nproj(:)=0 259 nproj(1:lmax+1)=nproj_tmp(1:lmax+1) 260 if(pspso==2) then 261 nproj(mpsang+1:mpsang+lmax)=nproj_tmp(lmax+2:2*lmax+1) 262 end if 263 264 !Can now allocate grids, potentials and projectors 265 ABI_ALLOCATE(rad,(mmax)) 266 ABI_ALLOCATE(vloc,(mmax)) 267 ABI_ALLOCATE(vpspll,(mmax,lnmax)) 268 269 !Will now proceed at the reading of pots and projectors 270 271 !rad(:)=radial grid r(i) 272 !vpspll(:,1),...,vpspll(:,lnmax)=nonlocal projectors 273 !vloc(:)=local potential 274 275 !Read Vanderbilt-Kleinman-Bylander energies and projectors for each l 276 !or read local potential for l=lloc. 277 !Also get rad array (actually read more than once) 278 ll_err=0 279 iln0=0 280 do nn=1,nso 281 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 282 ll=ipsang-(nn-1)*lmax-1 283 if (nproj_tmp(ipsang)>0) then 284 read(tmp_unit,*, err=10, iomsg=errmsg) llin,ekb(iln0+1:iln0+nproj_tmp(ipsang)) 285 if(llin/=ll) then 286 ll_err=ipsang 287 exit 288 end if 289 do irad=1,mmax 290 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vpspll(irad,iln0+1:iln0+nproj_tmp(ipsang)) 291 end do 292 iln0=iln0+nproj_tmp(ipsang) 293 elseif(ll==lloc .and. nn==1) then 294 read(tmp_unit,*, err=10, iomsg=errmsg) llin 295 if(llin/=ll) then 296 ll_err=ipsang 297 exit 298 end if 299 do irad=1,mmax 300 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad) 301 end do 302 end if 303 end do !ipsang 304 !Provision for general local potential /= any angular momentum potential 305 if(nn==1 .and. lloc>lmax) then 306 read(tmp_unit,*, err=10, iomsg=errmsg) llin 307 if(llin==lloc) then 308 do irad=1,mmax 309 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad) 310 end do 311 else 312 ll_err=lloc+1 313 exit 314 end if 315 end if 316 end do !nn 317 318 if(ll_err>0) then 319 write(message, '(5a,i4,a,i4,a,a)' )& 320 & 'Pseudopotential input file does not have angular momenta in order',ch10,& 321 & 'or has inconsistent general local potential index',ch10,& 322 & 'Expected',ll_err-1,' , got',ll,ch10,& 323 & 'Action: check your pseudopotential input file.' 324 MSG_ERROR(message) 325 end if 326 327 !Check that rad grid is linear starting at zero 328 amesh=rad(2)-rad(1) 329 damesh=zero 330 do irad=2,mmax-1 331 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 332 end do 333 if(damesh>tol8 .or. rad(1)/=zero) then 334 write(message, '(5a)' )& 335 & 'Pseudopotential input file requires linear radial mesh',ch10,& 336 & 'starting at zero.',ch10,& 337 & 'Action: check your pseudopotential input file.' 338 MSG_ERROR(message) 339 end if 340 341 !Get core charge function and derivatives, if needed 342 if(fchrg>1.0d-15)then 343 call psp8cc(mmax,n1xccc,rchrg,xccc1d) 344 ! The core charge function for pspcod=8 345 ! becomes zero beyond rchrg. Thus xcccrc must be set 346 ! equal to rchrg. 347 xcccrc=rchrg 348 else 349 xccc1d(:,:) = zero 350 xcccrc = zero 351 fchrg = zero 352 qchrg = zero 353 end if 354 355 maxrad = rad(mmax) 356 357 !! DEBUG 358 ! write(std_out,*)' xcccrc = ', xcccrc, rchrg 359 ! write(std_out,*) 360 ! write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc 361 ! do ii = 1, n1xccc 362 ! write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),& 363 ! & xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6) 364 ! enddo 365 ! write(std_out,*) 366 ! stop 367 !! ENDDEBUG 368 369 370 !-------------------------------------------------------------------- 371 !Carry out calculations for local (lloc) pseudopotential. 372 !Obtain Fourier transform (1-d sine transform) 373 !to get q^2 V(q). 374 375 call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,& 376 & vlspl(:,1),rad,vloc,yp1,ypn,zion) 377 378 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 379 ABI_ALLOCATE(work_spl,(mqgrid)) 380 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 381 vlspl(:,2)=work_spl(:) 382 ABI_DEALLOCATE(work_spl) 383 384 !! DEBUG 385 ! write(std_out,*)'# Vlocal = ' 386 ! write(std_out,*)' amesh = ', amesh 387 ! write(std_out,*)' epsatm = ', epsatm 388 ! write(std_out,*)' mmax = ', mmax 389 ! write(std_out,*)' mqgrid = ', mqgrid 390 ! do ir = 1, mqgrid 391 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 392 ! enddo 393 ! do ir = 1, mqgrid 394 ! write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2) 395 ! enddo 396 ! write(std_out,*) 397 ! do ir = 1, mmax 398 ! write(std_out,*)' rad = ', rad(ir), vloc(ir) 399 ! enddo 400 ! write(std_out,*) 401 ! write(std_out,*)' yp1 = ', yp1 402 ! write(std_out,*)' ypn = ', ypn 403 ! write(std_out,*)' zion = ', zion 404 ! stop 405 !! ENDDEBUG 406 407 408 !-------------------------------------------------------------------- 409 !Take care of non-local part 410 411 !Allow for option of no nonlocal corrections (lloc=lmax=0) 412 if (lloc==0.and.lmax==0) then 413 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 414 call wrtout(ab_out,message,'COLL') 415 call wrtout(std_out,message,'COLL') 416 else 417 418 ! ---------------------------------------------------------------------- 419 ! Compute Vanderbilt-KB form factors and fit splines 420 421 call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,& 422 & mqgrid,qgrid,rad,vpspll) 423 424 end if 425 426 !! DEBUG 427 ! write(std_out,*)'# KB Projectors = ' 428 ! write(std_out,*)' amesh = ', amesh 429 ! do ir = 1, mqgrid 430 ! do il = 1, lnmax 431 ! write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il) 432 ! enddo 433 ! enddo 434 ! do il = 1, lmnmax 435 ! write(std_out,*)' indlmn = ', il, indlmn(:,il) 436 ! enddo 437 ! write(std_out,*)' lmax = ', lmax 438 ! write(std_out,*)' lmnmax = ', lmnmax 439 ! write(std_out,*)' lnmax = ', lnmax 440 ! write(std_out,*)' mmax = ', mmax 441 ! write(std_out,*)' mqgrid = ', mqgrid 442 ! do ir = 1, mqgrid 443 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 444 ! enddo 445 ! do il = 1, lnmax 446 ! write(std_out,*) 447 ! write(std_out,*)'# il = ', il 448 ! do ir = 1, mmax 449 ! write(std_out,*)' rad = ', rad(ir), vpspll(ir,il) 450 ! enddo 451 ! enddo 452 ! stop 453 !! ENDDEBUG 454 455 ! Read pseudo valence charge in real space on the linear mesh 456 ! and transform it to reciprocal space on a regular grid. Use vloc as workspace. 457 vloc(:) = zero 458 if (has_tvale) then 459 do irad=1,mmax 460 read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),vloc(irad) 461 vloc(irad) = vloc(irad) / four_pi 462 end do 463 464 !! DEBUG 465 ! do irad = 1, mmax 466 ! write(std_out,*)' Valence Charge = ', rad(irad), vloc(irad) 467 ! enddo 468 ! stop 469 !! ENDDEBUG 470 471 472 ! Check that rad grid is linear starting at zero 473 amesh=rad(2)-rad(1) 474 damesh=zero 475 do irad=2,mmax-1 476 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 477 end do 478 if(damesh>tol8 .or. rad(1)/=zero) then 479 write(message, '(5a)' )& 480 & 'Pseudopotential input file requires linear radial mesh',ch10,& 481 & 'starting at zero.',ch10,& 482 & 'Action: check your pseudopotential input file.' 483 MSG_ERROR(message) 484 end if 485 486 ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space. 487 call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh) 488 call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl) 489 call pawrad_free(mesh) 490 end if 491 492 ABI_DEALLOCATE(vpspll) 493 ABI_DEALLOCATE(vloc) 494 ABI_DEALLOCATE(rad) 495 ABI_DEALLOCATE(nproj_tmp) 496 497 return 498 499 ! Handle IO error 500 10 continue 501 MSG_ERROR(errmsg) 502 503 end subroutine psp8in