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-2022 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 .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module m_psp_hgh
23 
24  use defs_basis
25  use m_splines
26  use m_abicore
27  use m_errors
28  use m_dtset
29 
30  use defs_datatypes, only : pseudopotential_type
31  use m_special_funcs,  only : abi_derfc
32 #if defined HAVE_BIGDFT
33  use BigDFT_API, only: atomic_info
34 #endif
35  use m_wvl_descr_psp,  only : wvl_descr_psp_fill
36 
37  implicit none
38 
39  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).

SOURCE

420 subroutine psp2lo(cc1,cc2,cc3,cc4,dvloc,epsatm,mqgrid,qgrid,q2vq,&
421 &  rloc,vlspl_recipSpace,yp1,ypn,zion)
422 
423 !Arguments ------------------------------------
424 !scalars
425  integer,intent(in) :: mqgrid
426  real(dp),intent(in) :: cc1,cc2,cc3,cc4,rloc,zion
427  real(dp),intent(out) :: epsatm,yp1,ypn
428  logical,intent(in) :: vlspl_recipSpace
429 !arrays
430  real(dp),intent(in) :: qgrid(mqgrid)
431  real(dp),intent(out) :: dvloc(mqgrid),q2vq(mqgrid)
432 
433 !Local variables-------------------------------
434 !scalars
435  integer :: iqgrid
436  real(dp) :: erfValue,gaussValue,polyValue,qmax,rq,rq2
437  character(len=500) :: message
438 
439 ! *************************************************************************
440 
441 !Compute epsatm = lim(q->0) [Vloc(q) + zion/(Pi*q^2)]
442  epsatm=2.d0*pi*rloc**2*zion+(2.d0*pi)**(1.5d0)*rloc**3*&
443 & (cc1+3.d0*cc2+15.d0*cc3+105.d0*cc4)
444 
445 !If vlspl_recipSpace is .true., we compute V(q)*q^2 in reciprocal space,
446 !else we compute V(r) in real space.
447  if (vlspl_recipSpace) then
448    write(message, '(a)' ) '-  Local part computed in reciprocal space.'
449    call wrtout(ab_out,message,'COLL')
450    call wrtout(std_out,  message,'COLL')
451 
452 !  d(q^2*V(q))/d(q) at q=0 and q=qmax
453    qmax=qgrid(mqgrid)
454    rq2=(2.d0*pi*qmax*rloc)**2
455    yp1=0.d0
456    ypn= (2.d0*pi*qmax*rloc**2)*exp(-0.5d0*rq2)* &
457 &   (2.d0*zion + sqrt(2.d0*pi)*rloc*&
458 &   (cc1*(2.d0-rq2) + cc2*(6.d0-7.d0*rq2+rq2**2) +&
459 &   cc3*(30.d0-55.d0*rq2+16.d0*rq2**2-rq2**3) +&
460 &   cc4*(210.d0-525.d0*rq2+231.d0*rq2**2-29.d0*rq2**3+rq2**4)))
461 !  ypn has been tested against Maple-derived expression.
462 
463 !  Compute q^2*vloc(q) on uniform grid
464    do iqgrid=1,mqgrid
465      rq2=(2.d0*pi*qgrid(iqgrid)*rloc)**2
466      q2vq(iqgrid)=exp(-0.5d0*rq2)*(-zion/pi+rq2*(rloc/sqrt(2.d0*pi)) *&
467 &     ( cc1 + cc2*(3.d0-rq2) + cc3*(15.d0-10.d0*rq2+rq2**2) +&
468 &     cc4*(105.d0-rq2*(105.d0-rq2*(21.d0-rq2)))  ))
469    end do
470  else
471    write(message, '(a)' ) '-  Local part computed in real space.'
472    call wrtout(ab_out,message,'COLL')
473    call wrtout(std_out,  message,'COLL')
474 
475 !  Compute derivatives for splines computations
476    yp1 = 0.d0
477    rq2 = (qgrid(mqgrid) / rloc) ** 2
478    erfValue = abi_derfc(sqrt(0.5d0 * rq2))
479    ypn = - 2.0d0 * zion / sqrt(2.d0 * pi) / qgrid(mqgrid) / rloc
480    ypn = ypn - rq2 * (cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3) / qgrid(mqgrid)
481    ypn = ypn + (2.d0 * cc2 * rq2 + 4.d0 * cc3 * rq2 ** 2 + 6.d0 * cc4 * rq2 ** 3) / qgrid(mqgrid)
482    ypn = ypn * exp(-0.5d0 * rq2)
483    ypn = ypn + zion / qgrid(mqgrid) ** 2 * erfValue
484 !  Note that ypn has been calculated on a full-proof a4 paper sheet.
485 
486 !  Compute local potential and its first derivatives.
487    do iqgrid = 1, mqgrid, 1
488      rq2 = (qgrid(iqgrid) / rloc) ** 2
489 !    Compute erf() part
490 !    Case r = 0
491      gaussValue = exp(-0.5d0 * rq2)
492      if (qgrid(iqgrid) == 0.d0) then
493        q2vq(iqgrid) = -zion / rloc * sqrt(2.d0 / pi)
494        dvloc(iqgrid) = 0.d0
495      else
496        erfValue = abi_derfc(sqrt(0.5d0 * rq2))
497        q2vq(iqgrid) = -zion / qgrid(iqgrid) * (1.0d0 - erfValue)
498        dvloc(iqgrid) = - sqrt(2.d0 / pi) * zion * gaussValue / (qgrid(iqgrid) * rloc) - &
499 &       q2vq(iqgrid) / qgrid(iqgrid)
500      end if
501 !    Add the gaussian part
502      polyValue = cc1 + cc2 * rq2 + cc3 * rq2 ** 2 + cc4 * rq2 ** 3
503      q2vq(iqgrid) = q2vq(iqgrid) + gaussValue * polyValue
504      rq = qgrid(iqgrid) / rloc
505      dvloc(iqgrid) = dvloc(iqgrid) - qgrid(iqgrid) / rloc ** 2 * gaussValue * polyValue + &
506 &     gaussValue * (2.0d0 * cc2 * rq / rloc + 3.0d0 * cc3 * rq ** 3 / rloc + &
507 &     6.0d0 * cc4 * rq ** 5 / rloc)
508    end do
509 
510    write(message, '(a,f12.7,a,a,f12.7,a,a,a,f12.7)' ) &
511 &   '  | dr spline step is : ', qgrid(2), ch10, &
512 &   '  | r > ', qgrid(mqgrid) ,' is set to 0.', ch10, &
513 &   '  | last non-nul potential value is : ', q2vq(mqgrid)
514    call wrtout(ab_out,message,'COLL')
515    call wrtout(std_out,  message,'COLL')
516  end if
517 
518 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

SOURCE

575 subroutine psp3in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion)
576 
577 !Arguments ------------------------------------
578 !scalars
579  integer,intent(in) :: ipsp,pspso
580  integer,intent(inout) :: lmax
581  real(dp),intent(in) :: zion
582  real(dp),intent(out) :: epsatm
583  type(dataset_type),intent(in) :: dtset
584  type(pseudopotential_type),intent(inout) :: psps
585 !arrays
586  integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang)
587  real(dp),intent(inout) :: ekb(psps%lnmax),ffspl(psps%mqgrid_ff,2,psps%lnmax)!vz_i
588  real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2)
589 
590 !Local variables-------------------------------
591 !scalars
592  integer :: iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nso
593  real(dp) :: cc1,cc2,cc3,cc4,h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s
594  real(dp) :: k11d,k11f,k11p,k22d,k22p,k33d,k33p,rloc
595  real(dp) :: rrd,rrf,rrp,rrs,yp1,ypn
596  character(len=500) :: message,errmsg
597 !arrays
598  real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:)
599  real(dp),allocatable :: ffspl_sr(:,:,:,:),work_space(:),work_spl(:)
600 
601 ! ***************************************************************************
602 
603 !Set various terms to 0 in case not defined below
604 !HGH values
605  rloc=zero ; rrs=zero  ; h11p=zero ; k33p=zero ; k11d=zero;
606  cc1=zero  ; h11s=zero ; h22p=zero ; rrd=zero  ; k22d=zero;
607  cc2=zero  ; h22s=zero ; h33p=zero ; h11d=zero ; k33d=zero;
608  cc3=zero  ; h33s=zero ; k11p=zero ; h22d=zero ; h11f=zero;
609  cc4=zero  ; rrp=zero  ; k22p=zero ; h33d=zero ; k11f=zero;
610  rrf=zero
611  nproj(1:psps%mpssoang)=0
612 
613 !Read and write different lines of the pseudopotential file
614 
615  read (tmp_unit,*,err=10,iomsg=errmsg) rloc,cc1,cc2,cc3,cc4
616  write(message, '(a,f12.7)' ) ' rloc=',rloc
617  call wrtout(ab_out,message,'COLL')
618  call wrtout(std_out,  message,'COLL')
619  write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' cc1 =',cc1,'; cc2 =',cc2,'; cc3 =',cc3,'; cc4 =',cc4
620  call wrtout(ab_out,message,'COLL')
621  call wrtout(std_out,  message,'COLL')
622 
623 !For the time being, the s state line must be present and is read,
624 !even for local pseudopotentials (zero must appear)
625  read (tmp_unit,*,err=10,iomsg=errmsg) rrs,h11s,h22s,h33s
626  write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrs =',rrs,'; h11s=',h11s,'; h22s=',h22s,'; h33s=',h33s
627  call wrtout(ab_out,message,'COLL')
628  call wrtout(std_out,  message,'COLL')
629 
630  if (lmax > 0) then
631 
632    read (tmp_unit,*,err=10,iomsg=errmsg) rrp,h11p,h22p,h33p
633    write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrp =',rrp,'; h11p=',h11p,'; h22p=',h22p,'; h33p=',h33p
634    call wrtout(ab_out,message,'COLL')
635    call wrtout(std_out,  message,'COLL')
636 
637    read (tmp_unit,*,err=10,iomsg=errmsg) k11p,k22p,k33p
638    write(message, '(a,f12.7,a,f12.7,a,f12.7)' )'                    k11p=',k11p,'; k22p=',k22p,'; k33p=',k33p
639    call wrtout(ab_out,message,'COLL')
640    call wrtout(std_out,  message,'COLL')
641 
642  end if
643 
644  if (lmax > 1) then
645    read (tmp_unit,*,err=10,iomsg=errmsg) rrd,h11d,h22d,h33d
646    write(message, '(a,f12.7,a,f12.7,a,f12.7,a,f12.7)' )' rrd =',rrd,'; h11d=',h11d,'; h22d=',h22d,'; h33d=',h33d
647    call wrtout(ab_out,message,'COLL')
648    call wrtout(std_out,  message,'COLL')
649 
650    read (tmp_unit,*,err=10,iomsg=errmsg) k11d,k22d,k33d
651    write(message, '(a,f12.7,a,f12.7,a,f12.7)' )'                    k11d=',k11d,'; k22d=',k22d,'; k33d=',k33d
652    call wrtout(ab_out,message,'COLL')
653    call wrtout(std_out,  message,'COLL')
654  end if
655 
656  if (lmax > 2) then
657    read (tmp_unit,*,err=10,iomsg=errmsg) rrf,h11f
658    write(message, '(a,f12.7,a,f12.7)' )' rrf =',rrf,'; h11f=',h11f
659    call wrtout(ab_out,message,'COLL')
660    call wrtout(std_out,  message,'COLL')
661 
662    read (tmp_unit,*,err=10,iomsg=errmsg) k11f
663    write(message, '(a,f12.7)' )'                    k11f=',k11f
664    call wrtout(ab_out,message,'COLL')
665    call wrtout(std_out,  message,'COLL')
666  end if
667 
668  if (abs(h11s)>1.d-08) nproj(1)=1
669  if (abs(h22s)>1.d-08) nproj(1)=2
670  if (abs(h33s)>1.d-08) nproj(1)=3
671 
672  if (abs(h11p)>1.d-08) then
673    nproj(2)=1
674    if (lmax<1) then
675      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
676 &     ' psp3in : COMMENT -',ch10,&
677 &     '  input lmax=',lmax,'  does not agree with input h11p=',h11p,ch10,&
678 &     '  setting lmax to 1'
679      call wrtout(ab_out,message,'COLL')
680      call wrtout(std_out,  message,'COLL')
681      lmax=1
682    end if
683  end if
684 
685  if (abs(h22p)>1.d-08) then
686    nproj(2)=2
687    if (lmax<1) then
688      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
689 &     ' psp3in : COMMENT -',ch10,&
690 &     '  input lmax=',lmax,' does not agree with input h22p=',h22p,ch10,&
691 &     '  setting lmax to 1'
692      call wrtout(ab_out,message,'COLL')
693      call wrtout(std_out,  message,'COLL')
694      lmax=1
695    end if
696  end if
697 
698  if (abs(h33p)>1.d-08) then
699    nproj(2)=3
700    if (lmax<1) then
701      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
702 &     ' psp3in : COMMENT -',ch10,&
703 &     '  input lmax=',lmax,' does not agree with input h33p=',h33p,ch10,&
704 &     '  setting lmax to 1'
705      call wrtout(ab_out,message,'COLL')
706      call wrtout(std_out,  message,'COLL')
707      lmax=1
708    end if
709  end if
710 
711  if (abs(h11d)>1.d-08) then
712    nproj(3)=1
713    if (lmax<2) then
714      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
715 &     ' psp3in : COMMENT -',ch10,&
716 &     '  input lmax=',lmax,'  does not agree with input h11d=',h11d,ch10,&
717 &     '  setting lmax to 2'
718      call wrtout(ab_out,message,'COLL')
719      call wrtout(std_out,  message,'COLL')
720      lmax=2
721    end if
722  end if
723 
724  if (abs(h22d)>1.d-08) then
725    nproj(3)=2
726    if (lmax<2) then
727      write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
728 &     ' psp3in : COMMENT -',ch10,&
729 &     '  input lmax=',lmax,'  does not agree with input h22d=',h22d,ch10,&
730 &     '  setting lmax to 2'
731      call wrtout(ab_out,message,'COLL')
732      call wrtout(std_out,  message,'COLL')
733      lmax=2
734    end if
735  end if
736 
737  if (abs(h33d)>1.d-08) then
738    nproj(3)=3
739    if (lmax<2) 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 h33d=',h33d,ch10,&
743 &     '  setting lmax to 2'
744      call wrtout(ab_out,message,'COLL')
745      call wrtout(std_out,  message,'COLL')
746      lmax=2
747    end if
748  end if
749 
750  if (abs(h11f)>1.d-08) then
751    nproj(4)=1
752    if (lmax<3) 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 h11f=',h11f,ch10,&
756 &     '  setting lmax to 3'
757      call wrtout(ab_out,message,'COLL')
758      call wrtout(std_out,  message,'COLL')
759      lmax=3
760    end if
761  end if
762 
763  if(pspso/=0) then
764 
765    if (abs(k11p)>1.d-08) then
766      nproj(psps%mpsang+1)=1
767      if (lmax<1) then
768        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
769 &       ' psp3in : COMMENT -',ch10,&
770 &       '  input lmax=',lmax,'  does not agree with input k11p=',k11p,ch10,&
771 &       '  setting lmax to 1'
772        call wrtout(ab_out,message,'COLL')
773        call wrtout(std_out,  message,'COLL')
774        lmax=1
775      end if
776    end if
777 
778    if (abs(k22p)>1.d-08) then
779      nproj(psps%mpsang+1)=2
780      if (lmax<1) then
781        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
782 &       ' psp3in : COMMENT -',ch10,&
783 &       '  input lmax=',lmax,' does not agree with input k22p=',k22p,ch10,&
784 &       '  setting lmax to 1'
785        call wrtout(ab_out,message,'COLL')
786        call wrtout(std_out,  message,'COLL')
787        lmax=1
788      end if
789    end if
790 
791 
792    if (abs(k33p)>1.d-08) then
793      nproj(psps%mpsang+1)=3
794      if (lmax<1) then
795        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
796 &       ' psp3in : COMMENT -',ch10,&
797 &       '  input lmax=',lmax,' does not agree with input k33p=',k33p,ch10,&
798 &       '  setting lmax to 1'
799        call wrtout(ab_out,message,'COLL')
800        call wrtout(std_out,  message,'COLL')
801        lmax=1
802      end if
803    end if
804 
805    if (abs(k11d)>1.d-08) then
806      nproj(psps%mpsang+2)=1
807      if (lmax<2) then
808        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
809 &       ' psp3in : COMMENT -',ch10,&
810 &       '  input lmax=',lmax,'  does not agree with input k11d=',k11d,ch10,&
811 &       '  setting lmax to 2'
812        call wrtout(ab_out,message,'COLL')
813        call wrtout(std_out,  message,'COLL')
814        lmax=2
815      end if
816    end if
817 
818    if (abs(k22d)>1.d-08) then
819      nproj(psps%mpsang+2)=2
820      if (lmax<2) then
821        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
822 &       ' psp3in : COMMENT -',ch10,&
823 &       '  input lmax=',lmax,'  does not agree with input k22d=',k22d,ch10,&
824 &       '  setting lmax to 2'
825        call wrtout(ab_out,message,'COLL')
826        call wrtout(std_out,  message,'COLL')
827        lmax=2
828      end if
829    end if
830 
831    if (abs(k33d)>1.d-08) then
832      nproj(psps%mpsang+2)=3
833      if (lmax<2) then
834        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
835 &       ' psp3in : COMMENT -',ch10,&
836 &       '  input lmax=',lmax,' does not agree with input k33d=',k33d,ch10,&
837 &       '  setting lmax to 2'
838        call wrtout(ab_out,message,'COLL')
839        call wrtout(std_out,  message,'COLL')
840        lmax=2
841      end if
842    end if
843 
844    if (abs(k11f)>1.d-08) then
845      nproj(psps%mpsang+3)=1
846      if (lmax<3) then
847        write(message, '(a,a,a,a,i5,a,e12.4,a,a)' ) ch10,&
848 &       ' psp3in : COMMENT -',ch10,&
849 &       '  input lmax=',lmax,' does not agree with input k11f=',k11f,ch10,&
850 &       '  setting lmax to 3'
851        call wrtout(ab_out,message,'COLL')
852        call wrtout(std_out,  message,'COLL')
853        lmax=3
854      end if
855    end if
856 
857  end if
858 
859 !Store the coefficients.
860  psps%gth_params%set(ipsp)          = .true.
861  psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc1, cc2, cc3, cc4, zero, zero /)
862  psps%gth_params%psppar(1, :, ipsp) = (/ rrs,  h11s, h22s, h33s, zero, zero, zero /)
863  psps%gth_params%psppar(2, :, ipsp) = (/ rrp,  h11p, h22p, h33p, zero, zero, zero /)
864  psps%gth_params%psppar(3, :, ipsp) = (/ rrd,  h11d, h22d, h33d, zero, zero, zero /)
865  psps%gth_params%psppar(4, :, ipsp) = (/ rrf,  h11f, zero, zero, zero, zero, zero /)
866 
867 !Store the k coefficients
868  psps%gth_params%psp_k_par(1, :, ipsp) = (/ zero, zero, zero /)
869  psps%gth_params%psp_k_par(2, :, ipsp) = (/ k11p, k22p, k33p /)
870  psps%gth_params%psp_k_par(3, :, ipsp) = (/ k11d, k22d, k33d /)
871  psps%gth_params%psp_k_par(4, :, ipsp) = (/ k11f, zero, zero /)
872 
873 !Additionnal wavelet parameters
874  if (dtset%usewvl == 1) then
875    call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
876  end if
877 
878 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn
879  nso=1;if(pspso/=0) nso=2
880  index=0;iln=0;indlmn(:,:)=0
881  do nn=1,nso
882    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
883      if (nproj(ipsang)>0) then
884        ll=ipsang-(nn-1)*lmax-1
885        do kk=1,nproj(ipsang)
886          iln=iln+1
887          do mm=1,2*ll*psps%useylm+1
888            index=index+1
889            indlmn(1,index)=ll
890            indlmn(2,index)=mm-ll*psps%useylm-1
891            indlmn(3,index)=kk
892            indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
893            indlmn(5,index)=iln
894            indlmn(6,index)=nn
895          end do
896        end do
897      end if
898    end do
899  end do
900 
901  ABI_MALLOC(dvlspl,(psps%mqgrid_ff))
902 !First, the local potential --  compute on q grid and fit spline
903  call psp2lo(cc1,cc2,cc3,cc4,dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,vlspl(:,1),rloc,.true.,yp1,ypn,zion)
904  ABI_FREE(dvlspl)
905 
906 !DEBUG
907 !write(std_out,*)' psp3in : after psp2lo '
908 !ENDDEBUG
909 
910 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
911  ABI_MALLOC(work_space,(psps%mqgrid_ff))
912  ABI_MALLOC(work_spl,(psps%mqgrid_ff))
913  call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl)
914  vlspl(:,2)=work_spl(:)
915  ABI_FREE(work_space)
916  ABI_FREE(work_spl)
917 
918 !Second, compute KB energies and form factors and fit splines
919  ekb(:)=zero
920 
921 !Check if any nonlocal projectors are being used
922  mproj=maxval(nproj)
923  if (mproj>0) then
924 
925    ABI_MALLOC(ekb_sr,(psps%mpsang,mproj))
926    ABI_MALLOC(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj))
927    ABI_MALLOC(ekb_so,(psps%mpsang,mproj))
928    ABI_MALLOC(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj))
929 
930    call psp3nl(ekb_sr,ffspl_sr,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
931 &   h33d,h11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs)
932    if(pspso/=0) then
933      call psp3nl(ekb_so,ffspl_so,zero,zero,zero,k11p,k22p,k33p,k11d,&
934 &     k22d,k33d,k11f,mproj,psps%mpsang,psps%mqgrid_ff,psps%qgrid_ff,rrd,rrf,rrp,rrs)
935    end if
936 
937 
938 !  Convert ekb and ffspl
939    iln0=0
940    do jj=1,psps%lmnmax
941      iln=indlmn(5,jj)
942      if (iln>iln0) then
943        iln0=iln
944        if (indlmn(6,jj)<=1) then
945          ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj))
946          ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj))
947        else
948          ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj))
949          ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj))
950        end if
951      end if
952    end do
953 
954    ABI_FREE(ekb_sr)
955    ABI_FREE(ffspl_sr)
956    ABI_FREE(ekb_so)
957    ABI_FREE(ffspl_so)
958  end if
959 
960  return
961 
962  ! Handle IO error
963  10 continue
964  ABI_ERROR(errmsg)
965 
966 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

SOURCE

1366 subroutine psp10in(dtset, ekb, epsatm, ffspl, indlmn, ipsp, lmax, nproj, psps, pspso, vlspl, zion)
1367 
1368 !Arguments ------------------------------------
1369 !scalars
1370  integer,intent(in) :: ipsp,pspso
1371  integer,intent(inout) :: lmax
1372  real(dp),intent(in) :: zion
1373  real(dp),intent(out) :: epsatm
1374  type(dataset_type),intent(in) :: dtset
1375  type(pseudopotential_type),intent(inout) :: psps
1376 !arrays
1377  integer,intent(out) :: indlmn(6,psps%lmnmax),nproj(psps%mpssoang)
1378  real(dp),intent(out) :: ekb(psps%lnmax) !vz_i
1379  real(dp),intent(inout) :: ffspl(psps%mqgrid_ff,2,psps%lnmax) !vz_i
1380  real(dp),intent(out) :: vlspl(psps%mqgrid_ff,2)
1381 
1382 !Local variables-------------------------------
1383 !scalars
1384  integer :: ii,iln,iln0,index,ipsang,jj,kk,ll,mm,mproj,nn,nnonloc,nprl,nso
1385  real(dp) :: rloc,yp1,ypn
1386  character(len=500) :: message,errmsg
1387 !arrays
1388  integer,allocatable :: dummy_nproj(:)
1389  real(dp) :: cc(4)
1390  real(dp),allocatable :: dvlspl(:),ekb_so(:,:),ekb_sr(:,:),ffspl_so(:,:,:,:)
1391  real(dp),allocatable :: ffspl_sr(:,:,:,:),hij(:,:,:),kij(:,:,:),rr(:)
1392  real(dp),allocatable :: work_space(:),work_spl(:)
1393 
1394 ! ***************************************************************************
1395 
1396 !Set various terms to 0 in case not defined below
1397 !HGH values
1398  rloc=zero ; cc(:)=zero
1399  nproj(1:psps%mpssoang)=0
1400 
1401 !Read and write different lines of the pseudopotential file
1402 
1403  read (tmp_unit,*,err=10,iomsg=errmsg) rloc,nn,(cc(jj),jj=1,nn)
1404  write(message, '(a,f12.7)' ) ' rloc=',rloc
1405  call wrtout(ab_out,message,'COLL')
1406  call wrtout(std_out,  message,'COLL')
1407  write(message, '(a,i1,a,4f12.7)' )' cc(1:',nn,')=',(cc(jj),jj=1,nn)
1408  call wrtout(ab_out,message,'COLL')
1409  call wrtout(std_out,  message,'COLL')
1410 
1411 !Read the number of the non-local projectors
1412  read (tmp_unit,*,err=10,iomsg=errmsg) nnonloc
1413  if (nnonloc/=lmax+1) then
1414    write(message, '(a,a,a,a,i5,a,i5,a,a,a,a,i5)' ) ch10,&
1415 &   ' psp10in : COMMENT -',ch10,&
1416 &   '  input lmax=',lmax,'  does not agree with input nnonloc=',nnonloc,ch10,&
1417 &   '  which has to be lmax+1.',ch10,&
1418 &   '  Setting lmax to ',nnonloc-1
1419    call wrtout(ab_out,message,'COLL')
1420    call wrtout(std_out,  message,'COLL')
1421    lmax=1
1422  end if
1423  ABI_MALLOC(rr,(0:lmax))
1424  ABI_MALLOC(hij,(0:lmax,3,3))
1425  ABI_MALLOC(kij,(0:lmax,3,3))
1426  rr(:)=zero; hij(:,:,:)=zero; kij(:,:,:)=zero
1427 
1428 !Read and echo the coefficients of non-local projectors
1429  prjloop: do ll=0,lmax
1430    read (tmp_unit,*,err=10,iomsg=errmsg) rr(ll),nprl,(hij(ll,1,jj),jj=1,nprl)
1431    do ii=2,nprl
1432      read (tmp_unit,*,err=10,iomsg=errmsg) (hij(ll,ii,jj),jj=ii,nprl)
1433    end do
1434    nproj(ll+1)=nprl
1435    write(message, '(a,i3,a,f12.7,2a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )&
1436 &   ' for angular momentum l =',ll,' r(l) =',rr(ll),ch10,&
1437 &   '   h11, h12, h13 =', (hij(ll,1,jj),jj=1,3),ch10,&
1438 &   '        h22, h23 =', (hij(ll,2,jj),jj=2,3),ch10,&
1439 &   '             h33 =', (hij(ll,3,jj),jj=3,3)
1440    call wrtout(ab_out,message,'COLL')
1441    call wrtout(std_out,  message,'COLL')
1442    if (ll==0) cycle
1443    do ii=1,nprl
1444      read (tmp_unit,*,err=10,iomsg=errmsg) (kij(ll,ii,jj),jj=ii,nprl)
1445    end do
1446    write(message, '(a,3f12.7,2a,12x,2f12.7,2a,24x,f12.7)' )&
1447 &   '   k11, k12, k13 =', (kij(ll,1,jj),jj=1,3),ch10,&
1448 &   '        k22, k23 =', (kij(ll,2,jj),jj=2,3),ch10,&
1449 &   '             k33 =', (kij(ll,3,jj),jj=3,3)
1450    call wrtout(ab_out,message,'COLL')
1451    call wrtout(std_out,  message,'COLL')
1452  end do prjloop
1453 
1454  if(pspso/=0) then
1455 
1456 !  MJV 10/2008: is this correct? For the normal HGH psp there are cases
1457 !  where there are more SO projectors than SR ones! e.g. Pb with 12 electrons.
1458    do ll=1,lmax
1459      nproj(psps%mpsang+ll)=nproj(ll+1)
1460    end do
1461 
1462  end if
1463 
1464 !Store the coefficients.
1465  psps%gth_params%set(ipsp)          = .true.
1466  psps%gth_params%psppar(0, :, ipsp) = (/ rloc, cc(1), cc(2), cc(3), cc(4), zero, zero /)
1467  do ii=1,4
1468    ll=ii-1
1469    if (ll>lmax) then
1470      psps%gth_params%psppar(ii,:,ipsp) = (/ zero, zero, zero, zero, zero, zero, zero /)
1471    else
1472      psps%gth_params%psppar(ii,:,ipsp) =&
1473 &     (/ 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) /)
1474    end if
1475  end do
1476 
1477 !Additionnal wavelet parameters
1478  if (dtset%usewvl == 1) then
1479    call wvl_descr_psp_fill(psps%gth_params, ipsp, 0, int(psps%zionpsp(ipsp)), int(psps%znuclpsp(ipsp)), tmp_unit)
1480  end if
1481 
1482 !Initialize array indlmn array giving l,m,n,ln,lm,s for i=lmn
1483  nso=1;if(pspso/=0) nso=2
1484  index=0;iln=0;indlmn(:,:)=0
1485  do nn=1,nso
1486    do ipsang=1+(nn-1)*(lmax+1),nn*lmax+1
1487      if (nproj(ipsang)>0) then
1488        ll=ipsang-(nn-1)*lmax-1
1489        do kk=1,nproj(ipsang)
1490          iln=iln+1
1491          do mm=1,2*ll*psps%useylm+1
1492            index=index+1
1493            indlmn(1,index)=ll
1494            indlmn(2,index)=mm-ll*psps%useylm-1
1495            indlmn(3,index)=kk
1496            indlmn(4,index)=ll*ll+(1-psps%useylm)*ll+mm
1497            indlmn(5,index)=iln
1498            indlmn(6,index)=nn
1499          end do
1500        end do
1501      end if
1502    end do
1503  end do
1504 
1505  ABI_MALLOC(dvlspl,(psps%mqgrid_ff))
1506 !First, the local potential --  compute on q grid and fit spline
1507  call psp2lo(cc(1),cc(2),cc(3),cc(4),dvlspl,epsatm,psps%mqgrid_ff,psps%qgrid_ff,&
1508 & vlspl(:,1),rloc,.true.,yp1,ypn,zion)
1509  ABI_FREE(dvlspl)
1510 
1511 !Fit spline to q^2 V(q) (Numerical Recipes subroutine)
1512  ABI_MALLOC(work_space,(psps%mqgrid_ff))
1513  ABI_MALLOC(work_spl,(psps%mqgrid_ff))
1514  call spline (psps%qgrid_ff,vlspl(:,1),psps%mqgrid_ff,yp1,ypn,work_spl)
1515  vlspl(:,2)=work_spl(:)
1516  ABI_FREE(work_space)
1517  ABI_FREE(work_spl)
1518 
1519 !Second, compute KB energies and form factors and fit splines
1520  ekb(:)=zero
1521 
1522 !Check if any nonlocal projectors are being used
1523  mproj=maxval(nproj)
1524 
1525  if (mproj>0) then
1526 
1527    ABI_MALLOC(ekb_sr,(psps%mpsang,mproj))
1528    ABI_MALLOC(ffspl_sr,(psps%mqgrid_ff,2,psps%mpsang,mproj))
1529    ABI_MALLOC(ekb_so,(psps%mpsang,mproj))
1530    ABI_MALLOC(ffspl_so,(psps%mqgrid_ff,2,psps%mpsang,mproj))
1531 
1532    call psp10nl(ekb_sr,ffspl_sr,hij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,&
1533 &   nproj,psps%qgrid_ff,rr)
1534    if(pspso/=0) then
1535      ABI_MALLOC(dummy_nproj,(psps%mpsang))
1536      dummy_nproj(1)=0
1537      do ll=1,lmax
1538        dummy_nproj(ll+1)=nproj(psps%mpsang+ll)
1539      end do
1540      call psp10nl(ekb_so,ffspl_so,kij,lmax,mproj,psps%mpsang,psps%mqgrid_ff,&
1541 &     dummy_nproj,psps%qgrid_ff,rr)
1542      ABI_FREE(dummy_nproj)
1543    end if
1544 
1545 !  Convert ekb and ffspl
1546    iln0=0
1547    do jj=1,psps%lmnmax
1548      iln=indlmn(5,jj)
1549      if (iln>iln0) then
1550        iln0=iln
1551        if (indlmn(6,jj)<=1) then
1552          ekb(iln)=ekb_sr(1+indlmn(1,jj),indlmn(3,jj))
1553          ffspl(:,:,iln)=ffspl_sr(:,:,1+indlmn(1,jj),indlmn(3,jj))
1554        else
1555          ekb(iln)=ekb_so(1+indlmn(1,jj),indlmn(3,jj))
1556          ffspl(:,:,iln)=ffspl_so(:,:,1+indlmn(1,jj),indlmn(3,jj))
1557        end if
1558      end if
1559    end do
1560 
1561    ABI_FREE(ekb_sr)
1562    ABI_FREE(ffspl_sr)
1563    ABI_FREE(ekb_so)
1564    ABI_FREE(ffspl_so)
1565  end if
1566 
1567  ABI_FREE(rr)
1568  ABI_FREE(hij)
1569  ABI_FREE(kij)
1570 
1571  return
1572 
1573  ! Handle IO error
1574  10 continue
1575  ABI_ERROR(errmsg)
1576 
1577 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

SOURCE

1607 subroutine psp10nl(ekb,ffspl,hij,lmax,mproj,mpsang,mqgrid,nproj,qgrid,rr)
1608 
1609 !Arguments ------------------------------------
1610 !scalars
1611  integer,intent(in) :: lmax,mproj,mpsang,mqgrid
1612 !arrays
1613  integer,intent(in) :: nproj(mpsang)
1614  real(dp),intent(in) :: hij(0:lmax,3,3),qgrid(mqgrid),rr(0:lmax)
1615  real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
1616 
1617 !Local variables-------------------------------
1618 !scalars
1619  integer :: info,ipack,iproj,iqgrid,jproj,ll,numproj
1620  real(dp) :: qmax,rrl
1621  character(len=500) :: message
1622  character :: jobz,uplo
1623 !arrays
1624  real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
1625  real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
1626 
1627 ! *************************************************************************
1628 
1629  ABI_MALLOC(ppspl,(mqgrid,2,mpsang,mproj))
1630  ABI_MALLOC(work,(mqgrid))
1631 
1632  qmax=qgrid(mqgrid)
1633  jobz='v'
1634  uplo='u'
1635  ekb(:,:)=zero
1636 
1637  lloop: do ll=0,lmax
1638    ap(:,:)=zero
1639    numproj=nproj(ll+1)
1640 
1641 !  Fill up the matrix in packed storage
1642    prjloop: do jproj=1,numproj
1643      priloop: do iproj=1,jproj
1644        ipack=iproj+(jproj-1)*jproj/2
1645        if(mod((jproj-1)*jproj,2)/=0) then
1646          ABI_ERROR("odd")
1647        end if
1648        ap(1,ipack)=hij(ll,iproj,jproj)
1649      end do priloop
1650    end do prjloop
1651 
1652    if(numproj/=0)then
1653 
1654      ABI_MALLOC(uu,(numproj,numproj))
1655      ABI_MALLOC(zz,(2,numproj,numproj))
1656 
1657      if (numproj > 1) then
1658        call ZHPEV(jobz,uplo,numproj,ap,ww,zz,numproj,work1,rwork1,info)
1659        uu(:,:)=zz(1,:,:)
1660      else
1661        ww(1)=hij(ll,1,1)
1662        uu(1,1)=one
1663      end if
1664 
1665 !    Initialization of ekb, and spline fitting
1666 
1667      if (ll==0) then ! s channel
1668 
1669        rrl=rr(0)
1670        do iproj=1,numproj
1671          ekb(1,iproj)=ww(iproj)*32.d0*(rrl**3)*(pi**2.5d0)/(4.d0*pi)**2
1672          if(iproj==1)then
1673            do iqgrid=1,mqgrid
1674              ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2)
1675            end do
1676            yp1j(1)=zero
1677            ypnj(1)=-(two_pi*rrl)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrl)**2)
1678          else if(iproj==2)then
1679            do iqgrid=1,mqgrid
1680              ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0)     &
1681 &             *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) &
1682 &             *( 3.d0-(two_pi*qgrid(iqgrid)*rrl)**2 )
1683            end do
1684            yp1j(2)=zero
1685            ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrl)**2*qmax &
1686 &           *exp(-0.5d0*(two_pi*qmax*rrl)**2) * (-5.d0+(two_pi*qmax*rrl)**2)
1687          else if(iproj==3)then
1688            do iqgrid=1,mqgrid
1689              ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
1690 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1691 &             (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrl)**2 + &
1692 &             (two_pi*qgrid(iqgrid)*rrl)**4)
1693            end do
1694            yp1j(3)=zero
1695            ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2) * &
1696 &           (two_pi*rrl)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrl)**2-(two_pi*qmax*rrl)**4)
1697          end if
1698          call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
1699 &         yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
1700        end do
1701 
1702      else if (ll==1) then ! p channel
1703 
1704        rrl=rr(1)
1705        do iproj=1,numproj
1706          ekb(2,iproj)=ww(iproj)*64.d0*(rrl**5)*(pi**2.5d0)/(4.d0*pi)**2
1707          if(iproj==1)then
1708            do iqgrid=1,mqgrid
1709              ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
1710 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))
1711            end do
1712            yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
1713            ypnj(1)=-two_pi*((two_pi*qmax*rrl)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)*&
1714 &           (1.0d0/sqrt(3.0d0))
1715          else if(iproj==2)then
1716            do iqgrid=1,mqgrid
1717              ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
1718 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1719 &             (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrl)**2)
1720            end do
1721            yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
1722            ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1723 &           (-8*(two_pi*qmax*rrl)**2 + (two_pi*qmax*rrl)**4 + 5.0d0)
1724          else if(iproj==3)then
1725            do iqgrid=1,mqgrid
1726              ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
1727 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1728 &             (two_pi*qgrid(iqgrid))*&
1729 &             (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrl)**2+(two_pi*qgrid(iqgrid)*rrl)**4)
1730            end do
1731            yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
1732            ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1733 &           (35.0d0-77.0d0*(two_pi*qmax*rrl)**2+19.0d0*(two_pi*qmax*rrl)**4 - &
1734 &           (two_pi*qmax*rrl)**6)
1735          end if
1736          call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
1737 &         yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
1738        end do
1739 
1740      else if (ll==2) then ! d channel
1741 
1742 !      If there is a third projector. Warning : only two projectors are allowed.
1743        if ( numproj>2 ) then
1744          write(message, '(3a)' )&
1745 &         ' only two d-projectors are allowed ',ch10,&
1746 &         ' Action: check your pseudopotential file.'
1747          ABI_ERROR(message)
1748        end if
1749 
1750        rrl=rr(2)
1751        do iproj=1,numproj
1752          ekb(3,iproj)=ww(iproj)*128.d0*(rrl**7)*(pi**2.5d0)/(4.d0*pi)**2
1753          if(iproj==1)then
1754            do iqgrid=1,mqgrid
1755              ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
1756 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * (two_pi*qgrid(iqgrid))**2
1757            end do
1758            yp1j(1)=zero
1759            ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
1760 &           exp(-0.5d0*(two_pi*qmax*rrl)**2)*qmax*(2d0-(two_pi*qmax*rrl)**2)
1761          else if(iproj==2)then
1762            do iqgrid=1,mqgrid
1763              ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
1764 &             exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2) * &
1765 &             ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrl)**2)
1766            end do
1767            yp1j(2)=zero
1768            ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrl)**2)* &
1769 &           qmax*(two_pi**2)*( (two_pi*qmax*rrl)**4 - 11.0d0*(two_pi*qmax*rrl)**2 + 14.0d0)
1770          end if
1771          call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
1772 &         yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
1773        end do
1774 
1775      else if (ll==3) then ! f channel
1776 
1777 !      If there is a second projector. Warning : only one projector is allowed.
1778        if ( numproj>1 ) then
1779          write(message, '(a,a,a)' )&
1780 &         'only one f-projector is allowed ',ch10,&
1781 &         'Action: check your pseudopotential file.'
1782          ABI_ERROR(message)
1783        end if
1784 
1785        rrl=rr(3)
1786        ekb(4,1)=ww(1)*(256.0d0/105.0d0)*(rrl**9)*(pi**2.5d0)/(4.d0*pi)**2
1787        do iqgrid=1,mqgrid
1788          ppspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
1789 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrl)**2)
1790        end do
1791 !      Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
1792        yp1j(1)=zero
1793        ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrl)**2)*&
1794 &       (3.0d0-(two_pi*qmax*rrl)**2)
1795 !      Fit spline to get second derivatives by spline fit
1796        call spline(qgrid,ppspl(:,1,4,1),mqgrid,&
1797 &       yp1j(1),ypnj(1),ppspl(:,2,4,1))
1798 
1799      else
1800        ABI_ERROR("lmax>3?")
1801      end if
1802 
1803 !    Linear combination using the eigenvectors
1804      ffspl(:,:,ll+1,:)=zero
1805      do jproj=1,numproj
1806        do iproj=1,numproj
1807          do iqgrid=1,mqgrid
1808            ffspl(iqgrid,1:2,ll+1,jproj)=ffspl(iqgrid,1:2,ll+1,jproj) &
1809 &           +uu(iproj,jproj)*ppspl(iqgrid,1:2,ll+1,iproj)
1810          end do
1811        end do
1812      end do
1813 
1814      ABI_FREE(uu)
1815      ABI_FREE(zz)
1816 
1817 !    End condition on numproj(/=0)
1818    end if
1819 
1820  end do lloop
1821 
1822  ABI_FREE(ppspl)
1823  ABI_FREE(work)
1824 
1825 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.

SOURCE

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

SOURCE

282 subroutine psp2nl(ekb,ffspl,h1p,h1s,h2s,lnmax,mqgrid,qgrid,rrp,rrs)
283 
284 !Arguments ------------------------------------
285 !scalars
286  integer,intent(in) :: lnmax,mqgrid
287  real(dp),intent(in) :: h1p,h1s,h2s,rrp,rrs
288 !arrays
289  real(dp),intent(in) :: qgrid(mqgrid)
290  real(dp),intent(inout) :: ekb(lnmax),ffspl(mqgrid,2,lnmax) !vz_i
291 
292 !Local variables-------------------------------
293 !scalars
294  integer :: iln,iqgrid
295  real(dp) :: qmax,yp1,ypn
296 !arrays
297  real(dp),allocatable :: work(:)
298 
299 ! *************************************************************************
300 
301  ABI_MALLOC(work,(mqgrid))
302 
303 !Kleinman-Bylander energies ekb were set to zero in calling program
304 
305 !Compute KB energies
306  iln=0
307  if (abs(h1s)>1.d-12) then
308    iln=iln+1
309    ekb(iln)=h1s*32.d0*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2)
310  end if
311  if (abs(h2s)>1.d-12) then
312    iln=iln+1
313    ekb(iln) =h2s*(128.d0/15.d0)*rrs**3*(pi**(2.5d0)/(4.d0*pi)**2)
314  end if
315  if (abs(h1p)>1.d-12) then
316    iln=iln+1
317    ekb(iln)=h1p*(64.d0/3.d0)*rrp**5*(pi**(2.5d0)/(4.d0*pi)**2)
318  end if
319 
320 !Compute KB form factor
321  iln=0
322 
323 !l=0 first projector
324  if (abs(h1s)>1.d-12) then
325    iln=iln+1
326    do iqgrid=1,mqgrid
327      ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
328    end do
329 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
330    yp1=0.d0
331    qmax=qgrid(mqgrid)
332    ypn=-4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2)
333 !  Fit spline to get second derivatives by spline fit
334    call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
335 !  else
336 !  or else put first projector nonlocal correction at l=0 to 0
337 !  ffspl(:,:,iln)=0.0d0
338  end if
339 
340 !l=0 second projector
341  if (abs(h2s)>1.d-12) then
342    iln=iln+1
343    do iqgrid=1,mqgrid
344      ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
345 &     (3.d0-(two_pi*qgrid(iqgrid)*rrs)**2)
346    end do
347 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
348    yp1=0.d0
349    qmax=qgrid(mqgrid)
350    ypn=4.d0*pi**2*qmax*rrs**2*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
351 &   (-5.d0+(two_pi*qmax*rrs)**2)
352 !  Fit spline to get second derivatives by spline fit
353    call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
354 !  else if(mproj>=2)then
355 !  or else put second projector nonlocal correction at l=0 to 0
356 !  ffspl(:,:,iln)=0.0d0
357  end if
358 
359 !l=1 first projector
360  if (abs(h1p)>1.d-12) then
361    iln=iln+1
362    do iqgrid=1,mqgrid
363      ffspl(iqgrid,1,iln)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
364 &     (two_pi*qgrid(iqgrid))
365    end do
366 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
367    yp1=two_pi
368    qmax=qgrid(mqgrid)
369    ypn=-two_pi*((two_pi*qmax*rrp)**2-1.d0) * exp(-0.5d0*(two_pi*qmax*rrp)**2)
370 !  Fit spline to get second derivatives by spline fit
371    call spline(qgrid,ffspl(:,1,iln),mqgrid,yp1,ypn,ffspl(:,2,iln))
372 !  else if(mpsang>=2)then
373 !  or else put first projector l=1 nonlocal correction to 0
374 !  ffspl(:,:,iln)=0.0d0
375  end if
376 
377  ABI_FREE(work)
378 
379 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

SOURCE

1005 subroutine psp3nl(ekb,ffspl,h11s,h22s,h33s,h11p,h22p,h33p,h11d,h22d,&
1006 &                  h33d,h11f,mproj,mpsang,mqgrid,qgrid,rrd,rrf,rrp,rrs)
1007 
1008 !Arguments ------------------------------------
1009 !scalars
1010  integer,intent(in) :: mproj,mpsang,mqgrid
1011  real(dp),intent(in) :: h11d,h11f,h11p,h11s,h22d,h22p,h22s,h33d,h33p,h33s,rrd
1012  real(dp),intent(in) :: rrf,rrp,rrs
1013 !arrays
1014  real(dp),intent(in) :: qgrid(mqgrid)
1015  real(dp),intent(out) :: ekb(mpsang,mproj),ffspl(mqgrid,2,mpsang,mproj)
1016 
1017 !Local variables-------------------------------
1018 !scalars
1019  integer :: info,iproj,iqgrid,ldz,mu,nproj,nu
1020  real(dp) :: qmax
1021  character(len=500) :: message
1022  character :: jobz,uplo
1023 !arrays
1024  real(dp) :: ap(2,9),rwork1(9),work1(2,9),ww(3),yp1j(3),ypnj(3)
1025  real(dp),allocatable :: ppspl(:,:,:,:),uu(:,:),work(:),zz(:,:,:)
1026 
1027 ! *************************************************************************
1028 
1029  ABI_MALLOC(ppspl,(mqgrid,2,mpsang,mproj))
1030  ABI_MALLOC(work,(mqgrid))
1031 
1032  qmax=qgrid(mqgrid)
1033  jobz='v'
1034  uplo='u'
1035  ekb(:,:)=0.0d0
1036 
1037 !---------------------------------------------------------------
1038 !Treat s channel
1039 
1040  nproj=0
1041  ap(:,:)=0.0d0
1042 !If there is at least one s-projector
1043  if  ( abs(h11s) >= 1.0d-8 ) then
1044    nproj=1 ; ldz=1 ; ap(1,1)=h11s
1045  end if
1046  nproj=1
1047 !If there is a second projector
1048  if  ( abs(h22s) >= 1.0d-8 ) then
1049    nproj=2 ; ldz=2 ; ap(1,3)=h22s
1050    ap(1,2)=-0.5d0*sqrt(0.6d0)*h22s
1051  end if
1052 !If there is a third projector
1053  if ( abs(h33s) >= 1.0d-8 ) then
1054    nproj=3 ; ldz=3 ; ap(1,6)=h33s
1055    ap(1,4)=0.5d0*sqrt(5.d0/21.d0)*h33s
1056    ap(1,5)=-0.5d0*sqrt(100.d0/63.d0)*h33s
1057  end if
1058 
1059  if(nproj/=0)then
1060 
1061    ABI_MALLOC(uu,(nproj,nproj))
1062    ABI_MALLOC(zz,(2,nproj,nproj))
1063 
1064    if (nproj > 1) then
1065      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1066      uu(:,:)=zz(1,:,:)
1067    else
1068      ww(1)=h11s
1069      uu(1,1)=1.0d0
1070    end if
1071 
1072 !  Initialization of ekb, and spline fitting
1073    do iproj=1,nproj
1074      ekb(1,iproj)=ww(iproj)*32.d0*(rrs**3)*(pi**2.5d0)/(4.d0*pi)**2
1075      if(iproj==1)then
1076        do iqgrid=1,mqgrid
1077          ppspl(iqgrid,1,1,1)=exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2)
1078        end do
1079        yp1j(1)=0.d0
1080        ypnj(1)=-(two_pi*rrs)**2*qmax*exp(-0.5d0*(two_pi*qmax*rrs)**2)
1081      else if(iproj==2)then
1082        do iqgrid=1,mqgrid
1083          ppspl(iqgrid,1,1,2)=2.0d0/sqrt(15.0d0)     &
1084 &         *exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) &
1085 &         *( 3.d0-(two_pi*qgrid(iqgrid)*rrs)**2 )
1086        end do
1087        yp1j(2)=0.0d0
1088        ypnj(2)=2.0d0/sqrt(15.0d0)*(two_pi*rrs)**2*qmax &
1089 &       *exp(-0.5d0*(two_pi*qmax*rrs)**2) * (-5.d0+(two_pi*qmax*rrs)**2)
1090      else if(iproj==3)then
1091        do iqgrid=1,mqgrid
1092          ppspl(iqgrid,1,1,3)=(4.0d0/3.0d0)/sqrt(105.0d0)*&
1093 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrs)**2) * &
1094 &         (15.0d0-10.0d0*(two_pi*qgrid(iqgrid)*rrs)**2 + &
1095 &         (two_pi*qgrid(iqgrid)*rrs)**4)
1096        end do
1097        yp1j(3)=0.0d0
1098        ypnj(3)=(4.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrs)**2) * &
1099 &       (two_pi*rrs)**2*qmax*(-35.0d0+14d0*(two_pi*qmax*rrs)**2-(two_pi*qmax*rrs)**4)
1100      end if
1101      call spline(qgrid,ppspl(:,1,1,iproj),mqgrid,&
1102 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,1,iproj))
1103    end do
1104 
1105 !  Linear combination using the eigenvectors
1106    ffspl(:,:,1,:)=0.0d0
1107    do mu=1,nproj
1108      do nu=1,nproj
1109        do iqgrid=1,mqgrid
1110          ffspl(iqgrid,1:2,1,mu)=ffspl(iqgrid,1:2,1,mu) &
1111 &         +uu(nu,mu)*ppspl(iqgrid,1:2,1,nu)
1112        end do
1113      end do
1114    end do
1115 
1116    ABI_FREE(uu)
1117    ABI_FREE(zz)
1118  end if !  End condition on nproj(/=0)
1119 
1120 !--------------------------------------------------------------------
1121 !Now treat p channel
1122 
1123  nproj=0
1124  ap(:,:)=0.0d0
1125 !If there is at least one projector
1126  if  ( abs(h11p) >= 1.0d-8 ) then
1127    nproj=1 ; ldz=1 ; ap(1,1)=h11p
1128  end if
1129 !If there is a second projector
1130  if  ( abs(h22p) >= 1.0d-8 ) then
1131    nproj=2 ; ldz=2 ; ap(1,3)=h22p
1132    ap(1,2)=-0.5d0*sqrt(5.d0/7.d0)*h22p
1133  end if
1134 !If there is a third projector
1135  if ( abs(h33p) >= 1.0d-8 ) then
1136    nproj=3 ; ldz=3 ; ap(1,6)=h33p
1137    ap(1,4)= (1.d0/6.d0)*sqrt(35.d0/11.d0)*h33p
1138    ap(1,5)=-(1.d0/6.d0)*(14.d0/sqrt(11.d0))*h33p
1139  end if
1140 
1141  if(nproj/=0)then
1142 
1143    ABI_MALLOC(uu,(nproj,nproj))
1144    ABI_MALLOC(zz,(2,nproj,nproj))
1145 
1146    if (nproj > 1) then
1147      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1148      uu(:,:)=zz(1,:,:)
1149    else
1150      ww(1)=h11p
1151      uu(1,1)=1.0d0
1152    end if
1153 
1154 !  Initialization of ekb, and spline fitting
1155    do iproj=1,nproj
1156      ekb(2,iproj)=ww(iproj)*64.d0*(rrp**5)*(pi**2.5d0)/(4.d0*pi)**2
1157      if(iproj==1)then
1158        do iqgrid=1,mqgrid
1159          ppspl(iqgrid,1,2,1)=(1.0d0/sqrt(3.0d0))* &
1160 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * (two_pi*qgrid(iqgrid))
1161        end do
1162        yp1j(1)=two_pi*(1.0d0/sqrt(3.0d0))
1163        ypnj(1)=-two_pi*((two_pi*qmax*rrp)**2-1.d0)*exp(-0.5d0*(two_pi*qmax*rrp)**2)*&
1164 &       (1.0d0/sqrt(3.0d0))
1165      else if(iproj==2)then
1166        do iqgrid=1,mqgrid
1167          ppspl(iqgrid,1,2,2)=(2.0d0/sqrt(105.0d0))* &
1168 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
1169 &         (two_pi*qgrid(iqgrid))*(5.0d0-(two_pi*qgrid(iqgrid)*rrp)**2)
1170        end do
1171        yp1j(2)=(5.0d0*two_pi)*(2.0d0/sqrt(105.0d0))
1172        ypnj(2)=(2.0d0/sqrt(105.0d0))*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
1173 &       (-8*(two_pi*qmax*rrp)**2 + (two_pi*qmax*rrp)**4 + 5.0d0)
1174      else if(iproj==3)then
1175        do iqgrid=1,mqgrid
1176          ppspl(iqgrid,1,2,3)=(4.0d0/3.0d0)/sqrt(1155d0)*&
1177 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrp)**2) * &
1178 &         (two_pi*qgrid(iqgrid))*&
1179 &         (35.0d0-14.0d0*(two_pi*qgrid(iqgrid)*rrp)**2+(two_pi*qgrid(iqgrid)*rrp)**4)
1180        end do
1181        yp1j(3)=(35.0d0*two_pi)*(4.0d0/3.0d0)/sqrt(1155.0d0)
1182        ypnj(3)=(4.0d0/3.0d0)/sqrt(1155.0d0)*two_pi*exp(-0.5d0*(two_pi*qmax*rrp)**2)* &
1183 &       (35.0d0-77.0d0*(two_pi*qmax*rrp)**2+19.0d0*(two_pi*qmax*rrp)**4 - &
1184 &       (two_pi*qmax*rrp)**6)
1185      end if
1186      call spline(qgrid,ppspl(:,1,2,iproj),mqgrid,&
1187 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,2,iproj))
1188    end do
1189 
1190 !  Linear combination using the eigenvectors
1191    ffspl(:,:,2,:)=0.0d0
1192    do mu=1,nproj
1193      do nu=1,nproj
1194        do iqgrid=1,mqgrid
1195          ffspl(iqgrid,1:2,2,mu)=ffspl(iqgrid,1:2,2,mu) &
1196 &         +uu(nu,mu)*ppspl(iqgrid,1:2,2,nu)
1197        end do
1198      end do
1199    end do
1200 
1201    ABI_FREE(uu)
1202    ABI_FREE(zz)
1203  end if !  End condition on nproj(/=0)
1204 
1205 !-----------------------------------------------------------------------
1206 !Now treat d channel.
1207 
1208  nproj=0
1209  ap(:,:)=0.0d0
1210 !If there is at least one projector
1211  if  ( abs(h11d) >= 1.0d-8 ) then
1212    nproj=1 ; ldz=1 ; ap(1,1)=h11d
1213  end if
1214 !If there is a second projector
1215  if  ( abs(h22d) >= 1.0d-8 ) then
1216    nproj=2 ; ldz=2 ; ap(1,3)=h22d
1217    ap(1,2)=-0.5d0*sqrt(7.d0/9.d0)*h22d
1218  end if
1219 !If there is a third projector. Warning : only two projectors are allowed.
1220  if ( abs(h33d) >= 1.0d-8 ) then
1221    write(message, '(a,a,a)' )&
1222 &   '  only two d-projectors are allowed ',ch10,&
1223 &   '  Action: check your pseudopotential file.'
1224    ABI_ERROR(message)
1225 !  nproj=3 ; ldz=3 ; ap(1,6)=h33d
1226 !  ap(1,4)= 0.5d0*sqrt(63.d0/143.d0)*h33d
1227 !  ap(1,5)= -0.5d0*(18.d0/sqrt(143.d0))*h33d
1228  end if
1229 
1230  if(nproj/=0)then
1231 
1232    ABI_MALLOC(uu,(nproj,nproj))
1233    ABI_MALLOC(zz,(2,nproj,nproj))
1234 
1235    if (nproj > 1) then
1236      call ZHPEV(jobz,uplo,nproj,ap,ww,zz,ldz,work1,rwork1,info)
1237      uu(:,:)=zz(1,:,:)
1238    else
1239      ww(1)=h11d
1240      uu(1,1)=1.0d0
1241    end if
1242 
1243 !  Initialization of ekb, and spline fitting
1244    do iproj=1,nproj
1245      ekb(3,iproj)=ww(iproj)*128.d0*(rrd**7)*(pi**2.5d0)/(4.d0*pi)**2
1246      if(iproj==1)then
1247        do iqgrid=1,mqgrid
1248          ppspl(iqgrid,1,3,1)=(1.0d0/sqrt(15.0d0))* &
1249 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * (two_pi*qgrid(iqgrid))**2
1250        end do
1251        yp1j(1)=0.0d0
1252        ypnj(1)=(1.0d0/sqrt(15.0d0))*(two_pi**2)*&
1253 &       exp(-0.5d0*(two_pi*qmax*rrd)**2)*qmax*(2d0-(two_pi*qmax*rrd)**2)
1254      else if(iproj==2)then
1255        do iqgrid=1,mqgrid
1256          ppspl(iqgrid,1,3,2)=(2.0d0/3.0d0)/sqrt(105.0d0)* &
1257 &         exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrd)**2) * &
1258 &         ((two_pi*qgrid(iqgrid))**2)*(7.0d0-(two_pi*qgrid(iqgrid)*rrd)**2)
1259        end do
1260        yp1j(2)=0.0d0
1261        ypnj(2)=(2.0d0/3.0d0)/sqrt(105.0d0)*exp(-0.5d0*(two_pi*qmax*rrd)**2)* &
1262 &       qmax*(two_pi**2)*( (two_pi*qmax*rrd)**4 - 11.0d0*(two_pi*qmax*rrd)**2 + 14.0d0)
1263      end if
1264      call spline(qgrid,ppspl(:,1,3,iproj),mqgrid,&
1265 &     yp1j(iproj),ypnj(iproj),ppspl(:,2,3,iproj))
1266    end do
1267 
1268 !  Linear combination using the eigenvectors
1269    ffspl(:,:,3,:)=0.0d0
1270    do mu=1,nproj
1271      do nu=1,nproj
1272        do iqgrid=1,mqgrid
1273          ffspl(iqgrid,1:2,3,mu)=ffspl(iqgrid,1:2,3,mu) &
1274 &         +uu(nu,mu)*ppspl(iqgrid,1:2,3,nu)
1275        end do
1276      end do
1277    end do
1278 
1279    ABI_FREE(uu)
1280    ABI_FREE(zz)
1281  end if !  End condition on nproj(/=0)
1282 
1283 !-----------------------------------------------------------------------
1284 !Treat now f channel (max one projector ! - so do not use ppspl)
1285 
1286 !l=3 first projector
1287  if (abs(h11f)>1.d-12) then
1288    ekb(4,1)=h11f*(256.0d0/105.0d0)*(rrf**9)*(pi**2.5d0)/(4.d0*pi)**2
1289    do iqgrid=1,mqgrid
1290      ffspl(iqgrid,1,4,1)=(two_pi*qgrid(iqgrid))**3* &
1291 &     exp(-0.5d0*(two_pi*qgrid(iqgrid)*rrf)**2)
1292    end do
1293 !  Compute yp1,ypn=derivatives of f(q) at q=0, q=qgrid(mqgrid)
1294    yp1j(1)=0d0
1295    ypnj(1)=(two_pi**3)*qmax**2*exp(-0.5d0*(two_pi*qmax*rrf)**2)*&
1296 &   (3.0d0-(two_pi*qmax*rrf)**2)
1297 !  Fit spline to get second derivatives by spline fit
1298    call spline(qgrid,ffspl(:,1,4,1),mqgrid,&
1299 &   yp1j(1),ypnj(1),ffspl(:,2,4,1))
1300  end if
1301 
1302 !-----------------------------------------------------------------------
1303 
1304  ABI_FREE(ppspl)
1305  ABI_FREE(work)
1306 
1307 end subroutine psp3nl