TABLE OF CONTENTS
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.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (JJ, MVer) 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
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
PARENTS
pspatm
CHILDREN
nctab_eval_tvalespl,pawrad_free,pawrad_init,ps_corecharge_get ps_destroy,ps_nonlocalprojectors_filter,ps_projector_get ps_provenance_get,ps_pseudoatomspec_get,ps_valenceconfiguration_get ps_valenceshell_get,psml_reader,psp8lo,psp8nl,psp9cc,spline,wrtout
SOURCE
76 #if defined HAVE_CONFIG_H 77 #include "config.h" 78 #endif 79 80 #include "abi_common.h" 81 82 subroutine psp9in(filpsp,ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,& 83 & mmax,mpsang,mpssoang,mqgrid,mqgrid_vl,nproj,n1xccc,pspso,qchrg,qgrid,qgrid_vl,& 84 & useylm,vlspl,xcccrc,xccc1d,zion,znucl,nctab,maxrad) 85 86 use defs_basis 87 use m_splines 88 use m_errors 89 use m_profiling_abi 90 91 use defs_datatypes, only : nctab_t 92 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 93 use m_psps, only : nctab_eval_tvalespl 94 95 #if defined HAVE_PSML 96 use m_psml 97 #endif 98 99 !This section has been created automatically by the script Abilint (TD). 100 !Do not modify the following lines by hand. 101 #undef ABI_FUNC 102 #define ABI_FUNC 'psp9in' 103 use interfaces_14_hidewrite 104 use interfaces_64_psp, except_this_one => psp9in 105 !End of the abilint section 106 107 implicit none 108 109 !Arguments ------------------------------------ 110 !scalars 111 integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mpsang,mpssoang,mqgrid,mqgrid_vl 112 integer,intent(in) :: pspso,n1xccc,useylm 113 integer,intent(out) :: mmax 114 real(dp),intent(in) :: zion,znucl 115 real(dp),intent(out) :: epsatm,qchrg,xcccrc,maxrad 116 type(nctab_t),intent(inout) :: nctab 117 character(len=fnlen),intent(in) :: filpsp 118 !arrays 119 integer,intent(out) :: indlmn(6,lmnmax),nproj(mpssoang) 120 real(dp),intent(in) :: qgrid(mqgrid),qgrid_vl(mqgrid_vl) 121 real(dp),intent(out) :: ekb(lnmax),ffspl(mqgrid,2,lnmax),vlspl(mqgrid,2) 122 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 123 124 !Local variables------------------------------- 125 !no_abirules 126 !scalars 127 #if defined HAVE_PSML 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_PSML 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_PSML 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 ! MSG_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 ! MSG_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_ALLOCATE( 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 MSG_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 MSG_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 MSG_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_ALLOCATE(rad,(mmax)) 369 ABI_ALLOCATE(vloc,(mmax)) 370 ABI_ALLOCATE(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(ab_out, '(a,i5,es16.6,es16.6)')' psp9in : mmax, amesh, rad(mmax) = ', mmax, amesh, rad(mmax) 383 384 !Check that rad grid is linear starting at zero 385 amesh=rad(2)-rad(1) 386 damesh=zero 387 do irad=2,mmax-1 388 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 389 end do 390 if(damesh>tol8 .or. rad(1)/=zero) then 391 write(message, '(5a)' )& 392 & 'Pseudopotential input file requires linear radial mesh',ch10,& 393 & 'starting at zero.',ch10,& 394 & 'Action: check your pseudopotential input file.' 395 MSG_ERROR(message) 396 end if 397 398 !Take care of the non-linear core corrections 399 !---------------------------------------------------------------------------- 400 ! xcccrc : XC core correction cutoff radius (bohr) 401 ! It is defined as the radius where the pseudo-core 402 ! charge density becomes zero 403 ! (here we have set up a tolerance of 1.d-12). 404 405 rmatch = zero 406 nders = 0 407 if (has_nlcc) then 408 409 ! In Abinit, at least for the Troullier-Martins pseudopotential, 410 ! the pseudocore charge density and its derivatives (xccc1d) 411 ! are introduced in a linear grid. 412 ! This grid is normalized, so the radial coordinates run between 413 ! from 0 and 1 (from 0 to xcccrc, where xcccrc is the radius 414 ! where the pseudo-core becomes zero). 415 416 call ps_CoreCharge_get(psxml, rc=rmatch, nderivs=nders) 417 write (message,'(1X,A,A,5X,A,1X,F8.3,A,5X,A,I8,A)') & 418 & "Reading pseudocore charge",ch10, & 419 & "- matching radius:",rmatch,ch10, & 420 & "- number of continuous derivatives",nders,ch10 421 call wrtout(std_out,message,'COLL') 422 423 !Get core charge function and derivatives, if needed 424 if(fchrg>1.0d-15)then 425 call psp9cc(psxml,mmax,n1xccc,rad,rchrg,xccc1d) 426 ! The core charge function for pspcod=9 427 ! becomes zero beyond rchrg. Thus xcccrc must be set 428 ! equal to rchrg. 429 xcccrc=rchrg 430 else 431 xccc1d(:,:) = zero 432 xcccrc = zero 433 fchrg = zero 434 qchrg = zero 435 end if 436 437 maxrad = rad(mmax) 438 439 end if ! has_nlcc 440 441 !! DEBUG 442 ! write(std_out,*)' xcccrc = ', xcccrc, rchrg 443 ! write(std_out,*) 444 ! write(std_out,*) '# psp8in NLCC data ', n1xccc, xcccrc 445 ! do ii = 1, n1xccc 446 ! write(std_out,'(7e20.8)')xcccrc*(ii-1.d0)/(n1xccc-1.d0),xccc1d(ii,1),& 447 ! & xccc1d(ii,2),xccc1d(ii,3),xccc1d(ii,4),xccc1d(ii,5),xccc1d(ii,6) 448 ! enddo 449 ! write(std_out,*) 450 ! stop 451 !! ENDDEBUG 452 453 454 !-------------------------------------------------------------------- 455 !Carry out calculations for local (lloc) pseudopotential. 456 !Obtain Fourier transform (1-d sine transform) 457 !to get q^2 V(q). 458 459 !Read and process vlocal: 460 !The local potential is given by a <radfunc> element under the <local-potential> 461 !element. 462 !After reading, this is a copy of the treatment to the 463 !local part carry out in psp8 464 !i.e. (as in Hamann pseudopotential) 465 ! 466 !Read the local component of the pseudopotential 467 vloc = zero 468 do ir = 1, mmax 469 vloc(ir) = ps_LocalPotential_Value(psxml, rad(ir)) 470 end do 471 472 call psp8lo(amesh,epsatm,mmax,mqgrid,qgrid,& 473 & vlspl(:,1),rad,vloc,yp1,ypn,zion) 474 475 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 476 ABI_ALLOCATE(work_spl,(mqgrid)) 477 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 478 vlspl(:,2)=work_spl(:) 479 ABI_DEALLOCATE(work_spl) 480 481 !! DEBUG 482 ! write(std_out,*)'# Vlocal = ' 483 ! write(std_out,*)' amesh = ', amesh 484 ! write(std_out,*)' epsatm = ', epsatm 485 ! write(std_out,*)' mmax = ', mmax 486 ! write(std_out,*)' mqgrid = ', mqgrid 487 ! do ir = 1, mqgrid 488 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 489 ! enddo 490 ! do ir = 1, mqgrid 491 ! write(std_out,'(a,i5,2f20.12)')' iq, vlspl = ', ir, vlspl(ir,1), vlspl(ir,2) 492 ! enddo 493 ! write(std_out,*) 494 ! do ir = 1, mmax 495 ! write(std_out,*)' rad = ', rad(ir), vloc(ir) 496 ! enddo 497 ! write(std_out,*) 498 ! write(std_out,*)' yp1 = ', yp1 499 ! write(std_out,*)' ypn = ', ypn 500 ! write(std_out,*)' zion = ', zion 501 ! stop 502 !! ENDDEBUG 503 504 505 !-------------------------------------------------------------------- 506 !Take care of non-local part 507 508 !Zero out all Kleinman-Bylander energies to initialize 509 do ii = 1, lmnmax ! loop over all possible projectors 510 if (indlmn(6,ii) == 1) then 511 call ps_Projector_Get(psxml, idx_sr(indlmn(5,ii)), ekb=ekb(indlmn(5,ii))) 512 else if (indlmn(6,ii) == 2) then 513 call ps_Projector_Get(psxml, idx_so(indlmn(5,ii)), ekb=ekb(indlmn(5,ii))) 514 end if 515 end do 516 517 !Read the KB projectors from the PSML file 518 !Note than in the PSML file the radial part of the projector is stored, 519 !while Abinit expects the radial part times the radii. 520 !We have to multiply by r after reading it. 521 !Note than in Hamann's format (psp8), Abinit directly reads r * radial_part_KB 522 vpspll = zero 523 do ii = 1, lmnmax 524 if (indlmn(6,ii) == 1) then 525 do ir = 1, mmax 526 vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_sr(indlmn(5,ii)), rad(ir)) 527 vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii)) 528 end do 529 else if (indlmn(6,ii) == 2) then 530 do ir = 1, mmax 531 vpspll(ir, indlmn(5,ii)) = ps_Projector_Value(psxml, idx_so(indlmn(5,ii)), rad(ir)) 532 vpspll(ir, indlmn(5,ii)) = rad(ir) * vpspll(ir, indlmn(5,ii)) 533 end do 534 end if 535 end do 536 537 !Allow for option of no nonlocal corrections (lloc=lmax=0) 538 if (lloc==0.and.lmax==0) then 539 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 540 call wrtout(ab_out,message,'COLL') 541 call wrtout(std_out,message,'COLL') 542 else 543 544 ! ---------------------------------------------------------------------- 545 ! Compute Vanderbilt-KB form factors and fit splines 546 547 call psp8nl(amesh,ffspl,indlmn,lmax,lmnmax,lnmax,mmax,& 548 & mqgrid,qgrid,rad,vpspll) 549 550 end if 551 552 !! DEBUG 553 ! write(std_out,*)'# KB Projectors = ' 554 ! write(std_out,*)' amesh = ', amesh 555 ! do ir = 1, mqgrid 556 ! do il = 1, lnmax 557 ! write(std_out,*)' iq, il, ffspl = ', ir, il, ffspl(ir,1,il), ffspl(ir,2,il) 558 ! enddo 559 ! enddo 560 ! do il = 1, lmnmax 561 ! write(std_out,*)' indlmn = ', il, indlmn(:,il) 562 ! enddo 563 ! write(std_out,*)' lmax = ', lmax 564 ! write(std_out,*)' lmnmax = ', lmnmax 565 ! write(std_out,*)' lnmax = ', lnmax 566 ! write(std_out,*)' mmax = ', mmax 567 ! write(std_out,*)' mqgrid = ', mqgrid 568 ! do ir = 1, mqgrid 569 ! write(std_out,*)' qgrid = ', ir, qgrid(ir) 570 ! enddo 571 ! do il = 1, lnmax 572 ! write(std_out,*) 573 ! write(std_out,*)'# il = ', il 574 ! do ir = 1, mmax 575 ! write(std_out,*)' rad = ', rad(ir), vpspll(ir,il) 576 ! enddo 577 ! enddo 578 ! stop 579 !! ENDDEBUG 580 581 ! Read pseudo valence charge in real space on the linear mesh 582 ! and transform it to reciprocal space on a regular grid. Use vloc as workspace. 583 vloc(:) = zero 584 if (has_tvale) then 585 do irad=1,mmax 586 vloc(irad) = ps_ValenceCharge_Value(psxml,rad(irad)) 587 vloc(irad) = vloc(irad) / four_pi 588 end do 589 590 !! DEBUG 591 ! do irad = 1, mmax 592 ! write(std_out,*)' Valence Charge = ', rad(irad), vloc(irad) 593 ! enddo 594 ! stop 595 !! ENDDEBUG 596 597 598 ! Check that rad grid is linear starting at zero 599 amesh=rad(2)-rad(1) 600 damesh=zero 601 do irad=2,mmax-1 602 damesh=max(damesh,abs(rad(irad)+amesh-rad(irad+1))) 603 end do 604 if(damesh>tol8 .or. rad(1)/=zero) then 605 write(message, '(5a)' )& 606 & 'Pseudopotential input file requires linear radial mesh',ch10,& 607 & 'starting at zero.',ch10,& 608 & 'Action: check your pseudopotential input file.' 609 MSG_ERROR(message) 610 end if 611 612 ! Evaluate spline-fit of the atomic pseudo valence charge in reciprocal space. 613 call pawrad_init(mesh,mesh_size=mmax,mesh_type=1,rstep=amesh) 614 call nctab_eval_tvalespl(nctab, zion, mesh, vloc, mqgrid_vl, qgrid_vl) 615 call pawrad_free(mesh) 616 end if 617 618 ABI_DEALLOCATE(vpspll) 619 ABI_DEALLOCATE(vloc) 620 ABI_DEALLOCATE(rad) 621 if (allocated(idx_sr)) then 622 ABI_FREE_NOCOUNT(idx_sr) 623 end if 624 if (allocated(idx_so)) then 625 ABI_FREE_NOCOUNT(idx_so) 626 end if 627 628 call ps_destroy(psxml) 629 630 !-------------------------------------------------------------------- 631 632 #else 633 ABI_UNUSED(mpsang) 634 ABI_UNUSED(pspso) 635 ABI_UNUSED(qgrid_vl) 636 ABI_UNUSED(nctab%mqgrid_vl) 637 !Initialize some arguments, for portability at compile time 638 indlmn=0 ; mmax=0 ; nproj=0 639 ekb=zero ; epsatm=zero ; ffspl=zero ; qchrg=zero ; vlspl=zero ; xcccrc=zero ; xccc1d=zero 640 641 if(.false.)write(std_out,*)filpsp ! Just to keep filpsp when HAVE_PSML is false 642 if(.false.)write(std_out,*)lloc ! Just to keep lloc when HAVE_PSML is false 643 if(.false.)write(std_out,*)lmax ! Just to keep lmax when HAVE_PSML is false 644 if(.false.)write(std_out,*)mpsang ! Just to keep mpsang when HAVE_PSML is false 645 if(.false.)write(std_out,*)pspso ! Just to keep pspso when HAVE_PSML is false 646 if(.false.)write(std_out,*)qgrid ! Just to keep qgrid when HAVE_PSML is false 647 if(.false.)write(std_out,*)qgrid_vl ! Just to keep qgrid_vl when HAVE_PSML is false 648 if(.false.)write(std_out,*)useylm ! Just to keep useylm when HAVE_PSML is false 649 if(.false.)write(std_out,*)zion ! Just to keep zion when HAVE_PSML is false 650 if(.false.)write(std_out,*)znucl ! Just to keep znucl when HAVE_PSML is false 651 #endif 652 653 end subroutine psp9in