TABLE OF CONTENTS


ABINIT/m_psp_hgh [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psp_hgh

FUNCTION

 Initialize pspcod=2, 3, 10 pseudopotentials (GTH)

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT, FD, PT)
  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_psp_hgh
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_splines
33  use m_abicore
34  use m_errors
35 
36  use m_special_funcs,  only : abi_derfc
37 #if defined HAVE_BIGDFT
38  use BigDFT_API, only: atomic_info
39 #endif
40  use m_wvl_descr_psp,  only : wvl_descr_psp_fill
41 
42  implicit none
43 
44  private

ABINIT/psp2lo [ Functions ]

[ Top ] [ Functions ]

NAME

 psp2lo

FUNCTION

 Treat local part of Goedecker-Teter-Hutter pseudopotentials (pspcod=2),
 as well as Hartwigsen-Goedecker-Hutter pseudopotentials (pspcod=3)

INPUTS

  cc1,2,3,4=parameters from analytic pseudopotential form
  mqgrid=number of grid points in q from 0 to qmax.
  qgrid(mqgrid)=values of q (or G) on grid from 0 to qmax (bohr^-1)
                if vlspl_recipSpace is .true. else values of r on grid from
                0 to 2pi / qmax * mqgrid_ff (bohr).
  rloc=local pseudopotential core radius (bohr)
  vlspl_recipSpace= .true. if computation of vlspl is done in reciprocal space
  zion=valence charge of atom
  parameters for local potential: rloc,c1,c2,c3,c4

OUTPUT

  dvloc(mqgrid)=dVloc(r)/dr (only allocated if vlspl_recipSpace is false).
  epsatm=$4\pi\int[r^2 (v(r)+\frac{Zv}{r} dr]$
  q2vq(mqgrid)&=&q^2 v(q) \nonumber \\
  &=&-Zv/\pi
   +q^2 4\pi\int[(\frac{\sin(2\pi qr)}{2\pi qr})(r^2 v(r)+r Zv)dr]\nonumber\\
  &=&\exp(-K^2*rloc^2/2) \nonumber \\
  &&   *(-\frac{zion}{\pi}+(\frac{K^2*rloc^3}{\sqrt{2*\pi}}*
       (c1+c2*(3-(rloc*K)^2) \nonumber \\
  &&    +c3*(15-10(rloc*K)^2+(rloc*K)^4) \nonumber \\
  &&    +c4*(105-105*(rloc*K)^2+21*(rloc*K)^4-(rloc*K)^6)) \nonumber
 for GTH vloc with $K=(2\pi q)$.
  yp1,ypn=derivative of q^2 v(q) wrt q at q=0 and q=qmax
   (needed for spline fitter).

PARENTS

      psp10in,psp2in,psp3in

CHILDREN

      wrtout

SOURCE

461 subroutine psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,mqgrid,qgrid,q2vq,&
462 &  rloc,vlspl_recipSpace,yp1,ypn,zion)
463 
464 
465 !This section has been created automatically by the script Abilint (TD).
466 !Do not modify the following lines by hand.
467 #undef ABI_FUNC
468 #define ABI_FUNC 'psp2lo'
469 !End of the abilint section
470 
471  implicit none
472 
473 !Arguments ------------------------------------
474 !scalars
475  integer,intent(in) :: mqgrid
476  real(dp),intent(in) :: cc1,cc2,cc3,cc4,rloc,zion
477  real(dp),intent(out) :: epsatm,yp1,ypn
478  logical,intent(in) :: vlspl_recipSpace
479 !arrays
480  real(dp),intent(in) :: qgrid(mqgrid)
481  real(dp),intent(out) :: dvloc(mqgrid),q2vq(mqgrid)
482 
483 !Local variables-------------------------------
484 !scalars
485  integer :: iqgrid
486  real(dp) :: erfValue,gaussValue,polyValue,qmax,rq,rq2
487  character(len=500) :: message
488 
489 ! *************************************************************************
490 
491 !Compute epsatm = lim(q->0) [Vloc(q) + zion/(Pi*q^2)]
492  epsatm=2.d0*pi*rloc**2*zion+(2.d0*pi)**(1.5d0)*rloc**3*&
493 & (cc1+3.d0*cc2+15.d0*cc3+105.d0*cc4)
494 
495 !If vlspl_recipSpace is .true., we compute V(q)*q^2 in reciprocal space,
496 !else we compute V(r) in real space.
497  if (vlspl_recipSpace) then
498    write(message, '(a)' ) '-  Local part computed in reciprocal space.'
499    call wrtout(ab_out,message,'COLL')
500    call wrtout(std_out,  message,'COLL')
501 
502 !  d(q^2*V(q))/d(q) at q=0 and q=qmax
503    qmax=qgrid(mqgrid)
504    rq2=(2.d0*pi*qmax*rloc)**2
505    yp1=0.d0
506    ypn= (2.d0*pi*qmax*rloc**2)*exp(-0.5d0*rq2)* &
507 &   (2.d0*zion + sqrt(2.d0*pi)*rloc*&
508 &   (cc1*(2.d0-rq2) + cc2*(6.d0-7.d0*rq2+rq2**2) +&
509 &   cc3*(30.d0-55.d0*rq2+16.d0*rq2**2-rq2**3) +&
510 &   cc4*(210.d0-525.d0*rq2+231.d0*rq2**2-29.d0*rq2**3+rq2**4)))
511 !  ypn has been tested against Maple-derived expression.
512 
513 !  Compute q^2*vloc(q) on uniform grid
514    do iqgrid=1,mqgrid
515      rq2=(2.d0*pi*qgrid(iqgrid)*rloc)**2
516      q2vq(iqgrid)=exp(-0.5d0*rq2)*(-zion/pi+rq2*(rloc/sqrt(2.d0*pi)) *&
517 &     ( cc1 + cc2*(3.d0-rq2) + cc3*(15.d0-10.d0*rq2+rq2**2) +&
518 &     cc4*(105.d0-rq2*(105.d0-rq2*(21.d0-rq2)))  ))
519    end do
520  else
521    write(message, '(a)' ) '-  Local part computed in real space.'
522    call wrtout(ab_out,message,'COLL')
523    call wrtout(std_out,  message,'COLL')
524 
525 !  Compute derivatives for splines computations
526    yp1 = 0.d0
527    rq2 = (qgrid(mqgrid) / rloc) ** 2
528    erfValue = abi_derfc(sqrt(0.5d0 * rq2))
529    ypn = - 2.0d0 * zion / sqrt(2.d0 * pi) / qgrid(mqgrid) / rloc
530    ypn = ypn - rq2 * (cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3) / qgrid(mqgrid)
531    ypn = ypn + (2.d0 * cc2 * rq2 + 4.d0 * cc3 * rq2 ** 2 + 6.d0 * cc4 * rq2 ** 3) / qgrid(mqgrid)
532    ypn = ypn * exp(-0.5d0 * rq2)
533    ypn = ypn + zion / qgrid(mqgrid) ** 2 * erfValue
534 !  Note that ypn has been calculated on a full-proof a4 paper sheet.
535 
536 !  Compute local potential and its first derivatives.
537    do iqgrid = 1, mqgrid, 1
538      rq2 = (qgrid(iqgrid) / rloc) ** 2
539 !    Compute erf() part
540 !    Case r = 0
541      gaussValue = exp(-0.5d0 * rq2)
542      if (qgrid(iqgrid) == 0.d0) then
543        q2vq(iqgrid) = -zion / rloc * sqrt(2.d0 / pi)
544        dvloc(iqgrid) = 0.d0
545      else
546        erfValue = abi_derfc(sqrt(0.5d0 * rq2))
547        q2vq(iqgrid) = -zion / qgrid(iqgrid) * (1.0d0 - erfValue)
548        dvloc(iqgrid) = - sqrt(2.d0 / pi) * zion * gaussValue / (qgrid(iqgrid) * rloc) - &
549 &       q2vq(iqgrid) / qgrid(iqgrid)
550      end if
551 !    Add the gaussian part
552      polyValue = cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3
553      q2vq(iqgrid) = q2vq(iqgrid) + gaussValue * polyValue
554      rq = qgrid(iqgrid) / rloc
555      dvloc(iqgrid) = dvloc(iqgrid) - qgrid(iqgrid) / rloc ** 2 * gaussValue * polyValue + &
556 &     gaussValue * (2.0d0 * cc2 * rq / rloc + 3.0d0 * cc3 * rq ** 3 / rloc + &
557 &     6.0d0 * cc4 * rq ** 5 / rloc)
558    end do
559 
560    write(message, '(a,f12.7,a,a,f12.7,a,a,a,f12.7)' ) &
561 &   '  | dr spline step is : ', qgrid(2), ch10, &
562 &   '  | r > ', qgrid(mqgrid) ,' is set to 0.', ch10, &
563 &   '  | last non-nul potential value is : ', q2vq(mqgrid)
564    call wrtout(ab_out,message,'COLL')
565    call wrtout(std_out,  message,'COLL')
566  end if
567 
568 end subroutine psp2lo

ABINIT/psp3in [ Functions ]

[ Top ] [ Functions ]

NAME

 psp3in

FUNCTION

 Initialize pspcod=3 pseudopotentials (HGH psps PRB58,3641(1998) [[cite:Hartwigsen1998]]):
 continue to read the file, then compute the corresponding
 local and non-local potentials.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  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 characteristics of the psp
   if =3 : this input requires HFN characteristics of the psp
  ipsp=id in the array of the pseudo-potential.
  zion=nominal valence 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
  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$ (hartree)
  ffspl(mqgrid_ff,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
  vlspl(mqgrid_ff,2)=q^2 Vloc(q) and second derivatives from spline fit

SIDE EFFECTS

  Input/output
  lmax : at input =value of lmax mentioned at the second line of the psp file
    at output= 1
  psps <type(pseudopotential_type)>=at output, values depending on the read
                                    pseudo are set.
   | lmnmax(IN)=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(IN)=max. number of (l,n) components over all type of psps
   |           angular momentum of nonlocal pseudopotential
   | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
   | mpssoang(IN)= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
   | mqgrid_ff(IN)=dimension of q (or G) grid for arrays.
   | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
   | useylm(IN)=governs the way the nonlocal operator is to be applied:
   |            1=using Ylm, 0=using Legendre polynomials

PARENTS

      pspatm

CHILDREN

      psp2lo,psp3nl,spline,wrtout,wvl_descr_psp_fill

SOURCE

 631 subroutine psp3in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion)
 632 
 633 
 634 !This section has been created automatically by the script Abilint (TD).
 635 !Do not modify the following lines by hand.
 636 #undef ABI_FUNC
 637 #define ABI_FUNC 'psp3in'
 638 !End of the abilint section
 639 
 640  implicit none
 641 
 642 !Arguments ------------------------------------
 643 !scalars
 644  integer,intent(in) :: ipsp,pspso
 645  integer,intent(inout) :: lmax
 646  real(dp),intent(in) :: zion
 647  real(dp),intent(out) :: epsatm
 648  type(dataset_type),intent(in) :: dtset
 649  type(pseudopotential_type),intent(inout) :: psps
 650 !arrays
 651  integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang)
 652  real(dp),intent(inout) :: ekb(psps%lnmax),ffspl(psps%mqgrid_ff,2,psps%lnmax)!vz_i
 653  real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2)
 654 
 655 !Local variables-------------------------------
 656 !scalars
 657  integer :: iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nso
 658  real(dp) :: cc1,cc2,cc3,cc4,h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s
 659  real(dp) :: k11d,k11f,k11p,k22d,k22p,k33d,k33p,rloc
 660  real(dp) :: rrd,rrf,rrp,rrs,yp1,ypn
 661  character(len=500) :: message,errmsg
 662 !arrays
 663  real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:)
 664  real(dp),allocatable :: ffspl_sr(:,:,:,:),work_space(:),work_spl(:)
 665 
 666 ! ***************************************************************************
 667 
 668 !Set various terms to 0 in case not defined below
 669 !HGH values
 670  rloc=zero ; rrs=zero  ; h11p=zero ; k33p=zero ; k11d=zero;
 671  cc1=zero  ; h11s=zero ; h22p=zero ; rrd=zero  ; k22d=zero;
 672  cc2=zero  ; h22s=zero ; h33p=zero ; h11d=zero ; k33d=zero;
 673  cc3=zero  ; h33s=zero ; k11p=zero ; h22d=zero ; h11f=zero;
 674  cc4=zero  ; rrp=zero  ; k22p=zero ; h33d=zero ; k11f=zero;
 675  rrf=zero
 676  nproj(1:psps%mpssoang)=0
 677 
 678 !Read and write different lines of the pseudopotential file
 679 
 680  read (tmp_unit,*,err=10,iomsg=errmsg) rloc,cc1,cc2,cc3,cc4
 681  write(message, '(a,f12.7)' ) ' rloc=',rloc
 682  call wrtout(ab_out,message,'COLL')
 683  call wrtout(std_out,  message,'COLL')
 684  write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' cc1 =',cc1,'; cc2 =',cc2,'; cc3 =',cc3,'; cc4 =',cc4
 685  call wrtout(ab_out,message,'COLL')
 686  call wrtout(std_out,  message,'COLL')
 687 
 688 !For the time being, the s state line must be present and is read,
 689 !even for local pseudopotentials (zero must appear)
 690  read (tmp_unit,*,err=10,iomsg=errmsg) rrs,h11s,h22s,h33s
 691  write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrs =',rrs,'; h11s=',h11s,'; h22s=',h22s,'; h33s=',h33s
 692  call wrtout(ab_out,message,'COLL')
 693  call wrtout(std_out,  message,'COLL')
 694 
 695  if (lmax > 0) then
 696 
 697    read (tmp_unit,*,err=10,iomsg=errmsg) rrp,h11p,h22p,h33p
 698    write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrp =',rrp,'; h11p=',h11p,'; h22p=',h22p,'; h33p=',h33p
 699    call wrtout(ab_out,message,'COLL')
 700    call wrtout(std_out,  message,'COLL')
 701 
 702    read (tmp_unit,*,err=10,iomsg=errmsg) k11p,k22p,k33p
 703    write(message, '(a,f12.7,a,f12.7,a,f12.7)' )'                    k11p=',k11p,'; k22p=',k22p,'; k33p=',k33p
 704    call wrtout(ab_out,message,'COLL')
 705    call wrtout(std_out,  message,'COLL')
 706 
 707  end if
 708 
 709  if (lmax > 1) then
 710    read (tmp_unit,*,err=10,iomsg=errmsg) rrd,h11d,h22d,h33d
 711    write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrd =',rrd,'; h11d=',h11d,'; h22d=',h22d,'; h33d=',h33d
 712    call wrtout(ab_out,message,'COLL')
 713    call wrtout(std_out,  message,'COLL')
 714 
 715    read (tmp_unit,*,err=10,iomsg=errmsg) k11d,k22d,k33d
 716    write(message, '(a,f12.7,a,f12.7,a,f12.7)' )'                    k11d=',k11d,'; k22d=',k22d,'; k33d=',k33d
 717    call wrtout(ab_out,message,'COLL')
 718    call wrtout(std_out,  message,'COLL')
 719  end if
 720 
 721  if (lmax > 2) then
 722    read (tmp_unit,*,err=10,iomsg=errmsg) rrf,h11f
 723    write(message, '(a,f12.7,a,f12.7)' )' rrf =',rrf,'; h11f=',h11f
 724    call wrtout(ab_out,message,'COLL')
 725    call wrtout(std_out,  message,'COLL')
 726 
 727    read (tmp_unit,*,err=10,iomsg=errmsg) k11f
 728    write(message, '(a,f12.7)' )'                    k11f=',k11f
 729    call wrtout(ab_out,message,'COLL')
 730    call wrtout(std_out,  message,'COLL')
 731  end if
 732 
 733  if (abs(h11s)>1.d-08) nproj(1)=1
 734  if (abs(h22s)>1.d-08) nproj(1)=2
 735  if (abs(h33s)>1.d-08) nproj(1)=3
 736 
 737  if (abs(h11p)>1.d-08) then
 738    nproj(2)=1
 739    if (lmax<1) then
 740      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 741 &     ' psp3in : COMMENT -',ch10,&
 742 &     '  input lmax=',lmax,'  does not agree with input h11p=',h11p,ch10,&
 743 &     '  setting lmax to 1'
 744      call wrtout(ab_out,message,'COLL')
 745      call wrtout(std_out,  message,'COLL')
 746      lmax=1
 747    end if
 748  end if
 749 
 750  if (abs(h22p)>1.d-08) then
 751    nproj(2)=2
 752    if (lmax<1) then
 753      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 754 &     ' psp3in : COMMENT -',ch10,&
 755 &     '  input lmax=',lmax,' does not agree with input h22p=',h22p,ch10,&
 756 &     '  setting lmax to 1'
 757      call wrtout(ab_out,message,'COLL')
 758      call wrtout(std_out,  message,'COLL')
 759      lmax=1
 760    end if
 761  end if
 762 
 763  if (abs(h33p)>1.d-08) then
 764    nproj(2)=3
 765    if (lmax<1) then
 766      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 767 &     ' psp3in : COMMENT -',ch10,&
 768 &     '  input lmax=',lmax,' does not agree with input h33p=',h33p,ch10,&
 769 &     '  setting lmax to 1'
 770      call wrtout(ab_out,message,'COLL')
 771      call wrtout(std_out,  message,'COLL')
 772      lmax=1
 773    end if
 774  end if
 775 
 776  if (abs(h11d)>1.d-08) then
 777    nproj(3)=1
 778    if (lmax<2) then
 779      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 780 &     ' psp3in : COMMENT -',ch10,&
 781 &     '  input lmax=',lmax,'  does not agree with input h11d=',h11d,ch10,&
 782 &     '  setting lmax to 2'
 783      call wrtout(ab_out,message,'COLL')
 784      call wrtout(std_out,  message,'COLL')
 785      lmax=2
 786    end if
 787  end if
 788 
 789  if (abs(h22d)>1.d-08) then
 790    nproj(3)=2
 791    if (lmax<2) then
 792      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 793 &     ' psp3in : COMMENT -',ch10,&
 794 &     '  input lmax=',lmax,'  does not agree with input h22d=',h22d,ch10,&
 795 &     '  setting lmax to 2'
 796      call wrtout(ab_out,message,'COLL')
 797      call wrtout(std_out,  message,'COLL')
 798      lmax=2
 799    end if
 800  end if
 801 
 802  if (abs(h33d)>1.d-08) then
 803    nproj(3)=3
 804    if (lmax<2) then
 805      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 806 &     ' psp3in : COMMENT -',ch10,&
 807 &     '  input lmax=',lmax,' does not agree with input h33d=',h33d,ch10,&
 808 &     '  setting lmax to 2'
 809      call wrtout(ab_out,message,'COLL')
 810      call wrtout(std_out,  message,'COLL')
 811      lmax=2
 812    end if
 813  end if
 814 
 815  if (abs(h11f)>1.d-08) then
 816    nproj(4)=1
 817    if (lmax<3) then
 818      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 819 &     ' psp3in : COMMENT -',ch10,&
 820 &     '  input lmax=',lmax,' does not agree with input h11f=',h11f,ch10,&
 821 &     '  setting lmax to 3'
 822      call wrtout(ab_out,message,'COLL')
 823      call wrtout(std_out,  message,'COLL')
 824      lmax=3
 825    end if
 826  end if
 827 
 828  if(pspso/=0) then
 829 
 830    if (abs(k11p)>1.d-08) then
 831      nproj(psps%mpsang+1)=1
 832      if (lmax<1) then
 833        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 834 &       ' psp3in : COMMENT -',ch10,&
 835 &       '  input lmax=',lmax,'  does not agree with input k11p=',k11p,ch10,&
 836 &       '  setting lmax to 1'
 837        call wrtout(ab_out,message,'COLL')
 838        call wrtout(std_out,  message,'COLL')
 839        lmax=1
 840      end if
 841    end if
 842 
 843    if (abs(k22p)>1.d-08) then
 844      nproj(psps%mpsang+1)=2
 845      if (lmax<1) then
 846        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 847 &       ' psp3in : COMMENT -',ch10,&
 848 &       '  input lmax=',lmax,' does not agree with input k22p=',k22p,ch10,&
 849 &       '  setting lmax to 1'
 850        call wrtout(ab_out,message,'COLL')
 851        call wrtout(std_out,  message,'COLL')
 852        lmax=1
 853      end if
 854    end if
 855 
 856 
 857    if (abs(k33p)>1.d-08) then
 858      nproj(psps%mpsang+1)=3
 859      if (lmax<1) then
 860        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 861 &       ' psp3in : COMMENT -',ch10,&
 862 &       '  input lmax=',lmax,' does not agree with input k33p=',k33p,ch10,&
 863 &       '  setting lmax to 1'
 864        call wrtout(ab_out,message,'COLL')
 865        call wrtout(std_out,  message,'COLL')
 866        lmax=1
 867      end if
 868    end if
 869 
 870    if (abs(k11d)>1.d-08) then
 871      nproj(psps%mpsang+2)=1
 872      if (lmax<2) then
 873        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 874 &       ' psp3in : COMMENT -',ch10,&
 875 &       '  input lmax=',lmax,'  does not agree with input k11d=',k11d,ch10,&
 876 &       '  setting lmax to 2'
 877        call wrtout(ab_out,message,'COLL')
 878        call wrtout(std_out,  message,'COLL')
 879        lmax=2
 880      end if
 881    end if
 882 
 883    if (abs(k22d)>1.d-08) then
 884      nproj(psps%mpsang+2)=2
 885      if (lmax<2) then
 886        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 887 &       ' psp3in : COMMENT -',ch10,&
 888 &       '  input lmax=',lmax,'  does not agree with input k22d=',k22d,ch10,&
 889 &       '  setting lmax to 2'
 890        call wrtout(ab_out,message,'COLL')
 891        call wrtout(std_out,  message,'COLL')
 892        lmax=2
 893      end if
 894    end if
 895 
 896    if (abs(k33d)>1.d-08) then
 897      nproj(psps%mpsang+2)=3
 898      if (lmax<2) then
 899        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 900 &       ' psp3in : COMMENT -',ch10,&
 901 &       '  input lmax=',lmax,' does not agree with input k33d=',k33d,ch10,&
 902 &       '  setting lmax to 2'
 903        call wrtout(ab_out,message,'COLL')
 904        call wrtout(std_out,  message,'COLL')
 905        lmax=2
 906      end if
 907    end if
 908 
 909    if (abs(k11f)>1.d-08) then
 910      nproj(psps%mpsang+3)=1
 911      if (lmax<3) then
 912        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
 913 &       ' psp3in : COMMENT -',ch10,&
 914 &       '  input lmax=',lmax,' does not agree with input k11f=',k11f,ch10,&
 915 &       '  setting lmax to 3'
 916        call wrtout(ab_out,message,'COLL')
 917        call wrtout(std_out,  message,'COLL')
 918        lmax=3
 919      end if
 920    end if
 921 
 922  end if
 923 
 924 !Store the coefficients.
 925  psps%gth_params%set(ipsp)          = .true.
 926  psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, zero, zero /)
 927  psps%gth_params%psppar(1, :, ipsp) = (/ rrs,  h11s, h22s, h33s, zero, zero, zero /)
 928  psps%gth_params%psppar(2, :, ipsp) = (/ rrp,  h11p, h22p, h33p, zero, zero, zero /)
 929  psps%gth_params%psppar(3, :, ipsp) = (/ rrd,  h11d, h22d, h33d, zero, zero, zero /)
 930  psps%gth_params%psppar(4, :, ipsp) = (/ rrf,  h11f, zero, zero, zero, zero, zero /)
 931 
 932 !Store the k coefficients
 933  psps%gth_params%psp_k_par(1, :, ipsp) = (/ zero, zero, zero /)
 934  psps%gth_params%psp_k_par(2, :, ipsp) = (/ k11p, k22p, k33p /)
 935  psps%gth_params%psp_k_par(3, :, ipsp) = (/ k11d, k22d, k33d /)
 936  psps%gth_params%psp_k_par(4, :, ipsp) = (/ k11f, zero, zero /)
 937 
 938 !Additionnal wavelet parameters
 939  if (dtset%usewvl == 1) then
 940    call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
 941  end if
 942 
 943 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn
 944  nso=1;if(pspso/=0) nso=2
 945  index=0;iln=0;indlmn(:,:)=0
 946  do nn=1,nso
 947    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
 948      if (nproj(ipsang)>0) then
 949        ll=ipsang-(nn-1)*lmax-1
 950        do kk=1,nproj(ipsang)
 951          iln=iln+1
 952          do mm=1,2*ll*psps%useylm+1
 953            index=index+1
 954            indlmn(1,index)=ll
 955            indlmn(2,index)=mm-ll*psps%useylm-1
 956            indlmn(3,index)=kk
 957            indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
 958            indlmn(5,index)=iln
 959            indlmn(6,index)=nn
 960          end do
 961        end do
 962      end if
 963    end do
 964  end do
 965 
 966  ABI_ALLOCATE(dvlspl,(psps%mqgrid_ff))
 967 !First, the local potential --  compute on q grid and fit spline
 968  call psp2lo(cc1,cc2,cc3,cc4,dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,vlspl(:,1),rloc,.true.,yp1,ypn,zion)
 969  ABI_DEALLOCATE(dvlspl)
 970 
 971 !DEBUG
 972 !write(std_out,*)' psp3in : after psp2lo '
 973 !ENDDEBUG
 974 
 975 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
 976  ABI_ALLOCATE(work_space,(psps%mqgrid_ff))
 977  ABI_ALLOCATE(work_spl,(psps%mqgrid_ff))
 978  call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl)
 979  vlspl(:,2)=work_spl(:)
 980  ABI_DEALLOCATE(work_space)
 981  ABI_DEALLOCATE(work_spl)
 982 
 983 !Second, compute KB energies and form factors and fit splines
 984  ekb(:)=zero
 985 
 986 !Check if any nonlocal projectors are being used
 987  mproj=maxval(nproj)
 988  if (mproj>0) then
 989 
 990    ABI_ALLOCATE(ekb_sr,(psps%mpsang,mproj))
 991    ABI_ALLOCATE(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj))
 992    ABI_ALLOCATE(ekb_so,(psps%mpsang,mproj))
 993    ABI_ALLOCATE(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj))
 994 
 995    call psp3nl(ekb_sr,ffspl_sr,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
 996 &   h33d,h11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs)
 997    if(pspso/=0) then
 998      call psp3nl(ekb_so,ffspl_so,zero,zero,zero,k11p,k22p,k33p,k11d,&
 999 &     k22d,k33d,k11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs)
1000    end if
1001 
1002 
1003 !  Convert ekb and ffspl
1004    iln0=0
1005    do jj=1,psps%lmnmax
1006      iln=indlmn(5,jj)
1007      if (iln>iln0) then
1008        iln0=iln
1009        if (indlmn(6,jj)<=1) then
1010          ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj))
1011          ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj))
1012        else
1013          ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj))
1014          ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj))
1015        end if
1016      end if
1017    end do
1018 
1019    ABI_DEALLOCATE(ekb_sr)
1020    ABI_DEALLOCATE(ffspl_sr)
1021    ABI_DEALLOCATE(ekb_so)
1022    ABI_DEALLOCATE(ffspl_so)
1023  end if
1024 
1025  return
1026 
1027  ! Handle IO error
1028  10 continue
1029  MSG_ERROR(errmsg)
1030 
1031 end subroutine psp3in

m_psp_hgh/psp10in [ Functions ]

[ Top ] [ m_psp_hgh ] [ Functions ]

NAME

 psp10in

FUNCTION

 Initialize pspcod=10 pseudopotentials (formalism is the same as in HGH psps
 PRB58,3641(1998) [[cite:Hartwigsen1998]], but the full h and k matrices are read, allowing for using
 also subsequent developments such as Theor. Chem. Acc. 114, 145 (2005) [[cite:Dolg2005]]:
 continue to read the file, then compute the corresponding
 local and non-local potentials.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  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 characteristics of the psp
   if =3 : this input requires HFN characteristics of the psp
  ipsp=id in the array of the pseudo-potential.
  zion=nominal valence 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
  epsatm=$ (4\pi)\int_0^\infty [r^2 (V(r)+\frac{Zv}{r}) dr]$ (hartree)
  ffspl(mqgrid_ff,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
  vlspl(mqgrid_ff,2)=q^2 Vloc(q) and second derivatives from spline fit

SIDE EFFECTS

  Input/output
  lmax : at input =value of lmax mentioned at the second line of the psp file
    at output= 1
  psps <type(pseudopotential_type)>=at output, values depending on the read
                                    pseudo are set.
   | lmnmax(IN)=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(IN)=max. number of (l,n) components over all type of psps
   |           angular momentum of nonlocal pseudopotential
   | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
   | mpssoang(IN)= 1+maximum (spin*angular momentum) for nonlocal pseudopotentials
   | mqgrid_ff(IN)=dimension of q (or G) grid for arrays.
   | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
   | useylm(IN)=governs the way the nonlocal operator is to be applied:
   |            1=using Ylm, 0=using Legendre polynomials

PARENTS

      pspatm

CHILDREN

      psp10nl,psp2lo,spline,wrtout,wvl_descr_psp_fill

SOURCE

1452 subroutine psp10in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion)
1453 
1454 
1455 !This section has been created automatically by the script Abilint (TD).
1456 !Do not modify the following lines by hand.
1457 #undef ABI_FUNC
1458 #define ABI_FUNC 'psp10in'
1459 !End of the abilint section
1460 
1461  implicit none
1462 
1463 !Arguments ------------------------------------
1464 !scalars
1465  integer,intent(in) :: ipsp,pspso
1466  integer,intent(inout) :: lmax
1467  real(dp),intent(in) :: zion
1468  real(dp),intent(out) :: epsatm
1469  type(dataset_type),intent(in) :: dtset
1470  type(pseudopotential_type),intent(inout) :: psps
1471 !arrays
1472  integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang)
1473  real(dp),intent(out) :: ekb(psps%lnmax) !vz_i
1474  real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i
1475  real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2)
1476 
1477 !Local variables-------------------------------
1478 !scalars
1479  integer :: ii,iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nnonloc,nprl,nso
1480  real(dp) :: rloc,yp1,ypn
1481  character(len=500) :: message,errmsg
1482 !arrays
1483  integer,allocatable :: dummy_nproj(:)
1484  real(dp) :: cc(4)
1485  real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:)
1486  real(dp),allocatable :: ffspl_sr(:,:,:,:),hij(:,:,:),kij(:,:,:),rr(:)
1487  real(dp),allocatable :: work_space(:),work_spl(:)
1488 
1489 ! ***************************************************************************
1490 
1491 !Set various terms to 0 in case not defined below
1492 !HGH values
1493  rloc=zero ; cc(:)=zero
1494  nproj(1:psps%mpssoang)=0
1495 
1496 !Read and write different lines of the pseudopotential file
1497 
1498  read (tmp_unit,*,err=10,iomsg=errmsg) rloc,nn,(cc(jj),jj=1,nn)
1499  write(message, '(a,f12.7)' ) ' rloc=',rloc
1500  call wrtout(ab_out,message,'COLL')
1501  call wrtout(std_out,  message,'COLL')
1502  write(message, '(a,i1,a,4f12.7)' )' cc(1:',nn,')=',(cc(jj),jj=1,nn)
1503  call wrtout(ab_out,message,'COLL')
1504  call wrtout(std_out,  message,'COLL')
1505 
1506 !Read the number of the non-local projectors
1507  read (tmp_unit,*,err=10,iomsg=errmsg) nnonloc
1508  if (nnonloc/=lmax+1) then
1509    write(message, '(a,a,a,a,i5,a,i5,a,a,a,a,i5)' ) ch10,&
1510 &   ' psp10in : COMMENT -',ch10,&
1511 &   '  input lmax=',lmax,'  does not agree with input nnonloc=',nnonloc,ch10,&
1512 &   '  which has to be lmax+1.',ch10,&
1513 &   '  Setting lmax to ',nnonloc-1
1514    call wrtout(ab_out,message,'COLL')
1515    call wrtout(std_out,  message,'COLL')
1516    lmax=1
1517  end if
1518  ABI_ALLOCATE(rr,(0:lmax))
1519  ABI_ALLOCATE(hij,(0:lmax,3,3))
1520  ABI_ALLOCATE(kij,(0:lmax,3,3))
1521  rr(:)=zero; hij(:,:,:)=zero; kij(:,:,:)=zero
1522 
1523 !Read and echo the coefficients of non-local projectors
1524  prjloop: do ll=0,lmax
1525    read (tmp_unit,*,err=10,iomsg=errmsg) rr(ll),nprl,(hij(ll,1,jj),jj=1,nprl)
1526    do ii=2,nprl
1527      read (tmp_unit,*,err=10,iomsg=errmsg) (hij(ll,ii,jj),jj=ii,nprl)
1528    end do
1529    nproj(ll+1)=nprl
1530    write(message, '(a,i3,a,f12.7,2a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )&
1531 &   ' for angular momentum l =',ll,' r(l) =',rr(ll),ch10,&
1532 &   '   h11, h12, h13 =', (hij(ll,1,jj),jj=1,3),ch10,&
1533 &   '        h22, h23 =', (hij(ll,2,jj),jj=2,3),ch10,&
1534 &   '             h33 =', (hij(ll,3,jj),jj=3,3)
1535    call wrtout(ab_out,message,'COLL')
1536    call wrtout(std_out,  message,'COLL')
1537    if (ll==0) cycle
1538    do ii=1,nprl
1539      read (tmp_unit,*,err=10,iomsg=errmsg) (kij(ll,ii,jj),jj=ii,nprl)
1540    end do
1541    write(message, '(a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )&
1542 &   '   k11, k12, k13 =', (kij(ll,1,jj),jj=1,3),ch10,&
1543 &   '        k22, k23 =', (kij(ll,2,jj),jj=2,3),ch10,&
1544 &   '             k33 =', (kij(ll,3,jj),jj=3,3)
1545    call wrtout(ab_out,message,'COLL')
1546    call wrtout(std_out,  message,'COLL')
1547  end do prjloop
1548 
1549  if(pspso/=0) then
1550 
1551 !  MJV 10/2008: is this correct? For the normal HGH psp there are cases
1552 !  where there are more SO projectors than SR ones! e.g. Pb with 12 electrons.
1553    do ll=1,lmax
1554      nproj(psps%mpsang+ll)=nproj(ll+1)
1555    end do
1556 
1557  end if
1558 
1559 !Store the coefficients.
1560  psps%gth_params%set(ipsp)          = .true.
1561  psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc(1), cc(2), cc(3), cc(4), zero, zero /)
1562  do ii=1,4
1563    ll=ii-1
1564    if (ll>lmax) then
1565      psps%gth_params%psppar(ii,:,ipsp) = (/ zero, zero, zero, zero, zero, zero, zero /)
1566    else
1567      psps%gth_params%psppar(ii,:,ipsp) =&
1568 &     (/ rr(ll), hij(ll,1,1), hij(ll,2,2), hij(ll,3,3), hij(ll,1,2), hij(ll,1,3), hij(ll,2,3) /)
1569    end if
1570  end do
1571 
1572 !Additionnal wavelet parameters
1573  if (dtset%usewvl == 1) then
1574    call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
1575  end if
1576 
1577 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn
1578  nso=1;if(pspso/=0) nso=2
1579  index=0;iln=0;indlmn(:,:)=0
1580  do nn=1,nso
1581    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
1582      if (nproj(ipsang)>0) then
1583        ll=ipsang-(nn-1)*lmax-1
1584        do kk=1,nproj(ipsang)
1585          iln=iln+1
1586          do mm=1,2*ll*psps%useylm+1
1587            index=index+1
1588            indlmn(1,index)=ll
1589            indlmn(2,index)=mm-ll*psps%useylm-1
1590            indlmn(3,index)=kk
1591            indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
1592            indlmn(5,index)=iln
1593            indlmn(6,index)=nn
1594          end do
1595        end do
1596      end if
1597    end do
1598  end do
1599 
1600  ABI_ALLOCATE(dvlspl,(psps%mqgrid_ff))
1601 !First, the local potential --  compute on q grid and fit spline
1602  call psp2lo(cc(1),cc(2),cc(3),cc(4),dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,&
1603 & vlspl(:,1),rloc,.true.,yp1,ypn,zion)
1604  ABI_DEALLOCATE(dvlspl)
1605 
1606 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
1607  ABI_ALLOCATE(work_space,(psps%mqgrid_ff))
1608  ABI_ALLOCATE(work_spl,(psps%mqgrid_ff))
1609  call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl)
1610  vlspl(:,2)=work_spl(:)
1611  ABI_DEALLOCATE(work_space)
1612  ABI_DEALLOCATE(work_spl)
1613 
1614 !Second, compute KB energies and form factors and fit splines
1615  ekb(:)=zero
1616 
1617 !Check if any nonlocal projectors are being used
1618  mproj=maxval(nproj)
1619 
1620  if (mproj>0) then
1621 
1622    ABI_ALLOCATE(ekb_sr,(psps%mpsang,mproj))
1623    ABI_ALLOCATE(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj))
1624    ABI_ALLOCATE(ekb_so,(psps%mpsang,mproj))
1625    ABI_ALLOCATE(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj))
1626 
1627    call psp10nl(ekb_sr,ffspl_sr,hij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,&
1628 &   nproj,psps%qgrid_ff,rr)
1629    if(pspso/=0) then
1630      ABI_ALLOCATE(dummy_nproj,(psps%mpsang))
1631      dummy_nproj(1)=0
1632      do ll=1,lmax
1633        dummy_nproj(ll+1)=nproj(psps%mpsang+ll)
1634      end do
1635      call psp10nl(ekb_so,ffspl_so,kij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,&
1636 &     dummy_nproj,psps%qgrid_ff,rr)
1637      ABI_DEALLOCATE(dummy_nproj)
1638    end if
1639 
1640 !  Convert ekb and ffspl
1641    iln0=0
1642    do jj=1,psps%lmnmax
1643      iln=indlmn(5,jj)
1644      if (iln>iln0) then
1645        iln0=iln
1646        if (indlmn(6,jj)<=1) then
1647          ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj))
1648          ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj))
1649        else
1650          ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj))
1651          ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj))
1652        end if
1653      end if
1654    end do
1655 
1656    ABI_DEALLOCATE(ekb_sr)
1657    ABI_DEALLOCATE(ffspl_sr)
1658    ABI_DEALLOCATE(ekb_so)
1659    ABI_DEALLOCATE(ffspl_so)
1660  end if
1661 
1662  ABI_DEALLOCATE(rr)
1663  ABI_DEALLOCATE(hij)
1664  ABI_DEALLOCATE(kij)
1665 
1666  return
1667 
1668  ! Handle IO error
1669  10 continue
1670  MSG_ERROR(errmsg)
1671 
1672 end subroutine psp10in

m_psp_hgh/psp10nl [ Functions ]

[ Top ] [ m_psp_hgh ] [ Functions ]

NAME

 psp10nl

FUNCTION

 Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998).
 Uses Gaussians for fully nonlocal form, analytic expressions.

INPUTS

  hij(0:lmax,3,3)=factor defining strength of (max 3) projectors for each
   angular momentum channel l among 0, 1, ..., lmax
  lmax=maximum angular momentum
  mproj=maximum number of projectors in any channel
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=number of grid points for qgrid
  nproj(1:lmax+1)=number of projectors in any channel
  qgrid(mqgrid)=array of |G| values
  rr(0:lmax)=core radius for each 0<l<lmax channel (bohr)

OUTPUT

  ekb(mpsang,mproj)=Kleinman-Bylander energies
  ffspl(mqgrid,2,mpssang,mproj)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum and
   each projectors

PARENTS

      psp10in

CHILDREN

      spline,zhpev

SOURCE

1708 subroutine psp10nl(ekb,ffspl,hij,lmax,mproj,mpsang,mqgrid,nproj,qgrid,rr)
1709 
1710 
1711 !This section has been created automatically by the script Abilint (TD).
1712 !Do not modify the following lines by hand.
1713 #undef ABI_FUNC
1714 #define ABI_FUNC 'psp10nl'
1715 !End of the abilint section
1716 
1717  implicit none
1718 
1719 !Arguments ------------------------------------
1720 !scalars
1721  integer,intent(in) :: lmax,mproj,mpsang,mqgrid
1722 !arrays
1723  integer,intent(in) :: nproj(mpsang)
1724  real(dp),intent(in) :: hij(0:lmax,3,3),qgrid(mqgrid),rr(0:lmax)
1725  real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
1726 
1727 !Local variables-------------------------------
1728 !scalars
1729  integer :: info,ipack,iproj,iqgrid,jproj,ll,numproj
1730  real(dp) :: qmax,rrl
1731  character(len=500) :: message
1732  character :: jobz,uplo
1733 !arrays
1734  real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
1735  real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
1736 
1737 ! *************************************************************************
1738 
1739  ABI_ALLOCATE(ppspl,(mqgrid,2,mpsang,mproj))
1740  ABI_ALLOCATE(work,(mqgrid))
1741 
1742  qmax=qgrid(mqgrid)
1743  jobz='v'
1744  uplo='u'
1745  ekb(:,:)=zero
1746 
1747  lloop: do ll=0,lmax
1748    ap(:,:)=zero
1749    numproj=nproj(ll+1)
1750 
1751 !  Fill up the matrix in packed storage
1752    prjloop: do jproj=1,numproj
1753      priloop: do iproj=1,jproj
1754        ipack=iproj+(jproj-1)*jproj/2
1755        if(mod((jproj-1)*jproj,2)/=0) then
1756          MSG_ERROR("odd")
1757        end if
1758        ap(1,ipack)=hij(ll,iproj,jproj)
1759      end do priloop
1760    end do prjloop
1761 
1762    if(numproj/=0)then
1763 
1764      ABI_ALLOCATE(uu,(numproj,numproj))
1765      ABI_ALLOCATE(zz,(2,numproj,numproj))
1766 
1767      if (numproj > 1) then
1768        call ZHPEV(jobz,uplo,numproj,ap,ww,zz,numproj,work1,rwork1,info)
1769        uu(:,:)=zz(1,:,:)
1770      else
1771        ww(1)=hij(ll,1,1)
1772        uu(1,1)=one
1773      end if
1774 
1775 !    Initialization of ekb, and spline fitting
1776 
1777      if (ll==0) then ! s channel
1778 
1779        rrl=rr(0)
1780        do iproj=1,numproj
1781          ekb(1,iproj)=ww(iproj)*32.d0*(rrl**3)*(pi**2.5d0)/(4.d0*pi)**2
1782          if(iproj==1)then
1783            do iqgrid=1,mqgrid
1784              ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2)
1785            end do
1786            yp1j(1)=zero
1787            ypnj(1)=-(two_pi*rrl)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrl)**2)
1788          else if(iproj==2)then
1789            do iqgrid=1,mqgrid
1790              ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0)     &
1791 &             *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) &
1792 &             *( 3.d0-(two_pi*qgrid(iqgrid)*rrl)**2 )
1793            end do
1794            yp1j(2)=zero
1795            ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrl)**2*qmax &
1796 &           *exp(-0.5d0*(two_pi*qmax*rrl)**2) * (-5.d0+(two_pi*qmax*rrl)**2)
1797          else if(iproj==3)then
1798            do iqgrid=1,mqgrid
1799              ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
1800 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1801 &             (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrl)**2 + &
1802 &             (two_pi*qgrid(iqgrid)*rrl)**4)
1803            end do
1804            yp1j(3)=zero
1805            ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2) * &
1806 &           (two_pi*rrl)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrl)**2-(two_pi*qmax*rrl)**4)
1807          end if
1808          call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
1809 &         yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
1810        end do
1811 
1812      else if (ll==1) then ! p channel
1813 
1814        rrl=rr(1)
1815        do iproj=1,numproj
1816          ekb(2,iproj)=ww(iproj)*64.d0*(rrl**5)*(pi**2.5d0)/(4.d0*pi)**2
1817          if(iproj==1)then
1818            do iqgrid=1,mqgrid
1819              ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
1820 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))
1821            end do
1822            yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
1823            ypnj(1)=-two_pi*((two_pi*qmax*rrl)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)*&
1824 &           (1.0d0/sqrt(3.0d0))
1825          else if(iproj==2)then
1826            do iqgrid=1,mqgrid
1827              ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
1828 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1829 &             (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrl)**2)
1830            end do
1831            yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
1832            ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1833 &           (-8*(two_pi*qmax*rrl)**2 + (two_pi*qmax*rrl)**4 + 5.0d0)
1834          else if(iproj==3)then
1835            do iqgrid=1,mqgrid
1836              ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
1837 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1838 &             (two_pi*qgrid(iqgrid))*&
1839 &             (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrl)**2+(two_pi*qgrid(iqgrid)*rrl)**4)
1840            end do
1841            yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
1842            ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1843 &           (35.0d0-77.0d0*(two_pi*qmax*rrl)**2+19.0d0*(two_pi*qmax*rrl)**4 - &
1844 &           (two_pi*qmax*rrl)**6)
1845          end if
1846          call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
1847 &         yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
1848        end do
1849 
1850      else if (ll==2) then ! d channel
1851 
1852 !      If there is a third projector. Warning : only two projectors are allowed.
1853        if ( numproj>2 ) then
1854          write(message, '(3a)' )&
1855 &         ' only two d-projectors are allowed ',ch10,&
1856 &         ' Action: check your pseudopotential file.'
1857          MSG_ERROR(message)
1858        end if
1859 
1860        rrl=rr(2)
1861        do iproj=1,numproj
1862          ekb(3,iproj)=ww(iproj)*128.d0*(rrl**7)*(pi**2.5d0)/(4.d0*pi)**2
1863          if(iproj==1)then
1864            do iqgrid=1,mqgrid
1865              ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
1866 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))**2
1867            end do
1868            yp1j(1)=zero
1869            ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
1870 &           exp(-0.5d0*(two_pi*qmax*rrl)**2)*qmax*(2d0-(two_pi*qmax*rrl)**2)
1871          else if(iproj==2)then
1872            do iqgrid=1,mqgrid
1873              ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
1874 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1875 &             ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrl)**2)
1876            end do
1877            yp1j(2)=zero
1878            ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1879 &           qmax*(two_pi**2)*( (two_pi*qmax*rrl)**4 - 11.0d0*(two_pi*qmax*rrl)**2 + 14.0d0)
1880          end if
1881          call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
1882 &         yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
1883        end do
1884 
1885      else if (ll==3) then ! f channel
1886 
1887 !      If there is a second projector. Warning : only one projector is allowed.
1888        if ( numproj>1 ) then
1889          write(message, '(a,a,a)' )&
1890 &         'only one f-projector is allowed ',ch10,&
1891 &         'Action: check your pseudopotential file.'
1892          MSG_ERROR(message)
1893        end if
1894 
1895        rrl=rr(3)
1896        ekb(4,1)=ww(1)*(256.0d0/105.0d0)*(rrl**9)*(pi**2.5d0)/(4.d0*pi)**2
1897        do iqgrid=1,mqgrid
1898          ppspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
1899 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2)
1900        end do
1901 !      Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
1902        yp1j(1)=zero
1903        ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrl)**2)*&
1904 &       (3.0d0-(two_pi*qmax*rrl)**2)
1905 !      Fit spline to get second derivatives by spline fit
1906        call spline(qgrid,ppspl(:,1,4,1),mqgrid,&
1907 &       yp1j(1),ypnj(1),ppspl(:,2,4,1))
1908 
1909      else
1910        MSG_ERROR("lmax>3?")
1911      end if
1912 
1913 !    Linear combination using the eigenvectors
1914      ffspl(:,:,ll+1,:)=zero
1915      do jproj=1,numproj
1916        do iproj=1,numproj
1917          do iqgrid=1,mqgrid
1918            ffspl(iqgrid,1:2,ll+1,jproj)=ffspl(iqgrid,1:2,ll+1,jproj) &
1919 &           +uu(iproj,jproj)*ppspl(iqgrid,1:2,ll+1,iproj)
1920          end do
1921        end do
1922      end do
1923 
1924      ABI_DEALLOCATE(uu)
1925      ABI_DEALLOCATE(zz)
1926 
1927 !    End condition on numproj(/=0)
1928    end if
1929 
1930  end do lloop
1931 
1932  ABI_DEALLOCATE(ppspl)
1933  ABI_DEALLOCATE(work)
1934 
1935 end subroutine psp10nl

m_psp_hgh/psp2in [ Functions ]

[ Top ] [ m_psp_hgh ] [ Functions ]

NAME

 psp2in

FUNCTION

 Initialize pspcod=2 pseudopotentials (GTH format):
 continue to read the file, then compute the corresponding
 local and non-local potentials.

INPUTS

  dtset <type(dataset_type)>=all input variables in this dataset
  ipsp=id in the array of the pseudo-potential.
  lmax=value of lmax mentioned at the second line of the psp file
  zion=nominal valence 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)+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
  vlspl(mqgrid_vl,2)=q^2 Vloc(q) and second derivatives from spline fit
  dvlspl(mqgrid_vl,2)=dVloc(r)/dr and second derivatives from spline fit (only
                      allocated if vlspl_recipSpace is false.

SIDE EFFECTS

  Input/output
  lmax : at input =value of lmax mentioned at the second line of the psp file
    at output= 1
  psps <type(pseudopotential_type)>=at output, values depending on the read
                                    pseudo are set.
   | lmnmax(IN)=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(IN)=max. number of (l,n) components over all type of psps
   |           angular momentum of nonlocal pseudopotential
   | mpsang(IN)= 1+maximum angular momentum for nonlocal pseudopotentials
   | mqgrid_ff(IN)=dimension of q (or G) grid for nl form factors (array ffspl)
   | mqgrid_vl(IN)=dimension of q (or G) grid or r grid (if vlspl_recipSpace = .false.)
   | qgrid_ff(mqgrid_ff)(IN)=values of q on grid from 0 to qmax (bohr^-1) for nl form factors
   | qgrid_vl(mqgrid_vl)(IN)=values of q on grid from 0 to qmax (bohr^-1) for Vloc
   |                         if vlspl_recipSpace is .true. else values of r on grid from
   |                         0 to 2pi / qmax * mqgrid_ff (bohr).
   | useylm(IN)=governs the way the nonlocal operator is to be applied:
   |            1=using Ylm, 0=using Legendre polynomials
   | vlspl_recipSpace(IN)=.true. if pseudo are expressed in reciprocal space.
   | gth_params(OUT)=store GTH coefficients and parameters.

PARENTS

      pspatm

CHILDREN

      psp2lo,psp2nl,spline,wrtout,wvl_descr_psp_fill

SOURCE

118 subroutine psp2in(dtset,ekb,epsatm,ffspl,indlmn,ipsp,lmax,nproj,psps,vlspl,dvlspl,zion)
119 
120 
121 !This section has been created automatically by the script Abilint (TD).
122 !Do not modify the following lines by hand.
123 #undef ABI_FUNC
124 #define ABI_FUNC 'psp2in'
125 !End of the abilint section
126 
127  implicit none
128 
129 !Arguments ------------------------------------
130 !scalars
131  integer,intent(in) :: ipsp,lmax
132  real(dp),intent(in) :: zion
133  real(dp),intent(out) :: epsatm
134  type(dataset_type),intent(in) :: dtset
135  type(pseudopotential_type),intent(inout) :: psps
136 !arrays
137  integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpsang)
138  real(dp),intent(out) :: dvlspl(psps%mqgrid_vl,2),ekb(psps%lnmax)
139  real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i
140  real(dp),intent(out) :: vlspl(psps%mqgrid_vl,2)
141 
142 !Local variables-------------------------------
143 !scalars
144  integer :: iln,index,ipsang,kk,ll,mm
145  real(dp) :: cc1,cc2,cc3,cc4,h1p,h1s,h2s,rloc,rrp,rrs
146  real(dp) :: yp1,ypn
147  character(len=500) :: message,errmsg
148 !arrays
149  real(dp),allocatable :: work_space(:),work_spl(:)
150  real(dp),allocatable :: dvloc(:)
151 
152 ! ***************************************************************************
153 
154 !Set various terms to 0 in case not defined below
155 !GTH values
156  rloc=0.d0
157  cc1=0.d0
158  cc2=0.d0
159  cc3=0.d0
160  cc4=0.d0
161  rrs=0.d0
162  h1s=0.d0
163  h2s=0.d0
164  rrp=0.d0
165  h1p=0.d0
166  nproj(1:psps%mpsang)=0
167 
168 !Read and write different lines of the pseudopotential file
169  read (tmp_unit,*, err=10, iomsg=errmsg) rloc,cc1,cc2,cc3,cc4
170  write(message, '(a,f12.7)' ) ' rloc=',rloc
171  call wrtout(ab_out,message,'COLL')
172  call wrtout(std_out,  message,'COLL')
173  write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )'  cc1=',cc1,'; cc2=',cc2,'; cc3=',cc3,'; cc4=',cc4
174  call wrtout(ab_out,message,'COLL')
175  call wrtout(std_out,  message,'COLL')
176 
177  read (tmp_unit,*, err=10, iomsg=errmsg) rrs,h1s,h2s
178  write(message, '(a,f12.7,a,f12.7,a,f12.7)' )'  rrs=',rrs,'; h1s=',h1s,'; h2s=',h2s
179  call wrtout(ab_out,message,'COLL')
180  call wrtout(std_out,  message,'COLL')
181 
182  read (tmp_unit,*, err=10, iomsg=errmsg) rrp,h1p
183  write(message, '(a,f12.7,a,f12.7)' )'  rrp=',rrp,'; h1p=',h1p
184  call wrtout(ab_out,message,'COLL')
185  call wrtout(std_out,  message,'COLL')
186 
187 !Store the coefficients.
188  psps%gth_params%set(ipsp)          = .true.
189  psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, 0.d0, 0.d0 /)
190  psps%gth_params%psppar(1, :, ipsp) = (/ rrs,  h1s, h2s, 0.d0, 0.d0, 0.d0, 0.d0 /)
191  psps%gth_params%psppar(2, :, ipsp) = (/ rrp,  h1p, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0 /)
192  if (dtset%usewvl == 1) then
193    call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
194  end if
195 
196  if (abs(h1s)>1.d-08) nproj(1)=1
197  if (abs(h2s)>1.d-08) nproj(1)=2
198 
199  if (abs(h1p)>1.d-08) then
200    if(psps%mpsang<2)then
201      write(message, '(a,es12.4,a,a,a,i2,a)' )&
202 &     'With non-zero h1p (=',h1p,'), mpsang should be at least 2,',ch10,&
203 &     'while mpsang=',psps%mpsang,'.'
204      MSG_ERROR(message)
205    end if
206    nproj(2)=1
207    if (lmax<1) then
208      write(message, '(a,i5,a,e12.4,a,a,a,a)' )&
209 &     'Input lmax=',lmax,' disagree with input h1p=',h1p,'.',&
210 &     'Your pseudopotential is incoherent.',ch10,&
211 &     'Action: correct your pseudopotential file.'
212      MSG_ERROR(message)
213    end if
214  end if
215 
216 !Initialize array indlmn array giving l,m,n,lm,ln,s for i=lmn
217  index=0;iln=0;indlmn(:,:)=0
218  do ipsang=1,lmax+1
219    if(nproj(ipsang)>0)then
220      ll=ipsang-1
221      do kk=1,nproj(ipsang)
222        iln=iln+1
223        do mm=1,2*ll*psps%useylm+1
224          index=index+1
225          indlmn(1,index)=ll
226          indlmn(2,index)=mm-ll*psps%useylm-1
227          indlmn(3,index)=kk
228          indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
229          indlmn(5,index)=iln
230          indlmn(6,index)=1
231        end do
232      end do
233    end if
234  end do
235 
236 !First, the local potential --
237 !compute q^2V(q) or V(r)
238 !MJV NOTE: psp2lo should never be called with dvspl unallocated, which
239 !is possible unless .not.psps%vlspl_recipSpace
240  ABI_ALLOCATE(dvloc,(psps%mqgrid_vl))
241  call psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,psps%mqgrid_vl,psps%qgrid_vl,&
242 & vlspl(:,1),rloc,psps%vlspl_recipSpace,yp1,ypn,zion)
243 
244 !Fit spline to (q^2)V(q) or V(r)
245  ABI_ALLOCATE(work_space,(psps%mqgrid_vl))
246  ABI_ALLOCATE(work_spl,(psps%mqgrid_vl))
247  call spline (psps%qgrid_vl,vlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl)
248  vlspl(:,2)=work_spl(:)
249  if (.not.psps%vlspl_recipSpace) then
250    dvlspl(:,1) = dvloc
251    call spline (psps%qgrid_vl,dvlspl(:,1),psps%mqgrid_vl,yp1,ypn,work_spl)
252    dvlspl(:,2)=work_spl(:)
253  end if
254 
255  ABI_DEALLOCATE(work_space)
256  ABI_DEALLOCATE(work_spl)
257  ABI_DEALLOCATE(dvloc)
258 
259 
260 !Second, compute KB energies and form factors and fit splines
261  ekb(:)=0.0d0
262 !First check if any nonlocal projectors are being used
263  if (maxval(nproj(1:lmax+1))>0) then
264    call psp2nl(ekb,ffspl,h1p,h1s,h2s,psps%lnmax,psps%mqgrid_ff,psps%qgrid_ff,rrp,rrs)
265  end if
266 
267  return
268 
269  ! Handle IO error
270  10 continue
271  MSG_ERROR(errmsg)
272 
273 end subroutine psp2in

m_psp_hgh/psp2nl [ Functions ]

[ Top ] [ m_psp_hgh ] [ Functions ]

NAME

 psp2nl

FUNCTION

 Goedecker-Teter-Hutter nonlocal pseudopotential (from preprint of 1996).
 Uses Gaussians for fully nonlocal form, analytic expressions.

INPUTS

  h1p=factor defining strength of 1st projector for l=1 channel
  h1s=factor defining strength of 1st projector for l=0 channel
  h2s=factor defining strength of 2nd projector for l=0 channel
  lnmax=max. number of (l,n) components over all type of psps
  mqgrid=number of grid points for qgrid
  qgrid(mqgrid)=array of |G| values
  rrp=core radius for p channel (bohr)
  rrs=core radius for s channel (bohr)

OUTPUT

  ekb(lnmax)=Kleinman-Bylander energy
  ffspl(mqgrid,2,lnmax)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum
   and each projector

PARENTS

      psp2in

CHILDREN

      spline

SOURCE

308 subroutine psp2nl(ekb,ffspl,h1p,h1s,h2s,lnmax,mqgrid,qgrid,rrp,rrs)
309 
310 
311 !This section has been created automatically by the script Abilint (TD).
312 !Do not modify the following lines by hand.
313 #undef ABI_FUNC
314 #define ABI_FUNC 'psp2nl'
315 !End of the abilint section
316 
317  implicit none
318 
319 !Arguments ------------------------------------
320 !scalars
321  integer,intent(in) :: lnmax,mqgrid
322  real(dp),intent(in) :: h1p,h1s,h2s,rrp,rrs
323 !arrays
324  real(dp),intent(in) :: qgrid(mqgrid)
325  real(dp),intent(inout) :: ekb(lnmax),ffspl(mqgrid,2,lnmax) !vz_i
326 
327 !Local variables-------------------------------
328 !scalars
329  integer :: iln,iqgrid
330  real(dp) :: qmax,yp1,ypn
331 !arrays
332  real(dp),allocatable :: work(:)
333 
334 ! *************************************************************************
335 
336  ABI_ALLOCATE(work,(mqgrid))
337 
338 !Kleinman-Bylander energies ekb were set to zero in calling program
339 
340 !Compute KB energies
341  iln=0
342  if (abs(h1s)>1.d-12) then
343    iln=iln+1
344    ekb(iln)=h1s*32.d0*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2)
345  end if
346  if (abs(h2s)>1.d-12) then
347    iln=iln+1
348    ekb(iln) =h2s*(128.d0/15.d0)*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2)
349  end if
350  if (abs(h1p)>1.d-12) then
351    iln=iln+1
352    ekb(iln)=h1p*(64.d0/3.d0)*rrp**5*(pi**(2.5d0)/(4.d0*pi)**2)
353  end if
354 
355 !Compute KB form factor
356  iln=0
357 
358 !l=0 first projector
359  if (abs(h1s)>1.d-12) then
360    iln=iln+1
361    do iqgrid=1,mqgrid
362      ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
363    end do
364 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
365    yp1=0.d0
366    qmax=qgrid(mqgrid)
367    ypn=-4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2)
368 !  Fit spline to get second derivatives by spline fit
369    call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
370 !  else
371 !  or else put first projector nonlocal correction at l=0 to 0
372 !  ffspl(:,:,iln)=0.0d0
373  end if
374 
375 !l=0 second projector
376  if (abs(h2s)>1.d-12) then
377    iln=iln+1
378    do iqgrid=1,mqgrid
379      ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
380 &     (3.d0-(two_pi*qgrid(iqgrid)*rrs)**2)
381    end do
382 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
383    yp1=0.d0
384    qmax=qgrid(mqgrid)
385    ypn=4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
386 &   (-5.d0+(two_pi*qmax*rrs)**2)
387 !  Fit spline to get second derivatives by spline fit
388    call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
389 !  else if(mproj>=2)then
390 !  or else put second projector nonlocal correction at l=0 to 0
391 !  ffspl(:,:,iln)=0.0d0
392  end if
393 
394 !l=1 first projector
395  if (abs(h1p)>1.d-12) then
396    iln=iln+1
397    do iqgrid=1,mqgrid
398      ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
399 &     (two_pi*qgrid(iqgrid))
400    end do
401 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
402    yp1=two_pi
403    qmax=qgrid(mqgrid)
404    ypn=-two_pi*((two_pi*qmax*rrp)**2-1.d0) * exp(-0.5d0*(two_pi*qmax*rrp)**2)
405 !  Fit spline to get second derivatives by spline fit
406    call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
407 !  else if(mpsang>=2)then
408 !  or else put first projector l=1 nonlocal correction to 0
409 !  ffspl(:,:,iln)=0.0d0
410  end if
411 
412  ABI_DEALLOCATE(work)
413 
414 end subroutine psp2nl

m_psp_hgh/psp3nl [ Functions ]

[ Top ] [ m_psp_hgh ] [ Functions ]

NAME

 psp3nl

FUNCTION

 Hartwigsen-Goedecker-Hutter nonlocal pseudopotential (from preprint of 1998).
 Uses Gaussians for fully nonlocal form, analytic expressions.

INPUTS

  h11s=factor defining strength of 1st projector for l=0 channel
  h22s=factor defining strength of 2nd projector for l=0 channel
  h33s=factor defining strength of 3rd projector for l=0 channel
  h11p=factor defining strength of 1st projector for l=1 channel
  h22p=factor defining strength of 2nd projector for l=1 channel
  h33p=factor defining strength of 2nd projector for l=1 channel
  h11d=factor defining strength of 1st projector for l=2 channel
  h22d=factor defining strength of 2nd projector for l=2 channel
  h33d=factor defining strength of 2nd projector for l=2 channel
  h11f=factor defining strength of 1st projector for l=3 channel
  mproj=maximum number of projectors in any channel
  mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
  mqgrid=number of grid points for qgrid
  qgrid(mqgrid)=array of |G| values
  rrd=core radius for d channel (bohr)
  rrf=core radius for f channel (bohr)
  rrp=core radius for p channel (bohr)
  rrs=core radius for s channel (bohr)

OUTPUT

  ekb(mpsang,mproj)=Kleinman-Bylander energies
  ffspl(mqgrid,2,mpssang,mproj)=Kleinman-Bylander form factor f_l(q) and
   second derivative from spline fit for each angular momentum and
   each projectors

PARENTS

      psp3in

CHILDREN

      spline,zhpev

SOURCE

1076 subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
1077 &                  h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs)
1078 
1079 
1080 !This section has been created automatically by the script Abilint (TD).
1081 !Do not modify the following lines by hand.
1082 #undef ABI_FUNC
1083 #define ABI_FUNC 'psp3nl'
1084 !End of the abilint section
1085 
1086  implicit none
1087 
1088 !Arguments ------------------------------------
1089 !scalars
1090  integer,intent(in) :: mproj,mpsang,mqgrid
1091  real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd
1092  real(dp),intent(in) :: rrf,rrp,rrs
1093 !arrays
1094  real(dp),intent(in) :: qgrid(mqgrid)
1095  real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
1096 
1097 !Local variables-------------------------------
1098 !scalars
1099  integer :: info,iproj,iqgrid,ldz,mu,nproj,nu
1100  real(dp) :: qmax
1101  character(len=500) :: message
1102  character :: jobz,uplo
1103 !arrays
1104  real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
1105  real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
1106 
1107 ! *************************************************************************
1108 
1109  ABI_ALLOCATE(ppspl,(mqgrid,2,mpsang,mproj))
1110  ABI_ALLOCATE(work,(mqgrid))
1111 
1112  qmax=qgrid(mqgrid)
1113  jobz='v'
1114  uplo='u'
1115  ekb(:,:)=0.0d0
1116 
1117 !---------------------------------------------------------------
1118 !Treat s channel
1119 
1120  nproj=0
1121  ap(:,:)=0.0d0
1122 !If there is at least one s-projector
1123  if  ( abs(h11s) >= 1.0d-8 ) then
1124    nproj=1 ; ldz=1 ; ap(1,1)=h11s
1125  end if
1126  nproj=1
1127 !If there is a second projector
1128  if  ( abs(h22s) >= 1.0d-8 ) then
1129    nproj=2 ; ldz=2 ; ap(1,3)=h22s
1130    ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s
1131  end if
1132 !If there is a third projector
1133  if ( abs(h33s) >= 1.0d-8 ) then
1134    nproj=3 ; ldz=3 ; ap(1,6)=h33s
1135    ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s
1136    ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s
1137  end if
1138 
1139  if(nproj/=0)then
1140 
1141    ABI_ALLOCATE(uu,(nproj,nproj))
1142    ABI_ALLOCATE(zz,(2,nproj,nproj))
1143 
1144    if (nproj > 1) then
1145      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1146      uu(:,:)=zz(1,:,:)
1147    else
1148      ww(1)=h11s
1149      uu(1,1)=1.0d0
1150    end if
1151 
1152 !  Initialization of ekb, and spline fitting
1153    do iproj=1,nproj
1154      ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2
1155      if(iproj==1)then
1156        do iqgrid=1,mqgrid
1157          ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
1158        end do
1159        yp1j(1)=0.d0
1160        ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2)
1161      else if(iproj==2)then
1162        do iqgrid=1,mqgrid
1163          ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0)     &
1164 &         *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) &
1165 &         *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 )
1166        end do
1167        yp1j(2)=0.0d0
1168        ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax &
1169 &       *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2)
1170      else if(iproj==3)then
1171        do iqgrid=1,mqgrid
1172          ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
1173 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
1174 &         (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + &
1175 &         (two_pi*qgrid(iqgrid)*rrs)**4)
1176        end do
1177        yp1j(3)=0.0d0
1178        ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
1179 &       (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4)
1180      end if
1181      call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
1182 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
1183    end do
1184 
1185 !  Linear combination using the eigenvectors
1186    ffspl(:,:,1,:)=0.0d0
1187    do mu=1,nproj
1188      do nu=1,nproj
1189        do iqgrid=1,mqgrid
1190          ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) &
1191 &         +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu)
1192        end do
1193      end do
1194    end do
1195 
1196    ABI_DEALLOCATE(uu)
1197    ABI_DEALLOCATE(zz)
1198  end if !  End condition on nproj(/=0)
1199 
1200 !--------------------------------------------------------------------
1201 !Now treat p channel
1202 
1203  nproj=0
1204  ap(:,:)=0.0d0
1205 !If there is at least one projector
1206  if  ( abs(h11p) >= 1.0d-8 ) then
1207    nproj=1 ; ldz=1 ; ap(1,1)=h11p
1208  end if
1209 !If there is a second projector
1210  if  ( abs(h22p) >= 1.0d-8 ) then
1211    nproj=2 ; ldz=2 ; ap(1,3)=h22p
1212    ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p
1213  end if
1214 !If there is a third projector
1215  if ( abs(h33p) >= 1.0d-8 ) then
1216    nproj=3 ; ldz=3 ; ap(1,6)=h33p
1217    ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p
1218    ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p
1219  end if
1220 
1221  if(nproj/=0)then
1222 
1223    ABI_ALLOCATE(uu,(nproj,nproj))
1224    ABI_ALLOCATE(zz,(2,nproj,nproj))
1225 
1226    if (nproj > 1) then
1227      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1228      uu(:,:)=zz(1,:,:)
1229    else
1230      ww(1)=h11p
1231      uu(1,1)=1.0d0
1232    end if
1233 
1234 !  Initialization of ekb, and spline fitting
1235    do iproj=1,nproj
1236      ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2
1237      if(iproj==1)then
1238        do iqgrid=1,mqgrid
1239          ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
1240 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid))
1241        end do
1242        yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
1243        ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*&
1244 &       (1.0d0/sqrt(3.0d0))
1245      else if(iproj==2)then
1246        do iqgrid=1,mqgrid
1247          ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
1248 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
1249 &         (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2)
1250        end do
1251        yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
1252        ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
1253 &       (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0)
1254      else if(iproj==3)then
1255        do iqgrid=1,mqgrid
1256          ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
1257 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
1258 &         (two_pi*qgrid(iqgrid))*&
1259 &         (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4)
1260        end do
1261        yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
1262        ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
1263 &       (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - &
1264 &       (two_pi*qmax*rrp)**6)
1265      end if
1266      call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
1267 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
1268    end do
1269 
1270 !  Linear combination using the eigenvectors
1271    ffspl(:,:,2,:)=0.0d0
1272    do mu=1,nproj
1273      do nu=1,nproj
1274        do iqgrid=1,mqgrid
1275          ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) &
1276 &         +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu)
1277        end do
1278      end do
1279    end do
1280 
1281    ABI_DEALLOCATE(uu)
1282    ABI_DEALLOCATE(zz)
1283  end if !  End condition on nproj(/=0)
1284 
1285 !-----------------------------------------------------------------------
1286 !Now treat d channel.
1287 
1288  nproj=0
1289  ap(:,:)=0.0d0
1290 !If there is at least one projector
1291  if  ( abs(h11d) >= 1.0d-8 ) then
1292    nproj=1 ; ldz=1 ; ap(1,1)=h11d
1293  end if
1294 !If there is a second projector
1295  if  ( abs(h22d) >= 1.0d-8 ) then
1296    nproj=2 ; ldz=2 ; ap(1,3)=h22d
1297    ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d
1298  end if
1299 !If there is a third projector. Warning : only two projectors are allowed.
1300  if ( abs(h33d) >= 1.0d-8 ) then
1301    write(message, '(a,a,a)' )&
1302 &   '  only two d-projectors are allowed ',ch10,&
1303 &   '  Action: check your pseudopotential file.'
1304    MSG_ERROR(message)
1305 !  nproj=3 ; ldz=3 ; ap(1,6)=h33d
1306 !  ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d
1307 !  ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d
1308  end if
1309 
1310  if(nproj/=0)then
1311 
1312    ABI_ALLOCATE(uu,(nproj,nproj))
1313    ABI_ALLOCATE(zz,(2,nproj,nproj))
1314 
1315    if (nproj > 1) then
1316      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1317      uu(:,:)=zz(1,:,:)
1318    else
1319      ww(1)=h11d
1320      uu(1,1)=1.0d0
1321    end if
1322 
1323 !  Initialization of ekb, and spline fitting
1324    do iproj=1,nproj
1325      ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2
1326      if(iproj==1)then
1327        do iqgrid=1,mqgrid
1328          ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
1329 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2
1330        end do
1331        yp1j(1)=0.0d0
1332        ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
1333 &       exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2)
1334      else if(iproj==2)then
1335        do iqgrid=1,mqgrid
1336          ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
1337 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * &
1338 &         ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2)
1339        end do
1340        yp1j(2)=0.0d0
1341        ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* &
1342 &       qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0)
1343      end if
1344      call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
1345 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
1346    end do
1347 
1348 !  Linear combination using the eigenvectors
1349    ffspl(:,:,3,:)=0.0d0
1350    do mu=1,nproj
1351      do nu=1,nproj
1352        do iqgrid=1,mqgrid
1353          ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) &
1354 &         +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu)
1355        end do
1356      end do
1357    end do
1358 
1359    ABI_DEALLOCATE(uu)
1360    ABI_DEALLOCATE(zz)
1361  end if !  End condition on nproj(/=0)
1362 
1363 !-----------------------------------------------------------------------
1364 !Treat now f channel (max one projector ! - so do not use ppspl)
1365 
1366 !l=3 first projector
1367  if (abs(h11f)>1.d-12) then
1368    ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2
1369    do iqgrid=1,mqgrid
1370      ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
1371 &     exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2)
1372    end do
1373 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
1374    yp1j(1)=0d0
1375    ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*&
1376 &   (3.0d0-(two_pi*qmax*rrf)**2)
1377 !  Fit spline to get second derivatives by spline fit
1378    call spline(qgrid,ffspl(:,1,4,1),mqgrid,&
1379 &   yp1j(1),ypnj(1),ffspl(:,2,4,1))
1380  end if
1381 
1382 !-----------------------------------------------------------------------
1383 
1384  ABI_DEALLOCATE(ppspl)
1385  ABI_DEALLOCATE(work)
1386 
1387 end subroutine psp3nl