TABLE OF CONTENTS
ABINIT/psp5in [ Functions ]
NAME
psp5in
FUNCTION
Initialize pspcod=5 ("Phoney pseudopotentials" with Hamman grid):
continue to read the corresponding file, then compute the
local and non-local potentials.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, FrD, FJ, MT) 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
mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
mpssoang= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
mqgrid=dimension of q (or G) grid for arrays.
n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
pspso= spin orbit signal
qgrid(mqgrid)=values of q (or |G|) on grid from 0 to qmax
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) if any, spin-orbit components begin at l=mpsang+1 ekb1(mpssoang)= Kleinman-Bylander energy from the psp file, for iproj=1 ekb2(mpssoang)= Kleinman-Bylander energy from the psp file, for iproj=2 epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r} dr]$ (hartree) epspsp(mpssoang)=values of epsatm for different angular momenta, from the psp file e990(mpssoang)=ecut at which 0.99 of the kinetic energy is recovered e999(mpssoang)=ecut at which 0.999 of the kinetic energy is recovered ffspl(mqgrid,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 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 the total (integrated) core charge rcpsp(mpssoang)=cut-off radius for each angular momentum rms(mpssoang)=root mean square of the KB psp vlspl(mqgrid,2)=q^2 Vloc(q) and second derivatives from spline fit xcccrc=XC core correction cutoff radius (bohr) from psp file xccc1d(n1xccc,6)=1D core charge function and five derivatives, from psp file
PARENTS
pspatm
CHILDREN
psp1cc,psp5lo,psp5nl,spline,wrtout
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 subroutine psp5in(ekb,ekb1,ekb2,epsatm,epspsp,e990,e999,ffspl,indlmn,& 78 & lloc,lmax,lmnmax,lnmax,mmax,mpsang,mpssoang,mqgrid,& 79 & nproj,n1xccc,pspso,qchrg,qgrid,rcpsp,rms,& 80 & useylm,vlspl,xcccrc,xccc1d,zion,znucl) 81 82 use defs_basis 83 use m_splines 84 use m_errors 85 use m_profiling_abi 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 'psp5in' 91 use interfaces_14_hidewrite 92 use interfaces_64_psp, except_this_one => psp5in 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 100 integer,intent(in) :: n1xccc,pspso,useylm 101 real(dp),intent(in) :: zion,znucl 102 real(dp),intent(out) :: epsatm,qchrg,xcccrc 103 !arrays 104 integer,intent(out) :: indlmn(6,lmnmax) !vz_i 105 integer,intent(inout) :: nproj(mpssoang) !vz_i 106 real(dp),intent(in) :: qgrid(mqgrid) 107 real(dp),intent(out) :: e990(mpssoang),e999(mpssoang),ekb(lnmax) 108 real(dp),intent(out) :: ekb1(mpssoang),ekb2(mpssoang),epspsp(mpssoang) 109 real(dp),intent(out) :: rcpsp(mpssoang),rms(mpssoang) !vz_i 110 real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i 111 real(dp),intent(out) :: vlspl(mqgrid,2) !vz_i 112 real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i 113 114 !Local variables------------------------------- 115 !scalars 116 integer :: i1,i2,ii,iln,index,ipsang,kk,lhigh,ll,mm,mproj,nn,nso,pspso0 117 real(dp) :: al,fchrg,r1,rchrg,yp1,ypn 118 logical :: test 119 character(len=500) :: message,errmsg 120 !arrays 121 real(dp),allocatable :: ekb_so(:),ekb_sr(:),ekb_tmp(:,:),ffspl_so(:,:,:) 122 real(dp),allocatable :: ffspl_sr(:,:,:),ffspl_tmp(:,:,:,:),rad(:),vloc(:) 123 real(dp),allocatable :: vpspll(:,:),vpspll_so(:,:),wfll(:,:),wfll_so(:,:) 124 real(dp),allocatable :: work_space(:),work_spl(:) 125 126 ! *************************************************************************** 127 128 !DEBUG 129 !write(std_out,*)' psp5in : enter ' 130 !ENDDEBUG 131 132 !File format of formatted Phoney psp input (the 3 first lines 133 !have already been read in calling -pspatm- routine) : 134 135 !(1) title (character) line 136 !(2) znucl,zion,pspdat 137 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well 138 !(4) r1,al,pspso 139 !For each angular momentum : 140 !(4) ll,e990(ll),e999(ll),nproj(ll),rcpsp(ll) 141 !(5) rms(ll),ekb1(ll),ekb2(ll),epspsp(ll) 142 !(6) rchrg,fchrg,qchrg 143 !(7) ll 144 !(8) (vpsp(j,ll),j=0,nmax) 145 !Then for iproj=1 to 2 146 !for ll=0,lmax 147 !(10) ll 148 !(11) ((upsp(j,ll,iproj),j=0,nmax) 149 150 !Read fourth line of the file ; parameter pspso is optional (but 151 !this is not treated correctly by all machines - problems with SGI ) 152 pspso0=1 153 read (tmp_unit,fmt=*,err=50,end=50) r1,al,pspso0 154 50 continue 155 156 if(pspso0/=1 .and. pspso0/=2)then 157 write(message, '(3a,i0,2a)' )& 158 & 'Problem reading the fourth line of pseudopotential file.',ch10,& 159 & 'The parameter pspso should be 1 or 2, but it is pspso= ',pspso0,ch10,& 160 & 'Action: check your pseudopotential input file.' 161 MSG_ERROR(message) 162 end if 163 164 write(message, '(2es16.6,t47,a)' ) r1,al,'r1 and al (Hamman grid)' 165 call wrtout(ab_out,message,'COLL') 166 call wrtout(std_out, message,'COLL') 167 168 if (pspso0/=1) then 169 write(message,'(a)') ' Pseudopotential is in spin-orbit format ' 170 call wrtout(ab_out,message,'COLL') 171 call wrtout(std_out, message,'COLL') 172 end if 173 174 if (pspso/=0.and.pspso0==1) then 175 write(message, '(a,a,a,a,a)' )& 176 & 'The treatment of spin-orbit interaction is required (pspso/=0)',ch10,& 177 & 'but pseudopotential file format cannot contain spin-orbit information !',ch10,& 178 & 'Action : check your pseudopotential input file.' 179 MSG_ERROR(message) 180 end if 181 182 nso=1;if (pspso0/=1) nso=2 183 do ipsang=1,(nso*lmax)+1 184 read (tmp_unit,*, err=10, iomsg=errmsg) ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang) 185 write(message, '(i5,2f8.3,i5,f12.7,t47,a)' ) & 186 & ll,e990(ipsang),e999(ipsang),nproj(ipsang),rcpsp(ipsang),'l,e99.0,e99.9,nproj,rcpsp' 187 call wrtout(ab_out,message,'COLL') 188 call wrtout(std_out, message,'COLL') 189 read (tmp_unit,*, err=10, iomsg=errmsg) rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang) 190 write(message, '(4f13.8,t55,a)' ) & 191 & rms(ipsang),ekb1(ipsang),ekb2(ipsang),epspsp(ipsang),' rms, ekb1, ekb2, epsatm' 192 call wrtout(ab_out,message,'COLL') 193 call wrtout(std_out, message,'COLL') 194 end do 195 196 !If pspso/=0 and nproj/=2, forces nproj to be 2 197 !(for compatibility with old psp-file format) 198 if (pspso/=0.and.lmax>0) then 199 test=.false. 200 do ipsang=1,(nso*lmax)+1 201 if (ipsang>1.and.nproj(ipsang)/=2) then 202 test=.true.;nproj(ipsang)=2 203 end if 204 end do 205 if (test) then 206 write(message, '(a,a,a,a,a)' )& 207 & 'Pseudopotential file is spin-orbit (pspso=2)',ch10,& 208 & 'and number of projector for l/=0 is not 2 !',ch10,& 209 & 'It has been forced to 2.' 210 call wrtout(std_out,message,'COLL') 211 MSG_WARNING(message) 212 end if 213 end if 214 215 !mproj=maxval(nproj(1:lmax+1)) 216 !mjv 10/2008: I believe this is correct. Perhaps unnecessary if the normal 217 !projectors are always more numerous, but should be conservative anyway with 218 !maxval. 219 mproj=maxval(nproj) 220 index=0;iln=0;indlmn(:,:)=0 221 do nn=1,nso 222 do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1 223 if (nproj(ipsang)>0) then 224 ll=ipsang-(nn-1)*lmax-1 225 do kk=1,nproj(ipsang) 226 iln=iln+1 227 do mm=1,2*ll*useylm+1 228 index=index+1 229 indlmn(1,index)=ll 230 indlmn(2,index)=mm-ll*useylm-1 231 indlmn(3,index)=kk 232 indlmn(4,index)=ll*ll+(1-useylm)*ll+mm 233 indlmn(5,index)=iln 234 indlmn(6,index)=nn 235 end do 236 end do 237 end if 238 end do 239 end do 240 241 read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg 242 write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg' 243 call wrtout(ab_out,message,'COLL') 244 call wrtout(std_out, message,'COLL') 245 246 !Generate core charge function and derivatives, if needed 247 xcccrc=zero 248 if(n1xccc>0)then 249 ! Use the revised expression of 5 Nov 1992, also used for format=1. 250 call psp1cc(fchrg,n1xccc,xccc1d) 251 xcccrc=3*rchrg 252 end if 253 254 !-------------------------------------------------------------------- 255 !Will now proceed at the reading of pots and wfs, as well as their treatment 256 257 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials 258 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg. 259 !rad(:)=radial grid r(i) 260 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions 261 ABI_ALLOCATE(vloc,(mmax)) 262 ABI_ALLOCATE(vpspll,(mmax,mpsang)) 263 264 !(1) Read atomic pseudopotential for each l, filling up array vpspll 265 !Note: put each l into vpspll(:,l+1) 266 267 if (pspso==0) then 268 269 ! --NON SPIN-ORBIT 270 do ipsang=1,lmax+1 271 read (tmp_unit,*, err=10, iomsg=errmsg) ll 272 read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax) 273 ! write(std_out,*) 'END OF READING PSP',ll,'OK' 274 end do 275 else 276 277 ! --SPIN-ORBIT 278 ABI_ALLOCATE(vpspll_so,(mmax,mpsang)) 279 read (tmp_unit,*, err=10, iomsg=errmsg) ll 280 read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll(ii,1),ii=1,mmax) 281 vpspll_so(:,1)=0.0d0 282 do ipsang=2,lmax+1 283 read (tmp_unit,*, err=10, iomsg=errmsg) ll 284 read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll(ii,ipsang),ii=1,mmax) 285 read (tmp_unit,*, err=10, iomsg=errmsg) ll 286 read (tmp_unit,*, err=10, iomsg=errmsg) (vpspll_so(ii,ipsang),ii=1,mmax) 287 end do 288 end if 289 290 !Copy appropriate nonlocal psp for use as local one 291 if (pspso==0) then 292 vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 ) 293 else 294 if(lloc<=0) then 295 vloc( 1:mmax ) = vpspll( 1:mmax , -lloc+1 ) 296 else 297 vloc( 1:mmax ) = vpspll_so( 1:mmax , lloc+1 ) 298 end if 299 end if 300 !DEBUG 301 !write(std_out,*) 'VLOC=',vloc(1),vloc(2),vloc(3) 302 !write(std_out,*) 'VLOC=',vloc(4),vloc(5),vloc(6) 303 !ENDDEBUG 304 305 306 !(2) Create radial grid, and associated quantities 307 308 !Now compute Hamman Grid 309 ABI_ALLOCATE(rad,(mmax)) 310 do ii=1,mmax 311 rad (ii)=r1*exp(dble(ii-1)*al) 312 end do 313 !DEBUG 314 !write(std_out,*) 'HAMMAN RADIAL GRID r1 and al',r1,al 315 !write(std_out,*) 'rad(1)=',rad(1) 316 !write(std_out,*) 'rad(10)=',rad(10) 317 !write(std_out,*) 'rad(100)=',rad(100) 318 !ENDDEBUG 319 320 321 !(3)Carry out calculations for local (lloc) pseudopotential. 322 !Obtain Fourier transform (1-d sine transform) 323 !to get q^2 V(q). 324 325 call psp5lo(al,epsatm,mmax,mqgrid,qgrid,& 326 & vlspl(:,1),rad,vloc,yp1,ypn,zion) 327 328 329 !Fit spline to q^2 V(q) (Numerical Recipes subroutine) 330 ABI_ALLOCATE(work_space,(mqgrid)) 331 ABI_ALLOCATE(work_spl,(mqgrid)) 332 call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl) 333 vlspl(:,2)=work_spl(:) 334 335 ABI_DEALLOCATE(work_space) 336 ABI_DEALLOCATE(work_spl) 337 338 !(4)Take care of non-local part 339 340 !DEBUG 341 !write(std_out,*)' psp5in : before nonlocal corrections ' 342 !write(std_out,*)' psp5in : lloc, lmax = ',lloc,lmax 343 !ENDDEBUG 344 345 !Zero out all Kleinman-Bylander energies to initialize 346 ekb(:)=0.0d0 347 348 !Allow for option of no nonlocal corrections (lloc=lmax=0) 349 if (lloc==0.and.lmax==0) then 350 351 write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl 352 call wrtout(ab_out,message,'COLL') 353 call wrtout(std_out, message,'COLL') 354 355 else 356 357 ! Proceed to make Kleinman-Bylander form factors for 358 ! each l up to lmax 359 360 ! Read wavefunctions for each l up to lmax 361 ABI_ALLOCATE( wfll,(mmax,mpsang)) 362 ! ----------------------------------------------------------------- 363 364 if (pspso==0) then 365 366 ! --NON SPIN-ORBIT 367 do ipsang=1,lmax+1 368 if (nproj(ipsang)/=0) then 369 read (tmp_unit,*, err=10, iomsg=errmsg) ll 370 if (ipsang/=ll+1) then 371 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 372 & 'Pseudopotential input file does not have',ch10,& 373 & 'angular momenta in order expected for first projection',& 374 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 375 & 'Action : check your pseudopotential input file.' 376 MSG_ERROR(message) 377 end if 378 read (tmp_unit,*, err=10, iomsg=errmsg) wfll(:,ipsang) 379 else 380 wfll(:,ipsang)=0.0d0 381 end if 382 end do 383 else 384 385 ! --SPIN-ORBIT 386 ABI_ALLOCATE(wfll_so,(mmax,mpsang)) 387 if (nproj(1)/=0) then 388 read (tmp_unit,*,err=10,iomsg=errmsg) ll 389 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,1) 390 else 391 wfll(:,1)=0.0d0 392 end if 393 wfll_so(:,1)=0.0d0 394 do ipsang=2,lmax+1 395 if (nproj(ipsang)/=0) then 396 read (tmp_unit,*,err=10,iomsg=errmsg) ll 397 read (tmp_unit,*,err=10,iomsg=errmsg) wfll(:,ipsang) 398 read (tmp_unit,*,err=10,iomsg=errmsg) ll 399 read (tmp_unit,*,err=10,iomsg=errmsg) wfll_so(:,ipsang) 400 else 401 wfll(:,ipsang)=0.0d0 402 wfll_so(:,ipsang)=0.0d0 403 end if 404 end do 405 end if 406 407 ! ---------------------------------------------------------------------- 408 ! Compute KB form factors and fit splines 409 ABI_ALLOCATE(ekb_tmp,(mpssoang,max(nso,mproj))) 410 ABI_ALLOCATE(ffspl_tmp,(mqgrid,2,mpssoang,max(nso,mproj))) 411 ekb_tmp(:,:)=0.d0 412 413 ABI_ALLOCATE(ekb_sr,(mpsang)) 414 ABI_ALLOCATE(ffspl_sr,(mqgrid,2,mpsang)) 415 call psp5nl(al,ekb_sr(:),ffspl_sr(:,:,:),lmax,mmax,mpsang,mqgrid,& 416 & qgrid,rad,vloc,vpspll,wfll) 417 ekb_tmp(1:mpsang,1)=ekb_sr(1:mpsang) 418 ffspl_tmp(:,:,1:mpsang,1)=ffspl_sr(:,:,1:mpsang) 419 420 if (pspso/=0) then 421 ABI_ALLOCATE(ekb_so,(mpsang)) 422 ABI_ALLOCATE(ffspl_so,(mqgrid,2,mpsang)) 423 call psp5nl(al,ekb_so,ffspl_so,lmax,mmax,mpsang,mqgrid,& 424 & qgrid,rad,vloc,vpspll_so,wfll_so) 425 ekb_tmp(mpsang+1:mpssoang,1)=ekb_so(2:mpsang) 426 do ipsang=2,lmax+1 427 if((ekb_sr(ipsang)*ekb_so(ipsang))<0.0) then 428 MSG_ERROR('BIG PROBLEM WITH THE SPIN ORBIT IN PSP5NL') 429 end if 430 end do 431 432 if(lloc<0) ekb_sr(-lloc+1)=ekb_so(-lloc+1) 433 if(lloc<0) ekb_tmp(-lloc+1,1)=ekb_tmp(-lloc+1+lmax,1) 434 if(lloc>0) ekb_so(lloc+1)=ekb_sr(lloc+1) 435 if(lloc>0) ekb_tmp(lmax+lloc+1,1)=ekb_tmp(lloc+1,1) 436 do ipsang=1,mpssoang 437 if(ekb_tmp(ipsang,1)>0) ekb_tmp(ipsang,1)= 1.d0 438 if(ekb_tmp(ipsang,1)<0) ekb_tmp(ipsang,1)=-1.d0 439 end do 440 441 ! v_ion is calculated in ffspl_tmp(:,:,1:mpsang,1) and v_so in 442 ! ffspl_tmp(:,:,mpsang+1:mpssoang,1) taking into account sqrt(ekb) 443 do i1=1,mqgrid 444 do i2=1,2 445 ffspl_tmp(i1,i2,1,1)=ffspl_sr(i1,i2,1)*sqrt(abs(ekb_sr(1))) 446 do ipsang=2,mpsang 447 ffspl_tmp(i1,i2,ipsang,1)=((ffspl_sr(i1,i2,ipsang)*& 448 & sqrt(abs(ekb_sr(ipsang)))*(ipsang-1))+& 449 & (ffspl_so(i1,i2,ipsang)*& 450 & sqrt(abs(ekb_so(ipsang)))*(ipsang)))& 451 & /(2.d0*ipsang-1) 452 ffspl_tmp(i1,i2,mpsang+ipsang-1,1)=(-ffspl_sr(i1,i2,ipsang)*& 453 & sqrt(abs(ekb_sr(ipsang)))+& 454 & ffspl_so(i1,i2,ipsang)*& 455 & sqrt(abs(ekb_so(ipsang))))*2.d0& 456 & /(2.d0*ipsang-1) 457 end do 458 end do 459 end do 460 ABI_DEALLOCATE(ekb_so) 461 ABI_DEALLOCATE(ffspl_so) 462 ABI_DEALLOCATE(vpspll_so) 463 ABI_DEALLOCATE(wfll_so) 464 465 ! The non local contribution is written as quadratic form of the vector 466 ! V=(v_ion,v_so) 467 ! t_V (Q1+Q2 L.S) V 468 ! with Q1= (1 0 ) et Q2=(0 1 ) 469 ! (0 l(l+1)/4) (1 -1/2) 470 ! The LS independent part is already diagonal. V is therefore built 471 ! putting v_so in the second projector of ffspl for the non spin-orbit 472 ! part and taking the eigenvalues of Q1 as new ekb (apart the sign) 473 do ipsang=2,mpsang 474 do i1=1,mqgrid 475 do i2=1,2 476 ffspl_tmp(i1,i2,ipsang,2)= ffspl_tmp(i1,i2,mpsang+ipsang-1,1) 477 end do 478 end do 479 ekb_tmp(ipsang,2)=ekb_tmp(mpsang+ipsang-1,1)*ipsang*(ipsang-1)*0.25d0 480 end do 481 482 ! For the spin orbit part, after diagonalisation of Q2, the eigenvectors 483 ! are: ((1-sqrt(17))/4 , 1) and ((1+sqrt(17))/4 ,1) 484 ! The passage matrix is therefore P=((1-sqrt(17))/4 (1+sqrt(17))/4) 485 ! ( 1 1 ) 486 ! t_P*Q2*P=( -sqrt(17)/2 0 ) 487 ! ( 0 sqrt(17)/2) 488 ! The diagonal values are the new ekb and the new ffspl are 489 ! P^-1 (v_ion) 490 ! (v_so ) 491 do ipsang=2,mpsang 492 do i1=1,mqgrid 493 do i2=1,2 494 ffspl_tmp(i1,i2,mpsang+ipsang-1,1)=-2.d0/sqrt(17.d0)*& 495 & (ffspl_tmp(i1,i2,ipsang,1)-& 496 & ((sqrt(17.d0)+1)*0.25d0)*& 497 ffspl_tmp(i1,i2,ipsang,2)) 498 ffspl_tmp(i1,i2,mpsang+ipsang-1,2)=2.d0/sqrt(17.d0)*& 499 & (ffspl_tmp(i1,i2,ipsang,1)+& 500 & ((sqrt(17.d0)-1)*0.25d0)*& 501 & ffspl_tmp(i1,i2,ipsang,2)) 502 end do 503 end do 504 ekb_tmp(mpsang+ipsang-1,1)=-(sqrt(17.d0)*0.5d0)*ekb_tmp(ipsang,1) 505 ekb_tmp(mpsang+ipsang-1,2)= (sqrt(17.d0)*0.5d0)*ekb_tmp(ipsang,1) 506 end do 507 508 end if 509 510 ABI_DEALLOCATE(ekb_sr) 511 ABI_DEALLOCATE(ffspl_sr) 512 513 ! FJ WARNING : No spin orbit if nproj>1 514 if (pspso==0) then 515 516 ! Read second wavefunction for second projection operator 517 ! (only read cases where nproj(ll)=2) 518 ! --also find highest l for which nproj(l)=2 519 lhigh=-1 520 do ipsang=1,min(lmax+1,mpsang) 521 if (nproj(ipsang)==2) then 522 lhigh=ipsang-1 523 read (tmp_unit,*, err=10, iomsg=errmsg) ll 524 if (ipsang/=ll+1) then 525 write(message, '(a,a,a,a,a,a,2i6,a,a)' )& 526 & 'Pseudopotential input file does not have',ch10,& 527 & 'angular momenta in order expected for second projection',& 528 & 'operator.',ch10,' Values are ',ipsang-1,ll,ch10,& 529 & 'Action : check your pseudopotential input file.' 530 MSG_ERROR(message) 531 end if 532 read (tmp_unit,*, err=10, iomsg=errmsg) wfll(:,ipsang) 533 ! DEBUG 534 ! write(std_out,*) 'WF second',ipsang-1,wfll(1,ipsang),wfll(2,ipsang),wfll(3,ipsang) 535 ! ENDDEBUG 536 else 537 wfll(:,ipsang)=0.0d0 538 end if 539 540 end do 541 542 ! Compute KB form factors and fit splines for second wf if any 543 if (lhigh>-1) then 544 call psp5nl(al,ekb_tmp(:,2),ffspl_tmp(:,:,:,2),lmax,& 545 & mmax,mpsang,mqgrid,qgrid,rad,vloc,vpspll,wfll) 546 end if 547 548 end if 549 550 ! Convert ekb and ffspl 551 iln=0 552 do ii=1,lmnmax 553 kk=indlmn(5,ii) 554 if (kk>iln) then 555 iln=kk 556 ll=indlmn(1,ii);nn=indlmn(3,ii) 557 if (indlmn(6,ii)==1) then 558 ekb(kk)=ekb_tmp(1+ll,nn) 559 ffspl(:,:,kk)=ffspl_tmp(:,:,1+ll,nn) 560 else 561 ekb(kk)=ekb_tmp(mpsang+ll,nn) 562 ffspl(:,:,kk)=ffspl_tmp(:,:,mpsang+ll,nn) 563 end if 564 end if 565 end do 566 567 ABI_DEALLOCATE(ekb_tmp) 568 ABI_DEALLOCATE(ffspl_tmp) 569 ABI_DEALLOCATE(wfll) 570 571 ! end of if concerning lloc 572 end if 573 574 ABI_DEALLOCATE(vpspll) 575 ABI_DEALLOCATE(rad) 576 ABI_DEALLOCATE(vloc) 577 578 !DEBUG 579 !write(std_out,*)' psp5in : vlspl(1,2)= ',vlspl(1,2) 580 !ENDDEBUG 581 582 return 583 584 ! Handle IO error 585 10 continue 586 MSG_ERROR(errmsg) 587 588 end subroutine psp5in