TABLE OF CONTENTS


ABINIT/m_psp5 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psp5

FUNCTION

 Initialize pspcod=5 ("Phoney pseudopotentials" with Hamman grid):

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_psp5
28 
29  use defs_basis
30  use m_splines
31  use m_errors
32  use m_abicore
33 
34  use m_psptk,           only : psp1cc, psp5lo, psp5nl
35 
36  implicit none
37 
38  private

ABINIT/psp5in [ Functions ]

[ Top ] [ 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

PARENTS

      pspatm

CHILDREN

      psp1cc,psp5lo,psp5nl,spline,wrtout

SOURCE

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