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