TABLE OF CONTENTS


ABINIT/m_pawpwij [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pawpwij

FUNCTION

  This module defines methods to calculate the onsite contribution of a plane wave in the PAW method:

    - pawpwff_t: Form factors used to calculate the onsite contributions of a plane wave.
    - pawpwij_t: Onsite matrix elements of a plane wave for a given atom type.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG,GKA)
 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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_pawpwij
27 
28  use defs_basis
29  use defs_datatypes
30  use defs_abitypes
31  use m_errors
32  use m_abicore
33  use m_fft
34 
35  use m_numeric_tools,  only : arth
36  use m_geometry,       only : metric
37  use m_paw_numeric,    only : paw_jbessel_4spline, paw_spline
38  use m_splines,        only : splfit
39  use m_pawang,         only : pawang_type, realgaunt
40  use m_pawrad,         only : pawrad_type, pawrad_init, pawrad_free, pawrad_copy, simp_gen
41  use m_pawtab,         only : pawtab_type
42  use m_pawcprj,        only : pawcprj_type
43  use m_mpinfo,         only : destroy_mpi_enreg, initmpi_seq
44  use m_initylmg,       only : initylmg
45 
46  implicit none
47 
48  private

ABINIT/paw_cross_rho_tw_g [ Functions ]

[ Top ] [ Functions ]

NAME

  paw_cross_rho_tw_g

FUNCTION

  Compute the cross term between the PAW onsite part and plane-wave part
  in rho_tw

INPUTS

OUTPUT

PARENTS

      calc_sigc_me,calc_sigx_me,cchi0,cchi0q0,prep_calc_ucrpa

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

1286 subroutine paw_cross_rho_tw_g(nspinor,npwvec,nr,ngfft,map2sphere,use_padfft,igfftg0,gbound,&
1287 & ur_ae1,ur_ae_onsite1,ur_ps_onsite1,i1,ktabr1,ktabp1,spinrot1,&
1288 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,i2,ktabr2,ktabp2,spinrot2,&
1289 & dim_rtwg,rhotwg)
1290 
1291 
1292 !This section has been created automatically by the script Abilint (TD).
1293 !Do not modify the following lines by hand.
1294 #undef ABI_FUNC
1295 #define ABI_FUNC 'paw_cross_rho_tw_g'
1296 !End of the abilint section
1297 
1298  implicit none
1299 
1300 !Arguments ------------------------------------
1301 !scalars
1302  integer,intent(in) :: i1,i2,npwvec,nr,nspinor,dim_rtwg,map2sphere,use_padfft
1303  complex(dpc),intent(in) :: ktabp1,ktabp2
1304 !arrays
1305  integer,intent(in) :: gbound(:,:)
1306  integer,intent(in) :: igfftg0(npwvec*map2sphere),ngfft(18)
1307  integer,intent(in) :: ktabr1(nr),ktabr2(nr)
1308  real(dp),intent(in) :: spinrot1(4),spinrot2(4)
1309  complex(gwpc),intent(in) :: ur_ae1(nr),ur_ae2(nr)
1310  complex(gwpc),intent(in) :: ur_ae_onsite1(nr),ur_ae_onsite2(nr)
1311  complex(gwpc),intent(in) :: ur_ps_onsite1(nr),ur_ps_onsite2(nr)
1312  complex(gwpc),intent(inout) :: rhotwg(npwvec*dim_rtwg)
1313 
1314 !Local variables-------------------------------
1315 !scalars
1316  integer,parameter :: ndat1=1
1317  integer :: ig,igfft,nx,ny,nz,ldx,ldy,ldz,mgfft,isprot1,isprot2
1318  type(fftbox_plan3_t) :: plan
1319 !arrays
1320  complex(dpc),allocatable :: usk(:),uu(:),rho(:)
1321 
1322 ! *************************************************************************
1323 
1324  SELECT CASE (nspinor)
1325 
1326  CASE (1) ! Collinear case.
1327 
1328    ABI_MALLOC(uu,(nr))
1329    ABI_MALLOC(usk,(nr))
1330    ABI_MALLOC(rho,(nr))
1331 
1332    uu  = ur_ae1(ktabr1)*ktabp1        - ur_ae_onsite1(ktabr1)*ktabp1; if (i1==1) uu  = CONJG(uu)
1333    usk = ur_ae_onsite2(ktabr2)*ktabp2 - ur_ps_onsite2(ktabr2)*ktabp2; if (i2==2) usk = CONJG(usk)
1334    rho = uu * usk
1335 
1336    uu  = ur_ae_onsite1(ktabr1)*ktabp1 - ur_ps_onsite1(ktabr1)*ktabp1; if (i1==1) uu  = CONJG(uu)
1337    usk = ur_ae2(ktabr2)*ktabp2        - ur_ae_onsite2(ktabr2)*ktabp2; if (i2==2) usk = CONJG(usk)
1338    rho = rho + uu * usk
1339 
1340    SELECT CASE (map2sphere)
1341 
1342    CASE (0) ! Need results on the full FFT box thus cannot use zero-padded FFT.
1343 
1344      call fftbox_plan3_many(plan,ndat1,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
1345      call fftbox_execute(plan,rho)
1346 
1347      rhotwg=rhotwg + rho
1348 
1349    CASE (1) ! Need results on the G-sphere. Call zero-padded FFT routines if required.
1350 
1351      if (use_padfft==1) then
1352        nx =ngfft(1); ny =ngfft(2); nz =ngfft(3); mgfft = MAXVAL(ngfft(1:3))
1353        ldx=nx      ; ldy=ny      ; ldz=nz
1354        call fftpad(rho,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat1,mgfft,-1,gbound)
1355      else
1356 
1357        call fftbox_plan3_many(plan,ndat1,ngfft(1:3),ngfft(1:3),ngfft(7),-1)
1358        call fftbox_execute(plan,rho)
1359      end if
1360 
1361      do ig=1,npwvec       ! Have to map FFT to G-sphere.
1362        igfft=igfftg0(ig)
1363        if (igfft/=0) then ! G-G0 belong to the FFT mesh.
1364          rhotwg(ig)=rhotwg(ig)+rho(igfft)
1365        end if
1366      end do
1367 
1368    CASE DEFAULT
1369      MSG_BUG("Wrong map2sphere")
1370    END SELECT
1371 
1372    RETURN
1373 
1374  CASE (2) ! Spinorial case.
1375 
1376    isprot1=spinrot1(1); isprot2=spinrot2(1) ! This is to bypass abirule
1377    MSG_ERROR("Spinorial case not implemented yet")
1378 
1379    SELECT CASE (map2sphere)
1380 
1381    CASE (0) ! Need results on the full FFT box thus cannot use zero-padded FFT.
1382    CASE (1) ! Need results on the G-sphere. Call zero-padded FFT routines if required.
1383    CASE DEFAULT
1384      MSG_BUG("Wrong map2sphere")
1385    END SELECT
1386 
1387    RETURN
1388 
1389  CASE DEFAULT
1390    MSG_BUG('Wrong nspinor')
1391  END SELECT
1392 
1393 end subroutine paw_cross_rho_tw_g

m_pawpwij/paw_mkrhox [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

 paw_mkrhox

FUNCTION

  Evaluate $<phj|e^{-i(q+G)}|phi>-<tphj|e^{-i(q+G)}|tphi>$
  for a fixed q-point and npw G vectors. Matrix elements are stored in packed storage mode.

INPUTS

  gmet(3,3)=reciprocal lattice metric tensor ($\textrm{Bohr}^{-2}$)
  gvec(3,npw)=G vectors in reduced coordinates
  npw=numper of G vectors
  Psps<pseudopotential_type>:
     %lnmax=Max. number of (l,n) components over all type of PAW datasets
  nq_spl=Number of points in the reciprocal space grid on which the radial functions pwff_spl are specified
  qgrid_spl(nq_spl)=values at which form factors have been evaluated
     %mpsang=1+maximum angular momentum
  qpt(3)= q-point in reduced coordinates
  ylm_q(npw,(2*Psps%mpsang-1)**2)=real spherical harmonics Ylm(q+G) for q-point qpt up to l=2*l_max
  pwff_spl(nq_spl,2,0:2*(Psps%mpsang-1),Psps%lnmax*(Psps%lnmax+1)/2))
  Pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
     %indklmn(8,lmn2_size)=array giving klm, kln, abs(il-jl) and (il+jl), ilm and jlm, ilmn and jlmn for each klmn=(ilmn,jlmn)

OUTPUT

  paw_rhox(2,npw,lmn2_size): $<phj|e^{-i(q+G).r}|phi>-<tphj|e^{-i(q+G).r}|tphi>$ in packed form for
    a type itypat (phase factor arising from atom position is not included)

PARENTS

      m_pawpwij

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

 922 subroutine paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,qgrid_spl,pwff_spl,&
 923 &  gmet,qpt,npw,gvec,ylm_q,Psps,Pawtab,paw_rhox)
 924 
 925 
 926 !This section has been created automatically by the script Abilint (TD).
 927 !Do not modify the following lines by hand.
 928 #undef ABI_FUNC
 929 #define ABI_FUNC 'paw_mkrhox'
 930 !End of the abilint section
 931 
 932  implicit none
 933 
 934 !Arguments ------------------------------------
 935 !scalars
 936  integer,intent(in) :: itypat,dim1,dim2,method,npw,nq_spl,lmn2_size
 937  type(Pseudopotential_type),intent(in) :: Psps
 938 !arrays
 939  integer,intent(in) :: gvec(3,npw)
 940  real(dp),intent(in) :: gmet(3,3)
 941  real(dp),intent(in) :: pwff_spl(nq_spl,2,0:dim1,dim2)
 942  real(dp),intent(in) :: qpt(3),ylm_q(npw,(2*Psps%mpsang-1)**2)
 943  real(dp),intent(out) :: paw_rhox(2,npw,lmn2_size)
 944  real(dp),intent(in) :: qgrid_spl(nq_spl)
 945  type(Pawtab_type),target,intent(in) :: Pawtab(Psps%ntypat)
 946 
 947 !Local variables-------------------------------
 948 !scalars
 949  integer :: ider,ig,ignt,il,ilm,ilm_G,ilmn,iln,im,ipow,jl,jlm
 950  integer :: jlmn,jln,jm,k0lm,k0lmn,k0ln,klm,klmn,kln,ll_G,mm_G,mpsang,ngnt
 951  real(dp) :: rgnt,dummy
 952  character(len=500) :: msg
 953 !arrays
 954  integer,allocatable :: gntselect(:,:)
 955  integer,ABI_CONTIGUOUS pointer :: indlmn(:,:)
 956  real(dp) :: mi_l(2,0:3),qpg(3)
 957  real(dp),allocatable :: derfun(:),newfun(:),qpg_norm(:),realgnt(:),wk_ffnl(:,:)
 958 
 959 ! *************************************************************************
 960 
 961  !write(std_out,*)itypat,dim1,dim2,method,npw,nq_spl,lmn2_size,Psps%mpsang
 962 
 963  mpsang=Psps%mpsang
 964  indlmn => Pawtab(itypat)%indlmn
 965  !
 966  ! === Pre-calculate (-i)^l ===
 967  mi_l(1,0)=one  ; mi_l(2,0)=zero
 968  mi_l(1,1)=zero ; mi_l(2,1)=-one
 969  mi_l(1,2)=-one ; mi_l(2,2)=zero
 970  mi_l(1,3)=zero ; mi_l(2,3)=one
 971  !
 972  ! === Calculate |q+G| ===
 973  ! * 2\pi is not included to be consistent with the spline.
 974  ABI_MALLOC(qpg_norm,(npw))
 975  do ig=1,npw
 976    qpg = qpt + gvec(:,ig)
 977    qpg_norm(ig)=SQRT(DOT_PRODUCT(qpg,MATMUL(gmet,qpg)))
 978  end do
 979 
 980  ! Check q-grid as %qgrid_spl must be large enoung.
 981  if (MAXVAL(qpg_norm)>MAXVAL(qgrid_spl)) then
 982    write(msg,'(3a,f8.4,a,f8.4,2a)')&
 983 &    ' Function values are being requested outside range of data. ',ch10,&
 984 &    ' Max qpg_norm = ',MAXVAL(qpg_norm),' Max qgrid_spl = ',MAXVAL(qgrid_spl),ch10,&
 985 &    ' Increase ecut(wfn), check qrid_ff and gsqcut '
 986    MSG_ERROR(msg)
 987  end if
 988 
 989  ABI_MALLOC(wk_ffnl,(nq_spl,2))
 990  ABI_MALLOC(newfun,(npw))
 991  ABI_MALLOC(derfun,(npw))
 992 
 993  SELECT CASE (method)
 994 
 995  CASE (PWIJ_ARNAUD)
 996    ! === Arnaud-Alouani exact expression ===
 997    ! * It does not describe the multipoles of the AE charge density
 998    ! * $ 4\pi \sum_{LM} (-i)^l Y_M^L(q+G) G_{\li\mi\lj\mj}^{\LM} ff^{aL}_{ij}(|q+G|) $
 999    !   where f has been calculated in paw_mkrhox_spl
1000    !
1001    ! === Re-evaluate Gaunt coefficients, just to be on the safe side ===
1002    ! * Note that gntselect is in packed form, thanks to invariance under permutation.
1003    ! * Could use Pawang% but size of gntselect depends on pawxcdev!
1004 
1005    ABI_MALLOC(  realgnt,((2*mpsang-1)**2*(mpsang)**4))
1006    ABI_MALLOC(gntselect,((2*mpsang-1)**2, mpsang**2*(mpsang**2+1)/2))
1007    call realgaunt(mpsang,ngnt,gntselect,realgnt)
1008    paw_rhox=zero
1009    !
1010    ! === Loop on (jl,jm,jn) channels for this atom ===
1011    do jlmn=1,Pawtab(itypat)%lmn_size
1012      jl =indlmn(1,jlmn)
1013      jm =indlmn(2,jlmn)
1014      jlm=indlmn(4,jlmn)
1015      jln=indlmn(5,jlmn)
1016 
1017      k0lmn=jlmn*(jlmn-1)/2
1018      k0lm =jlm *(jlm -1)/2
1019      k0ln =jln *(jln -1)/2
1020      !
1021      ! === Loop on (il,im,in) channels; klmn is index for packed form ===
1022      do ilmn=1,jlmn
1023        il =indlmn(1,ilmn)
1024        im =indlmn(2,ilmn)
1025        ilm=indlmn(4,ilmn)
1026        iln=indlmn(5,ilmn)
1027 
1028        klmn=k0lmn+ilmn
1029        klm =k0lm +ilm
1030        kln =k0ln +iln
1031        !
1032        ! === Summing over allowed (L,M), taking into account Gaunt selection rules ===
1033        do ll_G=ABS(jl-il),jl+il,2
1034          ipow=MOD(ll_G,4)
1035          ider=0
1036          wk_ffnl(:,:)=pwff_spl(:,:,ll_G,kln)
1037          call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw)
1038          do mm_G=-ll_G,ll_G
1039            ilm_G=1+ll_G**2+ll_G+mm_G
1040            ignt=gntselect(ilm_G,klm)
1041            if (ignt==0) CYCLE
1042            rgnt=realgnt(ignt)
1043            !
1044            ! === Evaluate matrix elements for each plane wave ===
1045            do ig=1,npw
1046              dummy = newfun(ig) * ylm_q(ig,ilm_G) * rgnt
1047              paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) &
1048 &               + ( dummy * mi_l(1,ipow) )
1049 
1050              paw_rhox(2,ig,klmn) = paw_rhox(2,ig,klmn) &
1051 &               + ( dummy * mi_l(2,ipow) )
1052            end do
1053          end do !mm_G
1054        end do !ll_G
1055 
1056      end do !ilmn
1057    end do !jlmn
1058 
1059    !
1060    ! * Multiply by 4\pi arising from the expansion of the plane wave
1061    paw_rhox = four_pi*paw_rhox
1062    ABI_FREE(realgnt)
1063    ABI_FREE(gntselect)
1064 
1065  CASE (PWIJ_SHISHKIN)
1066    ! === Shishkin-Kresse approximated expression ====
1067    ! * Better description of multipoles of AE charge,
1068    ! * Better results for energy degeneracies in GW band structure
1069    ! * $4\pi \sum_{LM} q_ij^{LM} Y_M^L(q+G) f^{aL}_{ij}(q+G)$ where f has been calculated in paw_mkrhox_spl
1070    !
1071    paw_rhox=zero
1072    !
1073    ! === Loop on (jl,jm,jn) channels for this atom ===
1074    !itypat=Cryst%typat(iatm)
1075    do jlmn=1,Pawtab(itypat)%lmn_size
1076      jl =indlmn(1,jlmn)
1077      jm =indlmn(2,jlmn)
1078      jlm=indlmn(4,jlmn)
1079      jln=indlmn(5,jlmn)
1080 
1081      k0lmn=jlmn*(jlmn-1)/2
1082      k0lm =jlm *(jlm -1)/2
1083      k0ln =jln *(jln -1)/2
1084      !
1085      ! === Loop on (il,im,in) channels; klmn is index for packed form ===
1086      do ilmn=1,jlmn
1087        il =indlmn(1,ilmn)
1088        im =indlmn(2,ilmn)
1089        ilm=indlmn(4,ilmn)
1090        iln=indlmn(5,ilmn)
1091 
1092        klmn=k0lmn+ilmn
1093        klm =k0lm +ilm
1094        kln =k0ln +iln
1095        !
1096        ! === Summing over allowed (l,m), taking into account Gaunt selection rules ===
1097        do ll_G=ABS(jl-il),jl+il,2
1098          ipow=MOD(ll_G,4)
1099          do mm_G=-ll_G,ll_G
1100            ! here I can move splfit before the loop over mm_G but I have to change paw_rhox_spl
1101            ilm_G=1+ll_G**2+ll_G+mm_G
1102            ider=0
1103            wk_ffnl(:,:)=pwff_spl(:,:,ilm_G-1,klmn)  ! Note klmn and ilm_G-1
1104            call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw)
1105            !
1106            ! === Evaluate matrix elements for each plane wave ===
1107            do ig=1,npw
1108              dummy = newfun(ig) * ylm_q(ig,ilm_G)
1109              paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) &
1110 &              + dummy * mi_l(1,ipow) !(ph3d(1,ig)*mi_l(1,ipow)-ph3d(2,ig)*mi_l(2,ipow))
1111 
1112              paw_rhox(2,ig,klmn) = paw_rhox(2,ig,klmn) &
1113 &              + dummy * mi_l(2,ipow) !(ph3d(1,ig)*mi_l(2,ipow)+ph3d(2,ig)*mi_l(1,ipow))
1114            end do
1115          end do !mm_G
1116        end do !ll_G
1117 
1118      end do !ilmn
1119    end do !jlmn
1120    !
1121    ! * Multiply by 4\pi arising from the expansion of the plane wave
1122    paw_rhox=four_pi*paw_rhox
1123 
1124  CASE DEFAULT
1125    write(msg,'(a,i3)')' Wrong value for method= ',method
1126    MSG_BUG(msg)
1127  END SELECT
1128 
1129  ABI_FREE(wk_ffnl)
1130  ABI_FREE(newfun)
1131  ABI_FREE(derfun)
1132  ABI_FREE(qpg_norm)
1133 
1134 end subroutine paw_mkrhox

m_pawpwij/paw_mkrhox_spl [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

 paw_mkrhox_spl

FUNCTION

  Evaluate PAW form factor ff^{aL}_{ij}(q) for each angular momentum L,
  each type of atom, a, and each [(in,il),(jn,jl)] channel. These quantities
  are used in paw_mkrhox to evaluate $<phi|e^{-i(q+G)}|phj>-<tphi|e^{-i(q+G)}|tphj>$
  for an arbitrary q+G vector.

INPUTS

  dim1
   = 2*(Pawtab(itypat)%l_size-1) if method=1
   = MAXVAL(Pawtab(:)%l_size)**2 if method=2
  dim2
   = Pawtab(itypat)%ij_size      if method=1
   = MAXVAL(Pawtab(:)%lmn2_size) if method=2
  method=integer flag defining the approach used:
   1 --> Expression based on the expansion on a plane wave in terms of Bessel functions
        and spherical harmonics (Arnaud-Alouani's methos, see PRB 62, 4464 [[cite:Arnaud2000]]
   2 --> Approximate expression with correct description of the multipoles. Eq. 9 in PRB 74, 035101 [[cite:Shishkin2006]]
  nq_spl=number of grid points in the q-mesh
  qgrid_spl(nq_spl)=values where form factors are returned
  ntypat=number of type of atoms
  Pawrad<type(Pawrad_type)>=datatype containing radial grid information
  Pawtab(ntypat)<type(pawtab_type)>=PAW tabulated starting data

OUTPUT

  pwff_spl(nq_spl,2,0:dim1,dim1_rhox2,ntypat)
   form factors ff^{aL}_{ij}(q) and second derivative in packed storage mode
  === if method=1 ===
    $ff_^{aL}_{ij}(q) =
      \int_0^{r_a} j_L(2\pi qr) [phi_{n_i,l_i}.phi_{n_j l_j}(r) - tphi_{n_i l_i}.tph_{n_j l_j}(r)] dr$
  === if method=2 ===
    $ff_^{aL}_{ij}(q) = q_ij^{LM} \int_0^{r_a} j_L(2\pi qr) g_L(r) r^2 dr$

NOTES

  * $j_L(2\pi q)$ is a spherical Bessel function
  * Output matrix elements are stored in packed storage mode
  * Inspired by pawpsp_nl

TODO

  One might save CPU time taking into account Gaunt selection rules!

PARENTS

      m_pawpwij

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

604 subroutine paw_mkrhox_spl(itypat,ntypat,method,dim1,dim2,nq_spl,qgrid_spl,Pawrad,Pawtab,pwff_spl)
605 
606 
607 !This section has been created automatically by the script Abilint (TD).
608 !Do not modify the following lines by hand.
609 #undef ABI_FUNC
610 #define ABI_FUNC 'paw_mkrhox_spl'
611 !End of the abilint section
612 
613  implicit none
614 
615 !Arguments ------------------------------------
616 !scalars
617  integer,intent(in) :: itypat,ntypat,method,dim2,dim1,nq_spl
618 !arrays
619  real(dp),intent(in) :: qgrid_spl(nq_spl)
620  real(dp),intent(out) :: pwff_spl(nq_spl,2,0:dim1,dim2)
621  type(Pawrad_type),intent(in) :: Pawrad(ntypat)
622  type(Pawtab_type),intent(in) :: Pawtab(ntypat)
623 
624 !Local variables-------------------------------
625 !scalars
626  integer :: mm,nlmn,jlmn,ilmn,klmn,ij_size,l_size,ider
627  integer :: iq,ir,ll,meshsz,mmax,iln,jln,nln,k0ln,kln,qlm
628  real(dp),parameter :: EPS=tol14**4,TOLJ=0.001_dp
629  real(dp) :: arg,argn,bes,besp,qr,yp1,ypn
630  character(len=500) :: msg
631  type(Pawrad_type) :: Tmpmesh
632 !arrays
633  real(dp),allocatable :: rrshape_l(:),shape_l(:),ff(:),gg(:),rr(:),rrdphi_ij(:)
634  real(dp),allocatable :: dphi_ij(:),tmp_spl(:,:,:,:),tmp_jgl(:,:,:)
635 
636 !*************************************************************************
637 
638  DBG_ENTER("COLL")
639 
640  pwff_spl=zero
641 
642  SELECT CASE (method)
643 
644  CASE (PWIJ_ARNAUD)
645    ! === Arnaud-Alouani exact expression PRB 62. 4464 [[cite:Arnaud2000]] ===
646    ! * $ff_^{aL}_{ij}(q) =
647    !    \int_0^{r_a} j_L(2\pi qr) [phi_{n_i,l_i}.phi_{n_j l_j}(r)-tphi_{n_i l_i}.tph_{n_j l_j}(r)]dr$
648    ! * It does not descrive correctly the multipoles of the AE charge density if low cutoff on G
649    write(msg,'(a,i3)')' paw_mkrhox_spl: Using Arnaud-Alouani expression for atom type: ',itypat
650    call wrtout(std_out,msg,'COLL')
651 
652    nln     = Pawtab(itypat)%basis_size
653    ij_size = Pawtab(itypat)%ij_size
654    l_size  = Pawtab(itypat)%l_size
655 
656    ABI_MALLOC(tmp_spl,(nq_spl,2,0:l_size-1,ij_size))
657 
658    ! Is mesh beginning with r=0 ?
659    if (ABS(Pawrad(itypat)%rad(1))>tol10) then
660      MSG_ERROR("Radial mesh starts with r/=0")
661    end if
662    !
663    ! === Initialize temporary arrays and variables ===
664    meshsz = Pawtab(itypat)%mesh_size ; mmax=meshsz
665 
666    ABI_MALLOC(dphi_ij,(meshsz))
667    ABI_MALLOC(rrdphi_ij,(meshsz))
668    ABI_MALLOC(ff,(meshsz))
669    ABI_CALLOC(gg,(meshsz))
670    ABI_MALLOC(rr,(meshsz))
671    rr(1:meshsz) = Pawrad(itypat)%rad(1:meshsz)
672    !
673    ! === Loop on (jln,iln) channels for this type. Packed form ===
674    tmp_spl=zero
675 
676    do jln=1,nln
677      k0ln=jln*(jln-1)/2
678      do iln=1,jln
679        kln=k0ln+iln
680 
681        dphi_ij(1:meshsz) = Pawtab(itypat)%phiphj(1:meshsz,kln)-Pawtab(itypat)%tphitphj(1:meshsz,kln)
682        rrdphi_ij(1:meshsz) = rr(1:meshsz)*dphi_ij(1:meshsz) ! will be used for first derivative
683 
684        ir=meshsz
685        do while (ABS(dphi_ij(ir))<EPS)
686          ir=ir-1
687        end do
688        mmax=MIN(ir+1,meshsz)
689        ! ir is equal to meshsz if no point are below EPS
690        if (mmax/=Pawrad(itypat)%int_meshsz) then  ! mmax=meshsz
691          call pawrad_init(Tmpmesh,mesh_size=Pawtab(itypat)%mesh_size,mesh_type=Pawrad(itypat)%mesh_type, &
692 &                         rstep=Pawrad(itypat)%rstep,lstep=Pawrad(itypat)%lstep,r_for_intg=rr(mmax))
693        else
694          call pawrad_copy(Pawrad(itypat),Tmpmesh)
695        end if
696        !
697        ! === Loop on l for Bessel function. Note the starting point ===
698        ! TODO Here I should loop only the moments allowed by Gaunt, I should use indklm!
699        ! and only on lmax for this atom
700        do ll=0,l_size-1
701          !
702          ! === Compute f_l(q=0) only if l=0, and first derivative fp_l(q=0) (nonzero only if ll==1) ===
703          tmp_spl(1,1,ll,kln)=zero; yp1=zero
704          if (ll==0) then
705            call simp_gen(tmp_spl(1,1,ll,kln),dphi_ij,Tmpmesh)
706          end if
707          if (ll==1) then
708            call simp_gen(yp1,rrdphi_ij,Tmpmesh)
709            yp1=yp1*two_pi*third
710          end if
711          !
712          ! === Compute f_l(0<q<qmax) ===
713          if (nq_spl>2) then
714            do iq=2,nq_spl-1
715              arg=two_pi*qgrid_spl(iq)
716              do ir=1,mmax
717                qr=arg*rr(ir)
718                call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ)
719                ff(ir)=bes*dphi_ij(ir)
720              end do
721              call simp_gen(tmp_spl(iq,1,ll,kln),ff,Tmpmesh)
722            end do
723          end if
724          !
725          ! === Compute f_l(q=qmax) and first derivative ===
726          if (nq_spl>1) then
727            argn=two_pi*qgrid_spl(nq_spl)
728            do ir=1,mmax
729              qr=argn*rr(ir)
730              call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ)
731              ff(ir)=bes *  dphi_ij(ir)
732              gg(ir)=besp*rrdphi_ij(ir)
733            end do
734            call simp_gen(tmp_spl(nq_spl,1,ll,kln),ff,Tmpmesh)
735            gg(:)=two_pi*gg(:) ! two_pi comes from 2\pi|q| in the Bessel function
736            call simp_gen(ypn,gg,Tmpmesh)
737          else
738            ypn=yp1
739          end if
740          !
741          ! === Compute second derivative of ff^{al}_{ij)(q) ===
742          !yp1=zero; ypn=zero
743          call paw_spline(qgrid_spl,tmp_spl(:,1,ll,kln),nq_spl,yp1,ypn,tmp_spl(:,2,ll,kln))
744        end do !ll
745 
746        call pawrad_free(Tmpmesh)
747 
748      end do !iln
749    end do !jln
750    !
751    ! === Save values for this atom type, each ll and kln channel ===
752    pwff_spl = tmp_spl
753 
754    ABI_FREE(dphi_ij)
755    ABI_FREE(rrdphi_ij)
756    ABI_FREE(ff)
757    ABI_FREE(gg)
758    ABI_FREE(rr)
759    ABI_FREE(tmp_spl)
760 
761    if (.FALSE.) then ! write form factors on file for plotting purpose.
762      ll=0
763      do iq=1,nq_spl
764        write(777+itypat,'(50(es16.8))')qgrid_spl(iq),((pwff_spl(iq,ider,ll,iln),ider=1,2),iln=1,dim2)
765      end do
766    end if
767 
768  CASE (PWIJ_SHISHKIN)
769    ! ==== Shishkin-Kresse approximated expression ====
770    ! $ff_^{aL}_{ij}(q) = q_ij^{LM} \int_0^{r_a} j_L(2\pi qr) g_L(r) r^2 dr$
771    ! * Better description of multipoles of AE charge,
772    ! * Better results for energy degeneracies in GW band structure
773    write(msg,'(a,i3)')' paw_mkrhox_spl: Using Shishkin-Kresse expression for atom type ',itypat
774    call wrtout(std_out,msg,'COLL')
775    l_size = Pawtab(itypat)%l_size
776    nlmn   = Pawtab(itypat)%lmn_size
777 
778    !allocate(tmp_jgl(nq_spl,2,0:2*(Psps%mpsang-1)))
779    ABI_MALLOC(tmp_jgl,(nq_spl,2,0:l_size-1))
780 
781    ! Is mesh beginning with r=0 ?
782    if (ABS(Pawrad(itypat)%rad(1))>tol10) then
783      MSG_ERROR("Radial mesh starts with r/=0")
784    end if
785    !
786    ! === Initialize temporary arrays and variables ===
787 
788    call pawrad_copy(Pawrad(itypat),Tmpmesh)
789    meshsz=Pawtab(itypat)%mesh_size ; mmax=meshsz
790    ABI_MALLOC(ff,(meshsz))
791    ABI_CALLOC(gg,(meshsz))
792    ABI_MALLOC(rr,(meshsz))
793    ABI_MALLOC(shape_l,(meshsz))
794    ABI_MALLOC(rrshape_l,(meshsz))
795    rr(1:meshsz)=Tmpmesh%rad(1:meshsz)
796 
797    tmp_jgl(:,:,:)=zero
798    rrshape_l(1)=zero
799      shape_l(1)=zero
800    !
801    ! TODO Here I should loop only the moments allowed by Gaunt, I should use indklm!
802    ! and only lmax for this atom
803    !do ll=0,2*(Psps%mpsang-1)
804    do ll=0,l_size-1
805        shape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**2
806      rrshape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**3
807      !
808      ! === Compute f_l(q=0) and first derivative fp_l(q=0) (only if ll==1) ===
809      tmp_jgl(1,1,ll)=zero ; yp1=zero
810      if (ll==0) then
811        call simp_gen(tmp_jgl(1,1,ll),shape_l,Tmpmesh)
812      end if
813      if (ll==1) then
814        call simp_gen(yp1,rrshape_l,Tmpmesh) !rr comes from d/dq
815        yp1=yp1*two_pi*third
816      end if
817      !
818      ! === Compute f_l(0<q<qmax) ===
819      if (nq_spl>2) then
820        do iq=2,nq_spl-1
821          arg=two_pi*qgrid_spl(iq)
822          do ir=1,mmax
823            qr=arg*rr(ir)
824            call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ)
825            ff(ir)=bes*shape_l(ir)
826          end do
827          call simp_gen(tmp_jgl(iq,1,ll),ff,Tmpmesh)
828        end do
829      end if
830      !
831      ! === Compute f_l(q=qmax) and first derivative ===
832      if (nq_spl>1) then
833        argn=two_pi*qgrid_spl(nq_spl)
834        do ir=1,mmax
835          qr=argn*rr(ir)
836          call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ)
837          ff(ir)=bes *  shape_l(ir)
838          gg(ir)=besp*rrshape_l(ir)
839        end do
840        call simp_gen(tmp_jgl(nq_spl,1,ll),ff,Tmpmesh)
841        gg(:)=two_pi*gg(:) !two_pi comes from 2\pi|q|
842        call simp_gen(ypn,gg,Tmpmesh)
843      else
844        ypn=yp1
845      end if
846      !
847      ! === Compute second derivative of ff_^{al}_{ij)(q) ===
848      call paw_spline(qgrid_spl,tmp_jgl(:,1,ll),nq_spl,yp1,ypn,tmp_jgl(:,2,ll))
849    end do !ll
850    !
851    ! === Save values for this type, each ll and ilmn,jlm channels ===
852    ! * Here we assembly q_{ij}^{lm} \int_0^{r_a} j_l(2\pi(q+G)r) g_l(r)r^2 dr
853    ! * Some of the contributions from qijl are zero due to Gaunt selection rules
854    do jlmn=1,nlmn
855      do ilmn=1,jlmn
856        klmn=ilmn+(jlmn-1)*jlmn/2
857        !do ll=0,2*(Psps%mpsang-1)
858        do ll=0,l_size-1
859          do mm=-ll,ll
860            qlm=1+ll**2+ll+mm
861            pwff_spl(:,:,qlm-1,klmn)=tmp_jgl(:,:,ll)*Pawtab(itypat)%qijl(qlm,klmn)
862          end do
863        end do
864      end do
865    end do
866 
867    call pawrad_free(Tmpmesh)
868    ABI_FREE(shape_l)
869    ABI_FREE(rrshape_l)
870    ABI_FREE(ff)
871    ABI_FREE(gg)
872    ABI_FREE(rr)
873    ABI_FREE(tmp_jgl)
874 
875  CASE DEFAULT
876    write(msg,'(a,i3)')' Called with wrong value for method ',method
877    MSG_BUG(msg)
878  END SELECT
879 
880  DBG_EXIT("COLL")
881 
882 end subroutine paw_mkrhox_spl

m_pawpwij/paw_rho_tw_g [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

 paw_rho_tw_g

FUNCTION

  Evaluates the PAW onsite contribution to the oscillator strengths:
  sum_{i,j} <\tpsi_{k-q,b1}|\cprj_i> <\cprj_j|\tpsi_{k,b2}>*
   \[ <\phi_i|e^{-i(q+G).r}|\phi_j> - <\tilde\phi_i|e^{-i(q+G).r}|\tilde\phi_j> \].

INPUTS

 dim_rtwg=Define the size of the array rhotwg
   === for nspinor==1 ===
    dim_rtwg=1
   === for nspinor==2 ===
    dim_rtwg=2 if only <up|up>, <dwn|dwn> matrix elements are required
    dim_rtwg=4 for <up|up>, <dwn|dwn>, <up|dwn> and <dwn|up>.
  nspinor=number of spinorial components.
  npw=number of plane waves for oscillator matrix elements
  natom=number of atoms
  Cprj_kmqb1(natom,nspinor),Cprj_kb2(natom,nspinor) <type(pawcprj_type)>=
   projected input wave functions <Proj_i|Cnk> with all NL projectors corresponding to
   wavefunctions (k-q,b1,s) and (k,b2,s), respectively.

SIDE EFFECTS

  rhotwg(npw*dim_rtwg)=Updated oscillator strengths with the on-site PAW contributions added.

PARENTS

      calc_sigc_me,calc_sigx_me,cchi0,cchi0q0,cchi0q0_intraband
      check_completeness,cohsex_me,exc_build_block,exc_build_ham,m_shirley
      prep_calc_ucrpa

SOURCE

1173 pure subroutine paw_rho_tw_g(npw,dim_rtwg,nspinor,natom,ntypat,typat,xred,gvec,Cprj_kmqb1,Cprj_kb2,Pwij,rhotwg)
1174 
1175 
1176 !This section has been created automatically by the script Abilint (TD).
1177 !Do not modify the following lines by hand.
1178 #undef ABI_FUNC
1179 #define ABI_FUNC 'paw_rho_tw_g'
1180 !End of the abilint section
1181 
1182  implicit none
1183 
1184 !Arguments ------------------------------------
1185 !scalars
1186  integer,intent(in) :: natom,ntypat,npw,nspinor,dim_rtwg
1187 !arrays
1188  integer,intent(in) :: gvec(3,npw),typat(natom)
1189  real(dp),intent(in) :: xred(3,natom)
1190  complex(gwpc),intent(inout) :: rhotwg(npw*dim_rtwg)
1191  type(pawcprj_type),intent(in) :: Cprj_kmqb1(natom,nspinor),Cprj_kb2(natom,nspinor)
1192  type(pawpwij_t),intent(in) :: Pwij(ntypat)
1193 
1194 !Local variables-------------------------------
1195 !scalars
1196  integer :: ig,iat,nlmn,ilmn,jlmn,k0lmn,klmn,iab,isp1,isp2,spad,itypat
1197  real(dp) :: fij,re_psp,im_psp,re_pw,im_pw,arg
1198 !arrays
1199  integer,parameter :: spinor_idxs(2,4)=RESHAPE((/1,1,2,2,1,2,2,1/),(/2,4/))
1200  real(dp) :: tmp(2),qpg(3),x0(3),ph3d(2)
1201 
1202 ! *************************************************************************
1203 
1204  ! === Loop over the four spinorial combinations ===
1205  do iab=1,dim_rtwg
1206    isp1=spinor_idxs(1,iab)
1207    isp2=spinor_idxs(2,iab)
1208    spad=npw*(iab-1)
1209 
1210    do ig=1,npw
1211      tmp(:)=zero
1212      do iat=1,natom
1213        itypat = typat(iat)
1214        nlmn   = Pwij(itypat)%lmn_size
1215        x0(:) = xred(:,iat)
1216 
1217        ! === Structure factor e^{-i(q+G)*xred} ===
1218        qpg(:)= Pwij(itypat)%qpt(:) + gvec(:,ig)
1219        arg=-two_pi*DOT_PRODUCT(qpg(:),x0)
1220        ph3d(1)=COS(arg)
1221        ph3d(2)=SIN(arg)
1222 
1223        ! === Loop on [(jl,jm,jn),(il,im,in)] channels. packed storage mode ===
1224        do jlmn=1,nlmn
1225          k0lmn=jlmn*(jlmn-1)/2
1226          do ilmn=1,jlmn
1227            re_psp =  Cprj_kmqb1(iat,isp1)%cp(1,ilmn) * Cprj_kb2(iat,isp2)%cp(1,jlmn) &
1228 &                   +Cprj_kmqb1(iat,isp1)%cp(2,ilmn) * Cprj_kb2(iat,isp2)%cp(2,jlmn) &
1229 &                   +Cprj_kmqb1(iat,isp1)%cp(1,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn) &
1230 &                   +Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(2,ilmn)
1231 
1232            im_psp =  Cprj_kmqb1(iat,isp1)%cp(1,ilmn) * Cprj_kb2(iat,isp2)%cp(2,jlmn) &
1233 &                   -Cprj_kmqb1(iat,isp1)%cp(2,ilmn) * Cprj_kb2(iat,isp2)%cp(1,jlmn) &
1234 &                   +Cprj_kmqb1(iat,isp1)%cp(1,jlmn) * Cprj_kb2(iat,isp2)%cp(2,ilmn) &
1235 &                   -Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn)
1236 
1237            klmn=k0lmn+ilmn; fij=one; if (jlmn==ilmn) fij=half
1238 
1239            ! Multiply by the phase due to the atom position.
1240            re_pw =  Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(1) &
1241                    -Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(2)
1242 
1243            im_pw =  Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(2) &
1244                    +Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(1)
1245 
1246            tmp(1)=tmp(1)+ fij * (re_pw*re_psp - im_pw*im_psp)
1247            tmp(2)=tmp(2)+ fij * (re_pw*im_psp + im_pw*re_psp)
1248          end do !ilmn
1249        end do !jlmn
1250      end do !iat
1251      !
1252      !if(ig==1) write(std_out,*) " TOTAL PW     osc str = ",rhotwg(ig+spad)
1253      !if(ig==1) write(std_out,*) " TOTAL PAW    osc str = ",tmp(1),tmp(2)
1254      ! Update input data using the appropriate index.
1255      rhotwg(ig+spad) = rhotwg(ig+spad) + CMPLX(tmp(1),tmp(2),kind=gwpc)
1256      !if(ig==1) write(std_out,*) " TOTAL PW+PAW osc str = ",rhotwg(ig+spad)
1257    end do !ig
1258 
1259  end do !dim_rtwg
1260 
1261 end subroutine paw_rho_tw_g

m_pawpwij/pawpwff_free [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

  pawpwff_free

FUNCTION

  Free the memory allocated in a structure of type pawpwff_t

SIDE EFFECTS

  Paw_pwff(:)=<pawpwff_t>=Object storing form factors for the spline of wf into PAW spheres

PARENTS

      bethe_salpeter,m_shirley,screening,sigma

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

272 subroutine pawpwff_free(Paw_pwff)
273 
274 
275 !This section has been created automatically by the script Abilint (TD).
276 !Do not modify the following lines by hand.
277 #undef ABI_FUNC
278 #define ABI_FUNC 'pawpwff_free'
279 !End of the abilint section
280 
281  implicit none
282 
283 !Arguments ------------------------------------
284 !scalars
285  type(pawpwff_t),intent(inout) :: Paw_pwff(:)
286 
287 !Local variables-------------------------------
288 !scalars
289  integer :: ii
290 
291 !************************************************************************
292 
293  !@pawpwff_t
294  do ii=1,SIZE(Paw_pwff)
295    if (allocated(Paw_pwff(ii)%qgrid_spl)) then
296      ABI_FREE(Paw_pwff(ii)%qgrid_spl)
297    end if
298    if (allocated(Paw_pwff(ii)%pwff_spl )) then
299      ABI_FREE(Paw_pwff(ii)%pwff_spl)
300    end if
301  end do
302 
303 end subroutine pawpwff_free

m_pawpwij/pawpwff_init [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

  pawpwff_init

FUNCTION

  Initialize the structure containing integrals used to evaluate the onsite
  matrix elements of a plane wave by means of a spline fit technique.

INPUTS

  method=1 for Arnaud-Alouani, 2 for Shishkin-Kresse.
  nq_spl(%ntypat)=Number of points in the mesh used for the spline.
  qmax(%ntypat)=Max |q| for the mesh
  gmet(3,3)=reciprocal space metric tensor in bohr**-2.
  Pawrad(%ntypat) <type(pawrad_type)>=paw radial mesh and related data....
  Pawtab(%ntypat) <type(pawtab_type)>=paw tabulated starting data.
    %lsize=1+maximum value of l leading to non zero Gaunt coeffs.
    %ij_size=Number of (i,j) elements for the symetric paw basis
    %lmn2_size=lmn_size*(lmn_size+1)/2
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %ntypat=Number of type of atoms.

OUTPUT

  Paw_pwff(%ntypat) <pawpwff_t>=Object storing the form factors
                                    for the spline used in pawpwij_init.

PARENTS

      bethe_salpeter,m_shirley,screening,sigma

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

180 subroutine pawpwff_init(Paw_pwff,method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps)
181 
182 
183 !This section has been created automatically by the script Abilint (TD).
184 !Do not modify the following lines by hand.
185 #undef ABI_FUNC
186 #define ABI_FUNC 'pawpwff_init'
187 !End of the abilint section
188 
189  implicit none
190 
191 !Arguments ------------------------------------
192 !scalars
193  integer,intent(in) :: method
194  type(Pseudopotential_type),intent(in) :: Psps
195 !arrays
196  integer,intent(in) :: nq_spl(Psps%ntypat)
197  real(dp),intent(in) :: gmet(3,3)
198  real(dp),intent(in) :: qmax(Psps%ntypat)
199  type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat)
200  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
201  type(pawpwff_t),intent(out) :: Paw_pwff(Psps%ntypat)
202 
203 !Local variables-------------------------------
204 !scalars
205  integer :: dim1,dim2,itypat,nq
206  real(dp) :: dq
207 
208 !************************************************************************
209 
210  !@pawpwff_t
211 
212 ! === Evaluate form factors for the radial part of phi.phj-tphi.tphj ===
213  do itypat=1,Psps%ntypat
214    Paw_pwff(itypat)%method = method
215 
216    select case (method)
217    case (PWIJ_ARNAUD)
218      dim1 = Pawtab(itypat)%l_size-1
219      dim2 = Pawtab(itypat)%ij_size
220    case (PWIJ_SHISHKIN)
221      dim1 = Pawtab(itypat)%l_size**2
222      dim2 = Pawtab(itypat)%lmn2_size
223    case default
224      MSG_BUG("Wrong method")
225    end select
226 
227    Paw_pwff(itypat)%dim1 = dim1
228    Paw_pwff(itypat)%dim2 = dim2
229 
230    Paw_pwff(itypat)%gmet = gmet
231 
232    ! === Setup of the q-mesh for spline ===
233    ! * It can be type-dependent.
234    nq = nq_spl(itypat)
235    dq = qmax(itypat)/(one*(nq-1))
236    !write(std_out,*)"nq,dq",nq,dq
237 
238    Paw_pwff(itypat)%nq_spl = nq
239    ABI_MALLOC(Paw_pwff(itypat)%qgrid_spl,(nq))
240    Paw_pwff(itypat)%qgrid_spl = arth(zero,dq,nq)
241    !
242    ! === Calculate form factors depending on method ===
243    ABI_MALLOC(Paw_pwff(itypat)%pwff_spl,(nq,2,0:dim1,dim2))
244 
245    call paw_mkrhox_spl(itypat,Psps%ntypat,method,dim1,dim2,nq,&
246 &    Paw_pwff(itypat)%qgrid_spl,Pawrad,Pawtab,Paw_pwff(itypat)%pwff_spl)
247  end do !itypat
248 
249 end subroutine pawpwff_init

m_pawpwij/pawpwff_t [ Types ]

[ Top ] [ m_pawpwij ] [ Types ]

NAME

 pawpwff_t

FUNCTION

 For PAW, form factors used to evaluate $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$

SOURCE

60  type,public :: pawpwff_t
61 
62   integer :: method
63   ! 1 For Arnaud-Alouani"s exact expression.
64   ! 2 For Shishkin-Kresse"s approximated expression.
65 
66   integer :: dim1
67   integer :: dim2
68   ! Dimensions of pwff_spl, depending on method.
69 
70   integer :: nq_spl
71    ! Number of points in the reciprocal space grid on which
72    ! the radial integrals are evaluated.
73 
74   real(dp) :: gmet(3,3)
75   ! Reciprocal space metric tensor in Bohr**-2
76 
77   real(dp),allocatable :: qgrid_spl(:)
78    ! qgrid_spl(nq_spl)
79    ! The coordinates of the points of the radial grid for the integrals used in the spline.
80 
81   real(dp),allocatable :: pwff_spl(:,:,:,:)
82   ! pwff_spl(nq_spl,2,0:dim1,dim2)
83   ! The different integrals on the radial |q| grid, for a given atom type.
84 
85  end type pawpwff_t
86 
87  public :: pawpwff_init    ! Initialize form factors for spline.
88  public :: pawpwff_free    ! Deallocate dynamic memory.

m_pawpwij/pawpwij_free_d1 [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

  pawpwij_free_d1

FUNCTION

  Free all memory allocated in a structure of type pawpwij_t

SIDE EFFECTS

  Paw_pwij(:)=<pawpwij_t>=Structure containing the onsite matrix elements of e^{-i(q+G).r}

PARENTS

      m_pawpwij

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

472 subroutine pawpwij_free_d1(Pwij)
473 
474 
475 !This section has been created automatically by the script Abilint (TD).
476 !Do not modify the following lines by hand.
477 #undef ABI_FUNC
478 #define ABI_FUNC 'pawpwij_free_d1'
479 !End of the abilint section
480 
481  implicit none
482 
483 !Arguments ------------------------------------
484 !scalars
485  type(pawpwij_t),intent(inout) :: Pwij(:)
486 
487 !Local variables-------------------------------
488 !scalars
489  integer :: ii
490 
491 !************************************************************************
492 
493  !@pawpwij_t
494  do ii=1,SIZE(Pwij)
495    if (allocated(Pwij(ii)%mqpgij))  then
496      ABI_FREE(Pwij(ii)%mqpgij)
497    end if
498  end do
499 
500 end subroutine pawpwij_free_d1

m_pawpwij/pawpwij_free_d2 [ Functions ]

[ Top ] [ m_pawpwij ] [ Functions ]

NAME

  pawpwij_free_d2

FUNCTION

  Free all memory allocated in a structure of type pawpwij_t

SIDE EFFECTS

  Paw_pwij(:)=<pawpwij_t>=Structure containing the onsite matrix elements of e^{-i(q+G).r}

PARENTS

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

522 subroutine pawpwij_free_d2(Pwij)
523 
524 
525 !This section has been created automatically by the script Abilint (TD).
526 !Do not modify the following lines by hand.
527 #undef ABI_FUNC
528 #define ABI_FUNC 'pawpwij_free_d2'
529 !End of the abilint section
530 
531  implicit none
532 
533 !Arguments ------------------------------------
534 !scalars
535  type(pawpwij_t),intent(inout) :: Pwij(:,:)
536 
537 !Local variables-------------------------------
538 !scalars
539  integer :: jj
540 
541 !************************************************************************
542 
543  do jj=1,SIZE(Pwij,DIM=2)
544    call pawpwij_free_d1(Pwij(:,jj))
545  end do
546 
547 end subroutine pawpwij_free_d2

m_pawpwij/pawpwij_t [ Types ]

[ Top ] [ m_pawpwij ] [ Types ]

NAME

 pawpwij_t

FUNCTION

 For PAW, object storing $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$
 for a given q-point, for a particular TYPE of atom. Therefore the phase factor
 e^{i(q+G).R_at} has to be considered to have the onsite contribution of a particular atom.

SOURCE

104  type,public :: pawpwij_t
105 
106   integer :: istpw
107   ! Storage mode (similar to istwfk), not used at present
108 
109   integer :: npw
110    ! The number of plane waves
111 
112   integer :: lmn_size
113    ! Number of (l,m,n) elements for the paw basis
114 
115   integer :: lmn2_size
116    ! lmn2_size=lmn_size*(lmn_size+1)/2
117    ! where lmn_size is the number of (l,m,n) elements for the paw basis
118 
119   real(dp) :: qpt(3)
120   ! The q-point in e^{-i(q+G)}.r}
121 
122   real(dp),allocatable :: mqpgij(:,:,:)
123   ! pwij(2,npw,lmn2_size)
124   ! $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$
125 
126  end type pawpwij_t
127 
128  public :: pawpwij_init       ! Calculate onsite matrix elements of a set of plane waves.
129  public :: pawpwij_free        ! Deallocate dynamic memory in the structure.
130  !public :: paw_pwij_bcast
131  public :: paw_rho_tw_g        ! Calculate the PAW contribution to the oscillator matrix element.
132  public :: paw_cross_rho_tw_g  ! Calculate the PAW cross term contribution to the oscillator matrix element.

m_pawrhox/pawpwij_init [ Functions ]

[ Top ] [ Functions ]

NAME

  pawpwij_init

FUNCTION

  Creation method for pawpwij_t. Calculates the onsite matrix elements
   $ <phj|e^{-i(q+G)}|phi> - <tphj|e^{-i(q+G)}|tphi> $
  for a given q and a set of G"s for a given __TYPE__ of atom.
  Phase factors arising from atom positions are therefore not included.

INPUTS

  npw=Number of plane waves
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %ntypat=Number of type of atoms,
  qpt_in(3)=The reduced components of the q-point.
  rprim(3,3)=dimensionless real space primitive translations
  Pawtab(%ntypat) <type(pawtab_type)>=paw tabulated starting data
  Paw_pwff(%ntypat) <pawpwff_t>=Object storing the form factors for the spline used in pawpwij_init.
  Psps <type(pseudopotential_type)>=variables related to pseudopotentials
    %ntypat

OUTPUT

  Pwij(%ntypat)<pawpwij_t>=Structure containing the onsite matrix elements of e^{-i(q+G).r}.
   Completely initialized in output.

PARENTS

      calc_sigc_me,calc_sigx_me,cchi0,cchi0q0,cchi0q0_intraband,cohsex_me
      exc_build_block,exc_build_ham,m_shirley,prep_calc_ucrpa

CHILDREN

      fftbox_execute,fftbox_plan3_many,fftpad

SOURCE

343 subroutine pawpwij_init(Pwij,npw,qpt_in,gvec,rprimd,Psps,Pawtab,Paw_pwff)
344 
345 
346 !This section has been created automatically by the script Abilint (TD).
347 !Do not modify the following lines by hand.
348 #undef ABI_FUNC
349 #define ABI_FUNC 'pawpwij_init'
350 !End of the abilint section
351 
352  implicit none
353 
354 !Arguments ------------------------------------
355 !scalars
356  integer,intent(in) :: npw
357  type(Pseudopotential_type),intent(in) :: Psps
358 !arrays
359  integer,intent(in) :: gvec(3,npw)
360  real(dp),intent(in) :: qpt_in(3),rprimd(3,3)
361  type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat)
362  type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat)
363  type(pawpwij_t),intent(out) :: Pwij(Psps%ntypat)
364 
365 !Local variables-------------------------------
366 !scalars
367  integer,parameter :: unkg0=0,unylm0=0
368  integer :: dim1,dim2,method
369  integer :: my_mqmem,my_nqpt,optder,two_lmaxp1,itypat
370  integer :: dummy_nsppol,lmn_size,lmn2_size,nq_spl,ierr
371  real(dp) :: ucvol
372  type(MPI_type) :: MPI_enreg_seq
373 !arrays
374  integer,allocatable :: npwarr(:),dummy_nband(:)
375  real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
376  real(dp),allocatable :: my_qtmp(:,:)
377  real(dp),allocatable :: ylm_q(:,:),ylmgr_q(:,:,:)
378 
379 ! *********************************************************************
380 
381  !@pawpwij_t
382 
383  ! ===============================================
384  ! === Get real spherical harmonics in G space ===
385  ! ===============================================
386 
387  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
388 
389  ! * Set up of REAL Ylm(q+G) up to 2*l_max for this q-point.
390  my_mqmem=1; two_lmaxp1=2*Psps%mpsang-1; optder=0
391 
392  ABI_MALLOC(ylm_q  ,(npw*my_mqmem,two_lmaxp1**2))
393  ABI_MALLOC(ylmgr_q,(npw*my_mqmem,3+6*(optder/2),two_lmaxp1**2))
394 
395  my_nqpt=1
396  ABI_MALLOC(my_qtmp,(3,my_nqpt))
397  my_qtmp(:,1)=qpt_in(:)
398  !
399  ! * dummy_nband and dummy_nsppol are not used in sequential mode.
400  dummy_nsppol=1
401  ABI_MALLOC(dummy_nband,(my_nqpt*dummy_nsppol))
402  dummy_nband=0
403  ABI_MALLOC(npwarr,(my_nqpt))
404  npwarr(:)=npw
405 
406  ! Fake MPI_type for sequential part.
407  call initmpi_seq(MPI_enreg_seq)
408 
409  call initylmg(gprimd,gvec,my_qtmp,my_mqmem,MPI_enreg_seq,two_lmaxp1,npw,dummy_nband,my_nqpt,&
410    npwarr,dummy_nsppol,optder,rprimd,ylm_q,ylmgr_q)
411 
412  call destroy_mpi_enreg(MPI_enreg_seq)
413 
414  ABI_FREE(my_qtmp)
415  ABI_FREE(dummy_nband)
416  ABI_FREE(npwarr)
417  ABI_FREE(ylmgr_q)
418 
419  ! Construct the Pwij structure
420  do itypat=1,Psps%ntypat
421 
422    Pwij(itypat)%istpw = 1
423    Pwij(itypat)%npw   = npw
424 
425    lmn_size  = Pawtab(itypat)%lmn_size
426    lmn2_size = lmn_size*(lmn_size+1)/2
427    Pwij(itypat)%lmn_size  = lmn_size
428    Pwij(itypat)%lmn2_size = lmn2_size
429 
430    Pwij(itypat)%qpt(:) = qpt_in(:)
431 
432    ! Prepare the call to paw_mkrhox
433    method = Paw_pwff(itypat)%method
434    dim1   = Paw_pwff(itypat)%dim1
435    dim2   = Paw_pwff(itypat)%dim2
436    nq_spl = Paw_pwff(itypat)%nq_spl
437    gmet   = Paw_pwff(itypat)%gmet
438 
439    ABI_STAT_MALLOC(Pwij(itypat)%mqpgij,(2,npw,lmn2_size), ierr)
440    ABI_CHECK(ierr==0, 'Out of memory in %mqpgij')
441 
442    ! Evaluate oscillator matrix elements mqpgij
443    call paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,Paw_pwff(itypat)%qgrid_spl,Paw_pwff(itypat)%pwff_spl,&
444 &    gmet,qpt_in,npw,gvec,ylm_q,Psps,Pawtab,Pwij(itypat)%mqpgij)
445  end do !itypat
446 
447  ABI_FREE(ylm_q)
448 
449 end subroutine pawpwij_init