TABLE OF CONTENTS


ABINIT/m_psp6 [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psp6

FUNCTION

 Initialize pspcod=6 (Pseudopotentials from the fhi98pp code):

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XG, AF, GJ,FJ,MT, DRH)
  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_psp6
28 
29  use defs_basis
30  use m_splines
31  use m_errors
32  use m_abicore
33 
34  use m_numeric_tools,  only : smooth, ctrap
35  use m_psptk,          only : psp5lo, psp5nl, cc_derivatives
36 
37  implicit none
38 
39  private

m_psp6/psden [ Functions ]

[ Top ] [ m_psp6 ] [ Functions ]

NAME

 psden

FUNCTION

 Calculate a pseudo-density from an original density on a radial grid (regular or logarithmic)

INPUTS

  ilog=1 if grid is logarithmic, else 0
  mesh= dimension of nc
  nc(mesh)= density to be pseudized
  rc= cut-off radius
  rad(mesh) = radial mesh

OUTPUT

  ff(mesh)= pseudized density

  Optional:
    ff1(mesh)= 1st derivative of pseudo density (only r<rc modified)
    ff2(mesh)= 2nd derivative of pseudo density (only r<rc modified)

NOTES

    ff=exp(-(a+b.r^2+c.r^4))

PARENTS

      psp6cc

CHILDREN

      ctrap

SOURCE

576 subroutine psden(ilog,ff,mesh,nc,rc,rad,ff1,ff2)
577 
578 
579 !This section has been created automatically by the script Abilint (TD).
580 !Do not modify the following lines by hand.
581 #undef ABI_FUNC
582 #define ABI_FUNC 'psden'
583 !End of the abilint section
584 
585  implicit none
586 
587 !Arguments ------------------------------------
588 !scalars
589  integer,intent(in) :: ilog,mesh
590  real(dp),intent(in) :: rc
591 !arrays
592  real(dp),intent(in) :: nc(mesh),rad(mesh)
593  real(dp),intent(out) :: ff(mesh)
594  real(dp),intent(inout),optional :: ff1(mesh),ff2(mesh)
595 
596 !Local variables-------------------------------
597 !scalars
598  integer :: ii,nc1
599  real(dp) :: aa,aa1,aa2,bb,cc,c1,c3,f0,f0p,norm1,norm2,rc1,step
600 !arrays
601  real(dp),allocatable :: fpir(:),gg(:)
602 
603 ! *************************************************************************
604 
605  rc1=rc/four
606 
607  ABI_ALLOCATE(fpir,(mesh))
608  fpir(1:mesh)=four_pi*rad(1:mesh)**2
609  if (ilog==1) fpir(1:mesh)=fpir(1:mesh)*rad(1:mesh)
610 
611  if (ilog==0) then
612    step=rad(2)-rad(1)
613    nc1=int(rc1/step)+1
614    rc1=(nc1-1)*step
615  else if (ilog==1) then
616    step=log(rad(2)/rad(1))
617    nc1=int(log(rc1/rad(1))/step)+1
618    rc1=rad(nc1)
619  end if
620  ff(1:nc1)=nc(1:nc1)*fpir(1:nc1)
621  call ctrap(nc1,ff(1:nc1),step,c3)
622  if (ilog==1) c3=c3+half*ff(1)
623  f0=nc(nc1);c1=-log(f0)
624  f0p=half*(nc(nc1+1)-nc(nc1-1))/step
625 
626  ii=0;aa1=zero;norm1=c3+one
627  do while (norm1>c3.and.ii<100)
628    ii=ii+1;aa1=aa1+one
629    aa=c1-aa1*rc1**4+rc1*(f0p/f0+four*aa1*rc1**3)*half
630    bb=-half*(f0p/f0+four*aa1*rc1**3)/rc1
631    ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa1*rad(1:nc1)**4)
632    call ctrap(nc1,ff(1:nc1),step,norm1)
633    if (ilog==1) norm1=norm1+half*ff(1)
634  end do
635  if (ii==100) then
636    MSG_ERROR('Big pb 1 in psden !')
637  end if
638 
639  ii=0;aa2=zero;norm2=c3-one
640  do while (norm2<c3.and.ii<100)
641    ii=ii+1;aa2=aa2-one
642    aa=c1-aa2*rc1**4+rc1*(f0p/f0+four*aa2*rc1**3)*half
643    bb=-half*(f0p/f0+four*aa2*rc1**3)/rc1
644    ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-aa2*rad(1:nc1)**4)
645    call ctrap(nc1,ff(1:nc1),step,norm2)
646    if (ilog==1) norm2=norm2+half*ff(1)
647  end do
648  if (ii==100) then
649    MSG_ERROR('Big pb 2 in psden !')
650  end if
651 
652  do while (abs(norm2-c3)>tol10)
653 
654    cc=(aa1+aa2)*half
655    aa=c1-cc*rc1**4+rc1*(f0p/f0+four*cc*rc1**3)*half
656    bb=-half*(f0p/f0+four*cc*rc1**3)/rc1
657    ff(1:nc1)=fpir(1:nc1)*exp(-aa-bb*rad(1:nc1)**2-cc*rad(1:nc1)**4)
658    call ctrap (nc1,ff(1:nc1),step,norm2)
659    if (ilog==1) norm2=norm2+half*ff(1)
660    if ((norm1-c3)*(norm2-c3)>zero) then
661      aa1=cc
662      norm1=norm2
663    else
664      aa2=cc
665    end if
666 
667  end do ! while
668 
669  ff(1)=exp(-aa);if (ilog==1) ff(1)=ff(1)*exp(-bb*rad(1)**2-cc*rad(1)**4)
670  ff(2:nc1)=ff(2:nc1)/fpir(2:nc1)
671  if (nc1<mesh) ff(nc1+1:mesh)=nc(nc1+1:mesh)
672  if (present(ff1)) ff1(1:nc1)=-(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)*ff(1:nc1)
673  if (present(ff2)) ff2(1:nc1)=-(two*bb+12.0_dp*cc*rad(1:nc1)**2)*ff(1:nc1) &
674 & +(two*bb*rad(1:nc1)+four*cc*rad(1:nc1)**3)**2*ff(1:nc1)
675 
676  ABI_ALLOCATE(gg,(mesh))
677  gg(1:mesh)=fpir(1:mesh)*ff(1:mesh)
678  call ctrap(mesh,gg(1:mesh),step,norm1)
679  if (ilog==1) norm1=norm1+half*gg(1)
680  !write(std_out,*) 'psden: tild_nc integral= ',norm1
681  ABI_DEALLOCATE(gg)
682 
683  ABI_DEALLOCATE(fpir)
684 
685 end subroutine psden

m_psp6/psp6cc [ Functions ]

[ Top ] [ m_psp6 ] [ Functions ]

NAME

 psp6cc

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for the format 6 of pseudopotentials.

INPUTS

  mmax=maximum number of points in real space grid in the psp file
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  rchrg=cut-off radius for the core density
  znucl=nuclear number of atom as specified in psp file

OUTPUT

  xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives
  Optional output:
    vh_tnzc(mmax) = Hartree potential induced by density tild_[n_Z+n_core]
                    (pseudized [n_Z+n_core], where n_Z=ions, n_core=core electrons)
                    using a simple pseudization scheme

PARENTS

      psp6in

CHILDREN

      psden,smooth,spline,splint,vhtnzc

SOURCE

371 subroutine psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,&
372 &                 vh_tnzc) ! optional argument
373 
374 
375 !This section has been created automatically by the script Abilint (TD).
376 !Do not modify the following lines by hand.
377 #undef ABI_FUNC
378 #define ABI_FUNC 'psp6cc'
379 !End of the abilint section
380 
381  implicit none
382 
383 !Arguments ------------------------------------
384 !scalars
385  integer,intent(in) :: mmax,n1xccc
386  real(dp),intent(in) :: rchrg,znucl
387 !arrays
388  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
389  real(dp),intent(out),optional :: vh_tnzc(mmax)
390 
391 !Local variables-------------------------------
392 !scalars
393  integer :: i1xccc,irad
394  real(dp) :: der1,dern
395  character(len=500) :: errmsg
396 !arrays
397  real(dp),allocatable :: ff(:),ff1(:),ff2(:),ff3(:),gg(:),gg1(:),gg2(:),gg3(:)
398  real(dp),allocatable :: gg4(:),nc(:),rad(:),work(:),xx(:)
399 
400 !**********************************************************************
401 
402  ABI_ALLOCATE(ff,(mmax))
403  ABI_ALLOCATE(ff1,(mmax))
404  ABI_ALLOCATE(ff2,(mmax))
405  ABI_ALLOCATE(ff3,(mmax))
406  ABI_ALLOCATE(rad,(mmax))
407  ABI_ALLOCATE(gg,(n1xccc))
408  ABI_ALLOCATE(gg1,(n1xccc))
409  ABI_ALLOCATE(gg2,(n1xccc))
410  ABI_ALLOCATE(gg3,(n1xccc))
411  ABI_ALLOCATE(gg4,(n1xccc))
412  ABI_ALLOCATE(work,(n1xccc))
413  ABI_ALLOCATE(xx,(n1xccc))
414 
415 !read from pp file the model core charge (ff) and first (ff1) and
416 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid
417 !the input functions contain the 4pi factor, it must be rescaled.
418 
419  do irad=1,mmax
420    read(tmp_unit,*, err=10, iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad)
421    ff(irad)=ff(irad)/four_pi
422    ff1(irad)=ff1(irad)/four_pi
423    ff2(irad)=ff2(irad)/four_pi
424  end do
425 
426 !Optional output: VHartree(tild_[n_Z+n_core])
427  if (present(vh_tnzc)) then
428    ABI_ALLOCATE(nc,(mmax))
429    nc=ff ! n_core
430    call psden(1,ff,mmax,nc,rchrg,rad,ff1=ff1,ff2=ff2)
431    call vhtnzc(ff,rchrg,vh_tnzc,mmax,rad,znucl)
432    ABI_DEALLOCATE(nc)
433  end if
434 
435  rad(1)=zero
436 
437 !calculate third derivative ff3 on logarithmic grid
438  der1=ff2(1)
439  dern=ff2(mmax)
440  call spline(rad,ff1,mmax,der1,dern,ff3)
441 
442 !generate uniform mesh xx in the box cut by rchrg:
443 
444  do i1xccc=1,n1xccc
445    xx(i1xccc)=(i1xccc-1)* rchrg/dble(n1xccc-1)
446  end do
447 
448 !now interpolate core charge and derivatives on the uniform grid
449 !core charge, input=ff,  output=gg
450  call splint(mmax,rad,ff,ff2,n1xccc,xx,gg)
451 
452 !first derivative input=ff1, output=gg1
453  call splint(mmax,rad,ff1,ff3,n1xccc,xx,gg1)
454 
455 !normalize gg1
456  gg1(:)=gg1(:)*rchrg
457 
458 !now calculate second to fourth derivative by forward differences
459 !to avoid numerical noise uses a smoothing function
460 
461  call smooth(gg1,n1xccc,10)
462 
463  gg2(n1xccc)=zero
464  do i1xccc=1,n1xccc-1
465    gg2(i1xccc)=(gg1(i1xccc+1)-gg1(i1xccc))*dble(n1xccc-1)
466  end do
467 
468  call smooth(gg2,n1xccc,10)
469 
470  gg3(n1xccc)=zero
471  do i1xccc=1,n1xccc-1
472    gg3(i1xccc)=(gg2(i1xccc+1)-gg2(i1xccc))*dble(n1xccc-1)
473  end do
474 
475  call smooth(gg3,n1xccc,10)
476 
477  gg4(n1xccc)=zero
478  do i1xccc=1,n1xccc-1
479    gg4(i1xccc)=(gg3(i1xccc+1)-gg3(i1xccc))*dble(n1xccc-1)
480  end do
481 
482  call smooth(gg4,n1xccc,10)
483 
484 !write on xcc1d
485  xccc1d(:,1)=gg(:)
486  xccc1d(:,2)=gg1(:)
487  xccc1d(:,3)=gg2(:)
488  xccc1d(:,4)=gg3(:)
489  xccc1d(:,5)=gg4(:)
490 
491 !WARNING : fifth derivative not yet computed
492  xccc1d(:,6)=zero
493 
494 !note: the normalization condition is the following:
495 !4pi rchrg /dble(n1xccc-1) sum xx^2 xccc1d(:,1) = qchrg
496 !
497 !norm=zero
498 !do i1xccc=1,n1xccc
499 !norm = norm + four_pi*rchrg/dble(n1xccc-1)*&
500 !&             xx(i1xccc)**2*xccc1d(i1xccc,1)
501 !end do
502 !write(std_out,*) ' norm=',norm
503 !
504 !write(std_out,*)' psp1cc : output of core charge density and derivatives '
505 !write(std_out,*)'   xx          gg           gg1  '
506 !do i1xccc=1,n1xccc
507 !write(10, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,1),xccc1d(i1xccc,2)
508 !end do
509 !write(std_out,*)'   xx          gg2          gg3  '
510 !do i1xccc=1,n1xccc
511 !write(11, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,3),xccc1d(i1xccc,4)
512 !end do
513 !write(std_out,*)'   xx          gg4          gg5  '
514 !do i1xccc=1,n1xccc
515 !write(12, '(3es14.6)' ) xx(i1xccc),xccc1d(i1xccc,5),xccc1d(i1xccc,6)
516 !end do
517 !write(std_out,*)' psp1cc : debug done, stop '
518 !stop
519 !ENDDEBUG
520 
521  ABI_DEALLOCATE(ff)
522  ABI_DEALLOCATE(ff1)
523  ABI_DEALLOCATE(ff2)
524  ABI_DEALLOCATE(ff3)
525  ABI_DEALLOCATE(gg)
526  ABI_DEALLOCATE(gg1)
527  ABI_DEALLOCATE(gg2)
528  ABI_DEALLOCATE(gg3)
529  ABI_DEALLOCATE(gg4)
530  ABI_DEALLOCATE(rad)
531  ABI_DEALLOCATE(work)
532  ABI_DEALLOCATE(xx)
533 
534  return
535 
536  ! Handle IO error
537  10 continue
538  MSG_ERROR(errmsg)
539 
540 end subroutine psp6cc

m_psp6/psp6cc_drh [ Functions ]

[ Top ] [ m_psp6 ] [ Functions ]

NAME

 psp6cc_drh

FUNCTION

 Compute the core charge density, for use in the XC core
 correction, following the function definition valid
 for the format 6 of pseudopotentials.
 Version modified by DHamann, with consistent treatment
 of the derivatives in this routine and the remaining of the code.

INPUTS

  mmax=maximum number of points in real space grid in the psp file
  n1xccc=dimension of xccc1d ; 0 if no XC core correction is used
  rchrg=cut-off radius for the core density

OUTPUT

  xccc1d(n1xccc,6)= 1D core charge function and its five first derivatives

PARENTS

      psp6in

CHILDREN

      cc_derivatives

NOTES

 Test version by DRH - requires very smooth model core charge

SOURCE

835 subroutine psp6cc_drh(mmax,n1xccc,rchrg,xccc1d)
836 
837 
838 !This section has been created automatically by the script Abilint (TD).
839 !Do not modify the following lines by hand.
840 #undef ABI_FUNC
841 #define ABI_FUNC 'psp6cc_drh'
842 !End of the abilint section
843 
844  implicit none
845 
846 !Arguments ------------------------------------
847 !scalars
848  integer,intent(in) :: mmax,n1xccc
849  real(dp),intent(in) :: rchrg
850 !arrays
851  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
852 
853 !Local variables-------------------------------
854 !scalars
855  integer :: irad
856  character(len=500) :: errmsg
857 !arrays
858  real(dp),allocatable :: ff(:),ff1(:),ff2(:),rad(:)
859 
860 !**********************************************************************
861 
862  ABI_ALLOCATE(ff,(mmax))
863  ABI_ALLOCATE(ff1,(mmax))
864  ABI_ALLOCATE(ff2,(mmax))
865  ABI_ALLOCATE(rad,(mmax))
866 
867 !
868 !read from pp file the model core charge (ff) and first (ff1) and
869 !second (ff2) derivative on logarithmic mesh mmax; rad is the radial grid
870 !the input functions contain the 4pi factor, it must be rescaled.
871 
872  !write(std_out,'(a,2i6)') 'drh:psp6cc_drh - mmax,n1xccc',mmax,n1xccc
873  do irad=1,mmax
874    read(tmp_unit,*,err=10,iomsg=errmsg) rad(irad),ff(irad),ff1(irad),ff2(irad)
875    ff(irad)=ff(irad)/4.d0/pi
876    ff1(irad)=ff1(irad)/4.d0/pi
877    ff2(irad)=ff2(irad)/4.d0/pi
878  end do
879  rad(1)=0.d0
880 
881  call cc_derivatives(rad,ff,ff1,ff2,mmax,n1xccc,rchrg,xccc1d)
882 
883  ABI_DEALLOCATE(ff)
884  ABI_DEALLOCATE(ff1)
885  ABI_DEALLOCATE(ff2)
886  ABI_DEALLOCATE(rad)
887 
888  return
889 
890  ! Handle IO error
891  10 continue
892  MSG_ERROR(errmsg)
893 
894 end subroutine psp6cc_drh

m_psp6/psp6in [ Functions ]

[ Top ] [ m_psp6 ] [ Functions ]

NAME

 psp6in

FUNCTION

 Initialize pspcod=6 (Pseudopotentials from the fhi98pp code):
 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
   angular momentum of nonlocal pseudopotential
  mpsang= 1+maximum 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
  optnlxccc=option for nl XC core correction (input variable)
  positron=0 if electron GS calculation
          1 if positron GS calculation
          2 if electron GS calculation in presence of the positron
  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)
  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(mpsang)=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

PARENTS

      pspatm

CHILDREN

      psp5lo,psp5nl,psp6cc,psp6cc_drh,spline,wrtout

SOURCE

105 subroutine psp6in(ekb,epsatm,ffspl,indlmn,lloc,lmax,lmnmax,lnmax,&
106 &                  mmax,mpsang,mqgrid,nproj,n1xccc,optnlxccc,positron,qchrg,qgrid,&
107 &                  useylm,vlspl,xcccrc,xccc1d,zion,znucl)
108 
109 
110 !This section has been created automatically by the script Abilint (TD).
111 !Do not modify the following lines by hand.
112 #undef ABI_FUNC
113 #define ABI_FUNC 'psp6in'
114 !End of the abilint section
115 
116  implicit none
117 
118 !Arguments ------------------------------------
119 !scalars
120  integer,intent(in) :: lloc,lmax,lmnmax,lnmax,mmax,mpsang,mqgrid,n1xccc
121  integer,intent(in) :: optnlxccc,positron,useylm
122  real(dp),intent(in) :: zion,znucl
123  real(dp),intent(out) :: epsatm,qchrg,xcccrc
124 !arrays
125  integer,intent(out) :: indlmn(6,lmnmax),nproj(mpsang)
126  real(dp),intent(in) :: qgrid(mqgrid)
127  real(dp),intent(out) :: ekb(lnmax),vlspl(mqgrid,2) !vz_i
128  real(dp),intent(inout) :: ffspl(mqgrid,2,lnmax) !vz_i
129  real(dp),intent(inout) :: xccc1d(n1xccc,6) !vz_i
130 
131 !Local variables-------------------------------
132 !scalars
133  integer :: ii,index,ipsang,irad,jj,jpsang,mm,mmax2
134  real(dp) :: al,al_announced,amesh,fchrg,ratio,rchrg,yp1,ypn
135  character(len=3) :: testxc
136  character(len=500) :: message,errmsg
137 !arrays
138  real(dp),allocatable :: ekb_tmp(:),ffspl_tmp(:,:,:),rad(:),vloc(:)
139 !real(dp),allocatable :: radbis
140  real(dp),allocatable :: vpspll(:,:),wfll(:,:),work_space(:),work_spl(:)
141 
142 ! ***************************************************************************
143 
144 !File format of formatted fhi psp input, as adapted for use
145 !by the ABINIT code (the 3 first lines
146 !have already been read in calling -pspatm- routine) :
147 
148 !(1) title (character) line
149 !(2) znucl,zion,pspdat
150 !(3) pspcod,pspxc,lmax,lloc,mmax,r2well
151 !(4) rchrg,fchrg,qchrg
152 !Note : prior to version 2.2, this 4th line started with  4--  ,
153 !and no core-correction was available.
154 !(5)-(18) -empty-
155 !(19) mmax, amesh ( mesh increment r(m+1)/r(m) )
156 !Then, for ll=0,lmax :
157 !for  irad=1,mmax  : irad, r(irad), upsp(irad,ll), vpsp(irad,ll)
158 
159  read (tmp_unit, '(a3)', err=10, iomsg=errmsg) testxc
160  if(testxc/='4--')then
161    backspace(tmp_unit)
162    read (tmp_unit,*, err=10, iomsg=errmsg) rchrg,fchrg,qchrg
163    write(message, '(3f20.14,t64,a)' ) rchrg,fchrg,qchrg,'rchrg,fchrg,qchrg'
164    call wrtout(ab_out,message,'COLL')
165    call wrtout(std_out,  message,'COLL')
166  else
167    write(message, '(a)' ) '  No XC core correction.'
168    call wrtout(ab_out,message,'COLL')
169    call wrtout(std_out,  message,'COLL')
170    rchrg=zero ; fchrg=zero ; qchrg=zero
171  end if
172  do ii=5,18
173    read(tmp_unit,*, err=10, iomsg=errmsg)
174  end do
175 
176  if (positron==1.and.abs(fchrg)<=tol14) then
177    write(message,'(5a)')&
178 &   'You can only perform positronic ground-state calculations (positron=1)',ch10,&
179 &   'using fhi pseudopotentials with a core density (fchrg>0)',ch10,&
180 &   'Action: change your psp file (add fchrg>0).'
181    MSG_ERROR(message)
182  end if
183 !--------------------------------------------------------------------
184 !Will now proceed at the reading of pots and wfs
185 
186 !rad(:)=radial grid r(i)
187 !vpspll(:,1),...,vpspll(:,4)=nonlocal pseudopotentials
188 !wfll(:,1),...,wfll(:,4)=reference config. wavefunctions
189 
190  ABI_ALLOCATE(rad,(mmax))
191  ABI_ALLOCATE(vpspll,(mmax,mpsang))
192  ABI_ALLOCATE(wfll,(mmax,mpsang))
193 
194 !Read atomic pseudopotential for each l, filling up arrays vpspll
195 !and wfll. Also set up rad array (actually read more than once)
196 !Note: put each l into vpspll(:,l+1)
197  do ipsang=1,lmax+1
198    nproj(ipsang)=1
199    read(tmp_unit,*, err=10, iomsg=errmsg)mmax2,amesh
200    if(ipsang==1)then
201      write(message, '(f10.6,t20,a)' ) amesh,' amesh (Hamman grid)'
202      al_announced=log(amesh)
203      call wrtout(ab_out,message,'COLL')
204      call wrtout(std_out,  message,'COLL')
205    end if
206    do irad=1,mmax
207      read(tmp_unit,*, err=10, iomsg=errmsg)jj,rad(irad),wfll(irad,ipsang),vpspll(irad,ipsang)
208    end do
209  end do
210 
211 
212 !Generate core charge function and derivatives, if needed
213  if(fchrg>tol14)then
214 
215    if (positron==1) then
216      call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl,vh_tnzc=vpspll(:,lloc+1))
217    else if(optnlxccc==1)then
218      call psp6cc(mmax,n1xccc,rchrg,xccc1d,znucl)
219    else if(optnlxccc==2)then
220      call psp6cc_drh(mmax,n1xccc,rchrg,xccc1d)
221    end if
222 
223 ! The core charge function for pspcod=6 becomes zero beyond rchrg.
224 ! Thus xcccrc must be set equal to rchrg.
225    xcccrc=rchrg
226  else
227    xccc1d(:,:)=zero
228    xcccrc=zero
229  end if
230 
231 !Compute in real(dp) al : the announced amesh is inaccurate.
232  ratio=rad(mmax)/rad(1)
233  al=log(ratio)/dble(mmax-1)
234 
235 !vloc(:)=Vlocal(r), lloc=0, 1, or 2 or -1 for avg.
236  ABI_ALLOCATE(vloc,(mmax))
237 !Copy appropriate nonlocal psp for use as local one
238  vloc( 1:mmax ) = vpspll( 1:mmax , lloc+1 )
239 
240 !--------------------------------------------------------------------
241 !Carry out calculations for local (lloc) pseudopotential.
242 !Obtain Fourier transform (1-d sine transform) to get q^2 V(q).
243 
244  call psp5lo(al,epsatm,mmax,mqgrid,qgrid,&
245 & vlspl(:,1),rad,vloc,yp1,ypn,zion)
246 
247 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
248  ABI_ALLOCATE(work_space,(mqgrid))
249  ABI_ALLOCATE(work_spl,(mqgrid))
250  call spline (qgrid,vlspl(:,1),mqgrid,yp1,ypn,work_spl)
251  vlspl(:,2)=work_spl(:)
252  ABI_DEALLOCATE(work_space)
253  ABI_DEALLOCATE(work_spl)
254 
255 !--------------------------------------------------------------------
256 !Take care of non-local part
257 
258  ABI_ALLOCATE(ekb_tmp,(mpsang))
259  ABI_ALLOCATE(ffspl_tmp,(mqgrid,2,mpsang))
260 
261 !Zero out all Kleinman-Bylander energies to initialize
262  ekb_tmp(:)=zero
263  ekb(:)=zero
264 
265 !Allow for option of no nonlocal corrections (lloc=lmax=0)
266  if (lloc==0.and.lmax==0) then
267    write(message, '(a,f5.1)' ) ' Note: local psp for atom with Z=',znucl
268    call wrtout(ab_out,message,'COLL')
269    call wrtout(std_out,  message,'COLL')
270  else
271 
272 !  ----------------------------------------------------------------------
273 !  Compute KB form factors and fit splines
274 
275    call psp5nl(al,ekb_tmp,ffspl_tmp,lmax,mmax,mpsang,mqgrid,qgrid,rad,vloc, vpspll,wfll)
276 
277  end if
278 
279  jj=0;index=0;indlmn(:,:)=0
280  do ipsang=1,lmax+1
281 !  nproj had been set at 1, by default
282    if(abs(ekb_tmp(ipsang))<tol10)then
283      nproj(ipsang)=0
284    end if
285 !  Possible values for nproj in this routine : 0 or 1.
286    if(nproj(ipsang)==1)then
287      if (useylm==1) then
288        jj=jj+1
289        do mm=1,2*ipsang-1
290          index=index+1
291          indlmn(1,index)=ipsang-1
292          indlmn(2,index)=mm-ipsang
293          indlmn(3,index)=1
294          indlmn(4,index)=mm+(ipsang-1)*(ipsang-1)
295          indlmn(5,index)=jj
296          indlmn(6,index)=1
297        end do
298      else
299        jj=jj+1
300        index=index+1
301        indlmn(1,index)=ipsang-1
302        indlmn(2,index)=0
303        indlmn(3,index)=1
304        indlmn(4,index)=ipsang+(ipsang-1)*(ipsang-1)
305        indlmn(5,index)=jj
306        indlmn(6,index)=1
307      end if
308    end if
309  end do
310 !Transfer ekb and ffspl to their definitive location
311  jpsang=1
312  do ipsang=1,lmax+1
313    if(nproj(ipsang)/=0)then
314      ekb(jpsang)=ekb_tmp(ipsang)
315      ffspl(:,:,jpsang)=ffspl_tmp(:,:,ipsang)
316      jpsang=jpsang+1
317      if(jpsang>lnmax)then
318        write(message,'(3a,2i6)')&
319 &       'Problem with the dimension of the ekb and ffspl arrays.',ch10,&
320 &       'ipsang,lnmax=',ipsang,lnmax
321      end if
322    end if
323  end do
324 
325  ABI_DEALLOCATE(ekb_tmp)
326  ABI_DEALLOCATE(ffspl_tmp)
327  ABI_DEALLOCATE(vpspll)
328  ABI_DEALLOCATE(rad)
329  ABI_DEALLOCATE(vloc)
330  ABI_DEALLOCATE(wfll)
331 
332  return
333 
334  ! Handle IO error
335  10 continue
336  MSG_ERROR(errmsg)
337 
338 end subroutine psp6in

m_psp6/vhtnzc [ Functions ]

[ Top ] [ m_psp6 ] [ Functions ]

NAME

 vhtnzc

FUNCTION

 Compute VHartree(tild[n_Z+n_core]) from input ncore

INPUTS

  mesh=dimension of radial mesh
  nc= core density (to be pseudized)
  rad(mesh)=radial mesh
  rc=cut-off radius
  znucl=nuclear number of atom as specified in psp file

OUTPUT

  vhtnzc(mesh) = hartree potential induced by density tild[n_Z+n_core] (pseudo core density + nucleus)

PARENTS

      psp6cc

CHILDREN

      ctrap

SOURCE

713 subroutine vhtnzc(nc,rc,vh_tnzc,mesh,rad,znucl)
714 
715 
716 !This section has been created automatically by the script Abilint (TD).
717 !Do not modify the following lines by hand.
718 #undef ABI_FUNC
719 #define ABI_FUNC 'vhtnzc'
720 !End of the abilint section
721 
722  implicit none
723 
724 !Arguments ------------------------------------
725 !scalars
726  integer,intent(in) :: mesh
727  real(dp),intent(in) :: znucl
728  real(dp),intent(in) :: rc
729 !arrays
730  real(dp),intent(in) :: nc(mesh),rad(mesh)
731  real(dp),intent(out) :: vh_tnzc(mesh)
732 
733 !Local variables-------------------------------
734 !scalars
735  integer :: ir,nc1
736  real(dp) :: gnorm,rc1,step,yp1,yp2,yp3
737 !arrays
738  real(dp),allocatable :: den1(:),den2(:),den3(:),den4(:),nzc(:),rvhn(:),shapefunc(:)
739 
740 ! *************************************************************************
741 
742  rc1=rc/four
743 
744  step=log(rad(2)/rad(1))
745  nc1=int(log(rc1/rad(1))/step)+1
746  rc1=rad(nc1)
747 
748  ABI_ALLOCATE(shapefunc,(mesh))
749  shapefunc(1)=one
750  shapefunc(2:nc1)=(sin(pi*rad(2:nc1)/rc1)/(pi*rad(2:nc1)/rc1))**2
751  if (nc1<mesh) shapefunc(nc1+1:mesh)=zero
752 
753  ABI_ALLOCATE(den1,(mesh))
754  den1(1:mesh)=four_pi*shapefunc(1:mesh)*rad(1:mesh)**3
755  call ctrap(mesh,den1,step,gnorm)
756  gnorm =one/gnorm
757  ABI_DEALLOCATE(den1)
758 
759  ABI_ALLOCATE(nzc,(mesh))
760  nzc(1:mesh)=four*pi*nc(1:mesh)*rad(1:mesh)**2-four_pi*shapefunc(1:mesh)*rad(1:mesh)**2*znucl*gnorm
761  ABI_DEALLOCATE(shapefunc)
762 
763  ABI_ALLOCATE(rvhn,(mesh))
764  rvhn(1)=zero
765 
766  ABI_ALLOCATE(den1,(mesh))
767  ABI_ALLOCATE(den2,(mesh))
768  ABI_ALLOCATE(den3,(mesh))
769  ABI_ALLOCATE(den4,(mesh))
770 
771  den1(1)=zero;den2(1)=zero
772  do ir=2,mesh
773    den1(ir)= rad(ir)*nzc(ir)
774    den2(ir)= den1(ir)/rad(ir)
775  end do
776 
777 !For first few points do stupid integral
778  den3(1)=zero;den4(1)=zero
779  do ir=2,mesh
780    call ctrap(ir,den1(1:ir),step,den3(ir))
781    call ctrap(ir,den2(1:ir),step,den4(ir))
782  end do
783 
784  do ir=1,mesh
785    rvhn(ir)=den3(ir)+rad(ir)*(den4(mesh)-den4(ir))
786  end do
787 
788  ABI_DEALLOCATE(den1)
789  ABI_DEALLOCATE(den2)
790  ABI_DEALLOCATE(den3)
791  ABI_DEALLOCATE(den4)
792 
793  vh_tnzc(2:mesh)=rvhn(2:mesh)/rad(2:mesh)
794  yp2=(vh_tnzc(3)-vh_tnzc(2))/(rad(3)-rad(2))
795  yp3=(vh_tnzc(4)-vh_tnzc(3))/(rad(4)-rad(3))
796  yp1=yp2+(yp2-yp3)*rad(2)/(rad(3)-rad(2))
797  vh_tnzc(1)=vh_tnzc(2)-(yp1+yp2)*rad(2)
798 
799  ABI_DEALLOCATE(nzc)
800  ABI_DEALLOCATE(rvhn)
801 
802 end subroutine vhtnzc