TABLE OF CONTENTS
- ABINIT/m_pawpwij
- ABINIT/paw_cross_rho_tw_g
- m_pawpwij/paw_mkrhox
- m_pawpwij/paw_mkrhox_spl
- m_pawpwij/paw_rho_tw_g
- m_pawpwij/pawpwff_free
- m_pawpwij/pawpwff_init
- m_pawpwij/pawpwff_t
- m_pawpwij/pawpwij_free_d1
- m_pawpwij/pawpwij_free_d2
- m_pawpwij/pawpwij_t
- m_pawrhox/pawpwij_init
ABINIT/m_pawpwij [ 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-2024 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
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 MODULE m_pawpwij 26 27 use defs_basis 28 use m_errors 29 use m_abicore 30 use m_fft 31 32 use m_fstrings, only : sjoin, itoa 33 use defs_datatypes, only : pseudopotential_type 34 use defs_abitypes, only : MPI_type 35 use m_numeric_tools, only : arth 36 use m_geometry, only : metric 37 use m_crystal, only : crystal_t 38 use m_paw_numeric, only : paw_jbessel_4spline, paw_spline 39 use m_splines, only : splfit 40 use m_pawang, only : pawang_type 41 use m_paw_sphharm, only : realgaunt 42 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free, pawrad_copy, simp_gen 43 use m_pawtab, only : pawtab_type 44 use m_pawcprj, only : pawcprj_type 45 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq 46 use m_initylmg, only : initylmg 47 48 implicit none 49 50 private
ABINIT/paw_cross_rho_tw_g [ 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
SOURCE
1119 subroutine paw_cross_rho_tw_g(nspinor,npwvec,nr,ngfft,map2sphere,use_padfft,igfftg0,gbound,& 1120 & ur_ae1,ur_ae_onsite1,ur_ps_onsite1,i1,ktabr1,ktabp1,spinrot1,& 1121 & ur_ae2,ur_ae_onsite2,ur_ps_onsite2,i2,ktabr2,ktabp2,spinrot2,& 1122 & dim_rtwg,rhotwg) 1123 1124 !Arguments ------------------------------------ 1125 !scalars 1126 integer,intent(in) :: i1,i2,npwvec,nr,nspinor,dim_rtwg,map2sphere,use_padfft 1127 complex(dpc),intent(in) :: ktabp1,ktabp2 1128 !arrays 1129 integer,intent(in) :: gbound(:,:) 1130 integer,intent(in) :: igfftg0(npwvec*map2sphere),ngfft(18) 1131 integer,intent(in) :: ktabr1(nr),ktabr2(nr) 1132 real(dp),intent(in) :: spinrot1(4),spinrot2(4) 1133 complex(gwpc),intent(in) :: ur_ae1(nr),ur_ae2(nr) 1134 complex(gwpc),intent(in) :: ur_ae_onsite1(nr),ur_ae_onsite2(nr) 1135 complex(gwpc),intent(in) :: ur_ps_onsite1(nr),ur_ps_onsite2(nr) 1136 complex(gwpc),intent(inout) :: rhotwg(npwvec*dim_rtwg) 1137 1138 !Local variables------------------------------- 1139 !scalars 1140 integer,parameter :: ndat1 = 1, fftcache0 = 0, gpu_option_0 = 0 1141 integer :: ig,igfft,nx,ny,nz,ldx,ldy,ldz,mgfft,isprot1,isprot2 1142 type(fftbox_plan3_t) :: plan 1143 !arrays 1144 complex(dpc),allocatable :: usk(:),uu(:),rho(:) 1145 1146 ! ************************************************************************* 1147 1148 SELECT CASE (nspinor) 1149 1150 CASE (1) ! Collinear case. 1151 1152 ABI_MALLOC(uu,(nr)) 1153 ABI_MALLOC(usk,(nr)) 1154 ABI_MALLOC(rho,(nr)) 1155 1156 uu = ur_ae1(ktabr1)*ktabp1 - ur_ae_onsite1(ktabr1)*ktabp1; if (i1==1) uu = CONJG(uu) 1157 usk = ur_ae_onsite2(ktabr2)*ktabp2 - ur_ps_onsite2(ktabr2)*ktabp2; if (i2==2) usk = CONJG(usk) 1158 rho = uu * usk 1159 1160 uu = ur_ae_onsite1(ktabr1)*ktabp1 - ur_ps_onsite1(ktabr1)*ktabp1; if (i1==1) uu = CONJG(uu) 1161 usk = ur_ae2(ktabr2)*ktabp2 - ur_ae_onsite2(ktabr2)*ktabp2; if (i2==2) usk = CONJG(usk) 1162 rho = rho + uu * usk 1163 1164 SELECT CASE (map2sphere) 1165 1166 CASE (0) 1167 ! Need results on the full FFT box thus cannot use zero-padded FFT. 1168 call plan%init(ndat1, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0) 1169 call plan%execute(rho, -1) 1170 call plan%free() 1171 1172 rhotwg=rhotwg + rho 1173 1174 CASE (1) 1175 ! Need results on the G-sphere. Call zero-padded FFT routines if required. 1176 if (use_padfft==1) then 1177 nx =ngfft(1); ny =ngfft(2); nz =ngfft(3); mgfft = MAXVAL(ngfft(1:3)) 1178 ldx=nx ; ldy=ny ; ldz=nz 1179 call fftpad(rho,ngfft,nx,ny,nz,ldx,ldy,ldz,ndat1,mgfft,-1,gbound) 1180 else 1181 call plan%init(ndat1, ngfft(1:3), ngfft(1:3), ngfft(7), fftcache0, gpu_option_0) 1182 call plan%execute(rho, -1) 1183 call plan%free() 1184 end if 1185 1186 ! Have to map FFT to G-sphere. 1187 do ig=1,npwvec 1188 igfft=igfftg0(ig) 1189 ! If G-G0 belong to the FFT mesh. 1190 if (igfft/=0) rhotwg(ig)=rhotwg(ig)+rho(igfft) 1191 end do 1192 1193 CASE DEFAULT 1194 ABI_BUG(sjoin("Wrong map2sphere:", itoa(map2sphere))) 1195 END SELECT 1196 1197 RETURN 1198 1199 CASE (2) 1200 ! Spinorial case. 1201 isprot1=spinrot1(1); isprot2=spinrot2(1) ! This is to bypass abirule 1202 ABI_ERROR("Spinorial case not implemented yet") 1203 1204 SELECT CASE (map2sphere) 1205 1206 CASE (0) ! Need results on the full FFT box thus cannot use zero-padded FFT. 1207 CASE (1) ! Need results on the G-sphere. Call zero-padded FFT routines if required. 1208 CASE DEFAULT 1209 ABI_BUG("Wrong map2sphere") 1210 END SELECT 1211 RETURN 1212 1213 CASE DEFAULT 1214 ABI_BUG(sjoin('Wrong nspinor:', itoa(nspinor))) 1215 END SELECT 1216 1217 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)
SOURCE
794 subroutine paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,qgrid_spl,pwff_spl,& 795 gmet,qpt,npw,gvec,ylm_q,Psps,Pawtab,paw_rhox) 796 797 !Arguments ------------------------------------ 798 !scalars 799 integer,intent(in) :: itypat,dim1,dim2,method,npw,nq_spl,lmn2_size 800 type(Pseudopotential_type),intent(in) :: Psps 801 !arrays 802 integer,intent(in) :: gvec(3,npw) 803 real(dp),intent(in) :: gmet(3,3) 804 real(dp),intent(in) :: pwff_spl(nq_spl,2,0:dim1,dim2) 805 real(dp),intent(in) :: qpt(3),ylm_q(npw,(2*Psps%mpsang-1)**2) 806 real(dp),intent(out) :: paw_rhox(2,npw,lmn2_size) 807 real(dp),intent(in) :: qgrid_spl(nq_spl) 808 type(Pawtab_type),target,intent(in) :: Pawtab(Psps%ntypat) 809 810 !Local variables------------------------------- 811 !scalars 812 integer :: ider,ig,ignt,il,ilm,ilm_G,ilmn,iln,im,ipow,jl,jlm 813 integer :: jlmn,jln,jm,k0lm,k0lmn,k0ln,klm,klmn,kln,ll_G,mm_G,mpsang,ngnt 814 real(dp) :: rgnt,dummy 815 character(len=500) :: msg 816 !arrays 817 integer,allocatable :: gntselect(:,:) 818 integer,ABI_CONTIGUOUS pointer :: indlmn(:,:) 819 real(dp) :: mi_l(2,0:3),qpg(3) 820 real(dp),allocatable :: derfun(:),newfun(:),qpg_norm(:),realgnt(:),wk_ffnl(:,:) 821 822 ! ************************************************************************* 823 824 !write(std_out,*)itypat,dim1,dim2,method,npw,nq_spl,lmn2_size,Psps%mpsang 825 mpsang = Psps%mpsang 826 indlmn => Pawtab(itypat)%indlmn 827 828 ! === Pre-calculate (-i)^l === 829 mi_l(1,0)=one ; mi_l(2,0)=zero 830 mi_l(1,1)=zero ; mi_l(2,1)=-one 831 mi_l(1,2)=-one ; mi_l(2,2)=zero 832 mi_l(1,3)=zero ; mi_l(2,3)=one 833 834 ! === Calculate |q+G| === 835 ! * 2\pi is not included to be consistent with the spline. 836 ABI_MALLOC(qpg_norm,(npw)) 837 do ig=1,npw 838 qpg = qpt + gvec(:,ig) 839 qpg_norm(ig)=SQRT(DOT_PRODUCT(qpg,MATMUL(gmet,qpg))) 840 end do 841 842 ! Check q-grid as %qgrid_spl must be large enoung. 843 if (MAXVAL(qpg_norm)>MAXVAL(qgrid_spl)) then 844 write(msg,'(3a,f8.4,a,f8.4,2a)')& 845 ' Function values are being requested outside range of data. ',ch10,& 846 ' Max qpg_norm = ',MAXVAL(qpg_norm),' Max qgrid_spl = ',MAXVAL(qgrid_spl),ch10,& 847 ' Increase ecut(wfn), check qrid_ff and gsqcut ' 848 ABI_ERROR(msg) 849 end if 850 851 ABI_MALLOC(wk_ffnl,(nq_spl,2)) 852 ABI_MALLOC(newfun,(npw)) 853 ABI_MALLOC(derfun,(npw)) 854 855 SELECT CASE (method) 856 CASE (PWIJ_ARNAUD) 857 ! === Arnaud-Alouani exact expression === 858 ! * It does not describe the multipoles of the AE charge density 859 ! * $ 4\pi \sum_{LM} (-i)^l Y_M^L(q+G) G_{\li\mi\lj\mj}^{\LM} ff^{aL}_{ij}(|q+G|) $ 860 ! where f has been calculated in paw_mkrhox_spl 861 ! 862 ! === Re-evaluate Gaunt coefficients, just to be on the safe side === 863 ! * Note that gntselect is in packed form, thanks to invariance under permutation. 864 ! * Could use Pawang% but size of gntselect depends on pawxcdev! 865 866 ABI_MALLOC( realgnt,((2*mpsang-1)**2*(mpsang)**4)) 867 ABI_MALLOC(gntselect,((2*mpsang-1)**2, mpsang**2*(mpsang**2+1)/2)) 868 call realgaunt(mpsang,ngnt,gntselect,realgnt) 869 paw_rhox=zero 870 871 ! Loop on (jl,jm,jn) channels for this atom 872 do jlmn=1,Pawtab(itypat)%lmn_size 873 jl =indlmn(1,jlmn) 874 jm =indlmn(2,jlmn) 875 jlm=indlmn(4,jlmn) 876 jln=indlmn(5,jlmn) 877 878 k0lmn=jlmn*(jlmn-1)/2 879 k0lm =jlm *(jlm -1)/2 880 k0ln =jln *(jln -1)/2 881 882 ! === Loop on (il,im,in) channels; klmn is index for packed form === 883 do ilmn=1,jlmn 884 il =indlmn(1,ilmn) 885 im =indlmn(2,ilmn) 886 ilm=indlmn(4,ilmn) 887 iln=indlmn(5,ilmn) 888 889 klmn=k0lmn+ilmn 890 klm =k0lm +ilm 891 kln =k0ln +iln 892 893 ! === Summing over allowed (L,M), taking into account Gaunt selection rules === 894 do ll_G=ABS(jl-il),jl+il,2 895 ipow=MOD(ll_G,4) 896 ider=0 897 wk_ffnl(:,:)=pwff_spl(:,:,ll_G,kln) 898 call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw) 899 do mm_G=-ll_G,ll_G 900 ilm_G=1+ll_G**2+ll_G+mm_G 901 ignt=gntselect(ilm_G,klm) 902 if (ignt==0) CYCLE 903 rgnt=realgnt(ignt) 904 905 ! === Evaluate matrix elements for each plane wave === 906 do ig=1,npw 907 dummy = newfun(ig) * ylm_q(ig,ilm_G) * rgnt 908 paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) + ( dummy * mi_l(1,ipow) ) 909 paw_rhox(2,ig,klmn) = paw_rhox(2,ig,klmn) + ( dummy * mi_l(2,ipow) ) 910 end do 911 end do !mm_G 912 end do !ll_G 913 914 end do !ilmn 915 end do !jlmn 916 917 ! Multiply by 4\pi arising from the expansion of the plane wave 918 paw_rhox = four_pi*paw_rhox 919 ABI_FREE(realgnt) 920 ABI_FREE(gntselect) 921 922 CASE (PWIJ_SHISHKIN) 923 ! === Shishkin-Kresse approximated expression ==== 924 ! * Better description of multipoles of AE charge, 925 ! * Better results for energy degeneracies in GW band structure 926 ! * $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 927 ! 928 paw_rhox=zero 929 ! 930 ! === Loop on (jl,jm,jn) channels for this atom === 931 !itypat=Cryst%typat(iatm) 932 do jlmn=1,Pawtab(itypat)%lmn_size 933 jl =indlmn(1,jlmn) 934 jm =indlmn(2,jlmn) 935 jlm=indlmn(4,jlmn) 936 jln=indlmn(5,jlmn) 937 938 k0lmn=jlmn*(jlmn-1)/2 939 k0lm =jlm *(jlm -1)/2 940 k0ln =jln *(jln -1)/2 941 ! 942 ! === Loop on (il,im,in) channels; klmn is index for packed form === 943 do ilmn=1,jlmn 944 il =indlmn(1,ilmn) 945 im =indlmn(2,ilmn) 946 ilm=indlmn(4,ilmn) 947 iln=indlmn(5,ilmn) 948 949 klmn=k0lmn+ilmn 950 klm =k0lm +ilm 951 kln =k0ln +iln 952 ! 953 ! === Summing over allowed (l,m), taking into account Gaunt selection rules === 954 do ll_G=ABS(jl-il),jl+il,2 955 ipow=MOD(ll_G,4) 956 do mm_G=-ll_G,ll_G 957 ! here I can move splfit before the loop over mm_G but I have to change paw_rhox_spl 958 ilm_G=1+ll_G**2+ll_G+mm_G 959 ider=0 960 wk_ffnl(:,:)=pwff_spl(:,:,ilm_G-1,klmn) ! Note klmn and ilm_G-1 961 call splfit(qgrid_spl,derfun,wk_ffnl,ider,qpg_norm,newfun,nq_spl,npw) 962 ! 963 ! === Evaluate matrix elements for each plane wave === 964 do ig=1,npw 965 dummy = newfun(ig) * ylm_q(ig,ilm_G) 966 paw_rhox(1,ig,klmn) = paw_rhox(1,ig,klmn) & 967 + dummy * mi_l(1,ipow) !(ph3d(1,ig)*mi_l(1,ipow)-ph3d(2,ig)*mi_l(2,ipow)) 968 paw_rhox(2,ig,klmn) = paw_rhox(2,ig,klmn) & 969 + dummy * mi_l(2,ipow) !(ph3d(1,ig)*mi_l(2,ipow)+ph3d(2,ig)*mi_l(1,ipow)) 970 end do 971 end do !mm_G 972 end do !ll_G 973 974 end do !ilmn 975 end do !jlmn 976 977 ! Multiply by 4\pi arising from the expansion of the plane wave 978 paw_rhox=four_pi*paw_rhox 979 980 CASE DEFAULT 981 ABI_BUG(sjoin('Wrong value for method:', itoa(method))) 982 END SELECT 983 984 ABI_FREE(wk_ffnl) 985 ABI_FREE(newfun) 986 ABI_FREE(derfun) 987 ABI_FREE(qpg_norm) 988 989 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!
SOURCE
493 subroutine paw_mkrhox_spl(itypat,ntypat,method,dim1,dim2,nq_spl,qgrid_spl,Pawrad,Pawtab,pwff_spl) 494 495 !Arguments ------------------------------------ 496 !scalars 497 integer,intent(in) :: itypat,ntypat,method,dim2,dim1,nq_spl 498 !arrays 499 real(dp),intent(in) :: qgrid_spl(nq_spl) 500 real(dp),intent(out) :: pwff_spl(nq_spl,2,0:dim1,dim2) 501 type(Pawrad_type),intent(in) :: Pawrad(ntypat) 502 type(Pawtab_type),intent(in) :: Pawtab(ntypat) 503 504 !Local variables------------------------------- 505 !scalars 506 integer :: mm,nlmn,jlmn,ilmn,klmn,ij_size,l_size !,ider 507 integer :: iq,ir,ll,meshsz,mmax,iln,jln,nln,k0ln,kln,qlm 508 real(dp),parameter :: EPS=tol14**4,TOLJ=0.001_dp 509 real(dp) :: arg,argn,bes,besp,qr,yp1,ypn 510 character(len=500) :: msg 511 type(Pawrad_type) :: Tmpmesh 512 !arrays 513 real(dp),allocatable :: rrshape_l(:),shape_l(:),ff(:),gg(:),rr(:),rrdphi_ij(:) 514 real(dp),allocatable :: dphi_ij(:),tmp_spl(:,:,:,:),tmp_jgl(:,:,:) 515 516 !************************************************************************* 517 518 DBG_ENTER("COLL") 519 520 pwff_spl=zero 521 522 SELECT CASE (method) 523 524 CASE (PWIJ_ARNAUD) 525 ! === Arnaud-Alouani exact expression PRB 62. 4464 [[cite:Arnaud2000]] === 526 ! * $ff_^{aL}_{ij}(q) = 527 ! \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$ 528 ! * It does not descrive correctly the multipoles of the AE charge density if low cutoff on G 529 write(msg,'(a,i3)')' paw_mkrhox_spl: Using Arnaud-Alouani expression for atom type: ',itypat 530 call wrtout(std_out,msg) 531 532 nln = Pawtab(itypat)%basis_size 533 ij_size = Pawtab(itypat)%ij_size 534 l_size = Pawtab(itypat)%l_size 535 536 ABI_MALLOC(tmp_spl,(nq_spl,2,0:l_size-1,ij_size)) 537 538 ! Is mesh beginning with r=0 ? 539 if (abs(Pawrad(itypat)%rad(1)) > tol10) then 540 ABI_ERROR("Radial mesh starts with r/=0") 541 end if 542 ! 543 ! === Initialize temporary arrays and variables === 544 meshsz = Pawtab(itypat)%mesh_size ; mmax=meshsz 545 546 ABI_MALLOC(dphi_ij,(meshsz)) 547 ABI_MALLOC(rrdphi_ij,(meshsz)) 548 ABI_MALLOC(ff,(meshsz)) 549 ABI_CALLOC(gg,(meshsz)) 550 ABI_MALLOC(rr,(meshsz)) 551 rr(1:meshsz) = Pawrad(itypat)%rad(1:meshsz) 552 ! 553 ! === Loop on (jln,iln) channels for this type. Packed form === 554 tmp_spl=zero 555 556 do jln=1,nln 557 k0ln=jln*(jln-1)/2 558 do iln=1,jln 559 kln=k0ln+iln 560 561 dphi_ij(1:meshsz) = Pawtab(itypat)%phiphj(1:meshsz,kln)-Pawtab(itypat)%tphitphj(1:meshsz,kln) 562 rrdphi_ij(1:meshsz) = rr(1:meshsz)*dphi_ij(1:meshsz) ! will be used for first derivative 563 564 ir=meshsz 565 do while (ABS(dphi_ij(ir))<EPS) 566 ir=ir-1 567 end do 568 mmax=MIN(ir+1,meshsz) 569 ! ir is equal to meshsz if no point are below EPS 570 if (mmax/=Pawrad(itypat)%int_meshsz) then ! mmax=meshsz 571 call pawrad_init(Tmpmesh,mesh_size=Pawtab(itypat)%mesh_size,mesh_type=Pawrad(itypat)%mesh_type, & 572 & rstep=Pawrad(itypat)%rstep,lstep=Pawrad(itypat)%lstep,r_for_intg=rr(mmax)) 573 else 574 call pawrad_copy(Pawrad(itypat),Tmpmesh) 575 end if 576 ! 577 ! === Loop on l for Bessel function. Note the starting point === 578 ! TODO Here I should loop only the moments allowed by Gaunt, I should use indklm! 579 ! and only on lmax for this atom 580 do ll=0,l_size-1 581 ! 582 ! === Compute f_l(q=0) only if l=0, and first derivative fp_l(q=0) (nonzero only if ll==1) === 583 tmp_spl(1,1,ll,kln)=zero; yp1=zero 584 if (ll==0) then 585 call simp_gen(tmp_spl(1,1,ll,kln),dphi_ij,Tmpmesh) 586 end if 587 if (ll==1) then 588 call simp_gen(yp1,rrdphi_ij,Tmpmesh) 589 yp1=yp1*two_pi*third 590 end if 591 ! 592 ! === Compute f_l(0<q<qmax) === 593 if (nq_spl>2) then 594 do iq=2,nq_spl-1 595 arg=two_pi*qgrid_spl(iq) 596 do ir=1,mmax 597 qr=arg*rr(ir) 598 call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ) 599 ff(ir)=bes*dphi_ij(ir) 600 end do 601 call simp_gen(tmp_spl(iq,1,ll,kln),ff,Tmpmesh) 602 end do 603 end if 604 ! 605 ! === Compute f_l(q=qmax) and first derivative === 606 if (nq_spl>1) then 607 argn=two_pi*qgrid_spl(nq_spl) 608 do ir=1,mmax 609 qr=argn*rr(ir) 610 call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ) 611 ff(ir)=bes * dphi_ij(ir) 612 gg(ir)=besp*rrdphi_ij(ir) 613 end do 614 call simp_gen(tmp_spl(nq_spl,1,ll,kln),ff,Tmpmesh) 615 gg(:)=two_pi*gg(:) ! two_pi comes from 2\pi|q| in the Bessel function 616 call simp_gen(ypn,gg,Tmpmesh) 617 else 618 ypn=yp1 619 end if 620 ! 621 ! === Compute second derivative of ff^{al}_{ij)(q) === 622 !yp1=zero; ypn=zero 623 call paw_spline(qgrid_spl,tmp_spl(:,1,ll,kln),nq_spl,yp1,ypn,tmp_spl(:,2,ll,kln)) 624 end do !ll 625 626 call pawrad_free(Tmpmesh) 627 628 end do !iln 629 end do !jln 630 ! 631 ! === Save values for this atom type, each ll and kln channel === 632 pwff_spl = tmp_spl 633 634 ABI_FREE(dphi_ij) 635 ABI_FREE(rrdphi_ij) 636 ABI_FREE(ff) 637 ABI_FREE(gg) 638 ABI_FREE(rr) 639 ABI_FREE(tmp_spl) 640 641 !if (.FALSE.) then ! write form factors on file for plotting purpose. 642 ! ll=0 643 ! do iq=1,nq_spl 644 ! write(777+itypat,'(50(es16.8))')qgrid_spl(iq),((pwff_spl(iq,ider,ll,iln),ider=1,2),iln=1,dim2) 645 ! end do 646 !end if 647 648 CASE (PWIJ_SHISHKIN) 649 ! ==== Shishkin-Kresse approximated expression ==== 650 ! $ff_^{aL}_{ij}(q) = q_ij^{LM} \int_0^{r_a} j_L(2\pi qr) g_L(r) r^2 dr$ 651 ! * Better description of multipoles of AE charge, 652 ! * Better results for energy degeneracies in GW band structure 653 write(msg,'(a,i3)')' paw_mkrhox_spl: Using Shishkin-Kresse expression for atom type ',itypat 654 call wrtout(std_out, msg) 655 l_size = Pawtab(itypat)%l_size 656 nlmn = Pawtab(itypat)%lmn_size 657 658 !allocate(tmp_jgl(nq_spl,2,0:2*(Psps%mpsang-1))) 659 ABI_MALLOC(tmp_jgl,(nq_spl,2,0:l_size-1)) 660 661 ! Is mesh beginning with r=0 ? 662 if (ABS(Pawrad(itypat)%rad(1))>tol10) then 663 ABI_ERROR("Radial mesh starts with r/=0") 664 end if 665 666 ! === Initialize temporary arrays and variables === 667 call pawrad_copy(Pawrad(itypat),Tmpmesh) 668 meshsz=Pawtab(itypat)%mesh_size ; mmax=meshsz 669 ABI_MALLOC(ff,(meshsz)) 670 ABI_CALLOC(gg,(meshsz)) 671 ABI_MALLOC(rr,(meshsz)) 672 ABI_MALLOC(shape_l,(meshsz)) 673 ABI_MALLOC(rrshape_l,(meshsz)) 674 rr(1:meshsz)=Tmpmesh%rad(1:meshsz) 675 676 tmp_jgl(:,:,:)=zero 677 rrshape_l(1)=zero 678 shape_l(1)=zero 679 ! 680 ! TODO Here I should loop only the moments allowed by Gaunt, I should use indklm! 681 ! and only lmax for this atom 682 !do ll=0,2*(Psps%mpsang-1) 683 do ll=0,l_size-1 684 shape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**2 685 rrshape_l(2:meshsz)=Pawtab(itypat)%shapefunc(2:meshsz,ll+1)*rr(2:meshsz)**3 686 ! 687 ! === Compute f_l(q=0) and first derivative fp_l(q=0) (only if ll==1) === 688 tmp_jgl(1,1,ll)=zero ; yp1=zero 689 if (ll==0) then 690 call simp_gen(tmp_jgl(1,1,ll),shape_l,Tmpmesh) 691 end if 692 if (ll==1) then 693 call simp_gen(yp1,rrshape_l,Tmpmesh) !rr comes from d/dq 694 yp1=yp1*two_pi*third 695 end if 696 ! 697 ! === Compute f_l(0<q<qmax) === 698 if (nq_spl>2) then 699 do iq=2,nq_spl-1 700 arg=two_pi*qgrid_spl(iq) 701 do ir=1,mmax 702 qr=arg*rr(ir) 703 call paw_jbessel_4spline(bes,besp,ll,0,qr,TOLJ) 704 ff(ir)=bes*shape_l(ir) 705 end do 706 call simp_gen(tmp_jgl(iq,1,ll),ff,Tmpmesh) 707 end do 708 end if 709 ! 710 ! === Compute f_l(q=qmax) and first derivative === 711 if (nq_spl>1) then 712 argn=two_pi*qgrid_spl(nq_spl) 713 do ir=1,mmax 714 qr=argn*rr(ir) 715 call paw_jbessel_4spline(bes,besp,ll,1,qr,TOLJ) 716 ff(ir)=bes * shape_l(ir) 717 gg(ir)=besp*rrshape_l(ir) 718 end do 719 call simp_gen(tmp_jgl(nq_spl,1,ll),ff,Tmpmesh) 720 gg(:)=two_pi*gg(:) !two_pi comes from 2\pi|q| 721 call simp_gen(ypn,gg,Tmpmesh) 722 else 723 ypn=yp1 724 end if 725 ! 726 ! === Compute second derivative of ff_^{al}_{ij)(q) === 727 call paw_spline(qgrid_spl,tmp_jgl(:,1,ll),nq_spl,yp1,ypn,tmp_jgl(:,2,ll)) 728 end do !ll 729 ! 730 ! === Save values for this type, each ll and ilmn,jlm channels === 731 ! * Here we assembly q_{ij}^{lm} \int_0^{r_a} j_l(2\pi(q+G)r) g_l(r)r^2 dr 732 ! * Some of the contributions from qijl are zero due to Gaunt selection rules 733 do jlmn=1,nlmn 734 do ilmn=1,jlmn 735 klmn=ilmn+(jlmn-1)*jlmn/2 736 !do ll=0,2*(Psps%mpsang-1) 737 do ll=0,l_size-1 738 do mm=-ll,ll 739 qlm=1+ll**2+ll+mm 740 pwff_spl(:,:,qlm-1,klmn)=tmp_jgl(:,:,ll)*Pawtab(itypat)%qijl(qlm,klmn) 741 end do 742 end do 743 end do 744 end do 745 746 call pawrad_free(Tmpmesh) 747 ABI_FREE(shape_l) 748 ABI_FREE(rrshape_l) 749 ABI_FREE(ff) 750 ABI_FREE(gg) 751 ABI_FREE(rr) 752 ABI_FREE(tmp_jgl) 753 754 CASE DEFAULT 755 ABI_BUG(sjoin('Called with wrong value for method:', itoa(method))) 756 END SELECT 757 758 DBG_EXIT("COLL") 759 760 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
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.
SOURCE
1022 pure subroutine paw_rho_tw_g(cryst, pwij, npw, dim_rtwg, nspinor, gvec, Cprj_kmqb1, Cprj_kb2, rhotwg) 1023 1024 !Arguments ------------------------------------ 1025 !scalars 1026 type(crystal_t),intent(in) :: cryst 1027 integer,intent(in) :: npw,nspinor,dim_rtwg 1028 !arrays 1029 integer,intent(in) :: gvec(3,npw) 1030 type(pawcprj_type),intent(in) :: Cprj_kmqb1(cryst%natom,nspinor), Cprj_kb2(cryst%natom,nspinor) 1031 type(pawpwij_t),intent(in) :: Pwij(cryst%ntypat) 1032 complex(gwpc),intent(inout) :: rhotwg(npw*dim_rtwg) 1033 1034 !Local variables------------------------------- 1035 !scalars 1036 integer :: ig,iat,nlmn,ilmn,jlmn,k0lmn,klmn,iab,isp1,isp2,spad,itypat 1037 real(dp) :: fij,re_psp,im_psp,re_pw,im_pw,arg 1038 !arrays 1039 integer,parameter :: spinor_idxs(2,4) = RESHAPE([1,1,2,2,1,2,2,1], [2, 4]) 1040 real(dp) :: tmp(2),qpg(3),x0(3),ph3d(2) 1041 1042 ! ************************************************************************* 1043 1044 ! Loop over the four spinorial combinations 1045 do iab=1,dim_rtwg 1046 isp1 = spinor_idxs(1,iab) 1047 isp2 = spinor_idxs(2,iab) 1048 spad = npw*(iab-1) 1049 1050 do ig=1,npw 1051 tmp(:)=zero 1052 do iat=1,cryst%natom 1053 itypat = cryst%typat(iat) 1054 nlmn = Pwij(itypat)%lmn_size 1055 x0(:) = cryst%xred(:,iat) 1056 1057 ! Structure factor e^{-i(q+G)*xred} 1058 qpg(:)= Pwij(itypat)%qpt(:) + gvec(:,ig) 1059 arg=-two_pi*DOT_PRODUCT(qpg(:),x0) 1060 ph3d(1)=COS(arg) 1061 ph3d(2)=SIN(arg) 1062 1063 ! Loop over [(jl,jm,jn), (il,im,in)] channels in packed storage mode. 1064 do jlmn=1,nlmn 1065 k0lmn=jlmn*(jlmn-1)/2 1066 do ilmn=1,jlmn 1067 re_psp = Cprj_kmqb1(iat,isp1)%cp(1,ilmn) * Cprj_kb2(iat,isp2)%cp(1,jlmn) & 1068 +Cprj_kmqb1(iat,isp1)%cp(2,ilmn) * Cprj_kb2(iat,isp2)%cp(2,jlmn) & 1069 +Cprj_kmqb1(iat,isp1)%cp(1,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn) & 1070 +Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(2,ilmn) 1071 1072 im_psp = Cprj_kmqb1(iat,isp1)%cp(1,ilmn) * Cprj_kb2(iat,isp2)%cp(2,jlmn) & 1073 -Cprj_kmqb1(iat,isp1)%cp(2,ilmn) * Cprj_kb2(iat,isp2)%cp(1,jlmn) & 1074 +Cprj_kmqb1(iat,isp1)%cp(1,jlmn) * Cprj_kb2(iat,isp2)%cp(2,ilmn) & 1075 -Cprj_kmqb1(iat,isp1)%cp(2,jlmn) * Cprj_kb2(iat,isp2)%cp(1,ilmn) 1076 1077 klmn=k0lmn+ilmn; fij=one; if (jlmn==ilmn) fij=half 1078 1079 ! Multiply by the phase due to the atom position. 1080 re_pw = Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(1) & 1081 -Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(2) 1082 1083 im_pw = Pwij(itypat)%mqpgij(1,ig,klmn) * ph3d(2) & 1084 +Pwij(itypat)%mqpgij(2,ig,klmn) * ph3d(1) 1085 1086 tmp(1)=tmp(1)+ fij * (re_pw*re_psp - im_pw*im_psp) 1087 tmp(2)=tmp(2)+ fij * (re_pw*im_psp + im_pw*re_psp) 1088 end do !ilmn 1089 end do !jlmn 1090 end do !iat 1091 1092 !if(ig==1) write(std_out,*) " TOTAL PW osc str = ",rhotwg(ig+spad) 1093 !if(ig==1) write(std_out,*) " TOTAL PAW osc str = ",tmp(1),tmp(2) 1094 ! Update input data using the appropriate index. 1095 rhotwg(ig+spad) = rhotwg(ig+spad) + CMPLX(tmp(1),tmp(2),kind=gwpc) 1096 !if(ig==1) write(std_out,*) " TOTAL PW+PAW osc str = ",rhotwg(ig+spad) 1097 end do !ig 1098 1099 end do ! dim_rtwg 1100 1101 end subroutine paw_rho_tw_g
m_pawpwij/pawpwff_free [ Functions ]
[ Top ] [ m_pawpwij ] [ Functions ]
NAME
pawpwff_free
FUNCTION
Free memory
SOURCE
246 subroutine pawpwff_free(Paw_pwff) 247 248 !Arguments ------------------------------------ 249 !scalars 250 type(pawpwff_t),intent(inout) :: Paw_pwff(:) 251 252 !Local variables------------------------------- 253 integer :: ii 254 255 !************************************************************************ 256 257 do ii=1,SIZE(Paw_pwff) 258 ABI_SFREE(Paw_pwff(ii)%qgrid_spl) 259 ABI_SFREE(Paw_pwff(ii)%pwff_spl) 260 end do 261 262 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.
SOURCE
174 subroutine pawpwff_init(Paw_pwff,method,nq_spl,qmax,gmet,Pawrad,Pawtab,Psps) 175 176 !Arguments ------------------------------------ 177 !scalars 178 integer,intent(in) :: method 179 type(Pseudopotential_type),intent(in) :: Psps 180 !arrays 181 integer,intent(in) :: nq_spl(Psps%ntypat) 182 real(dp),intent(in) :: gmet(3,3) 183 real(dp),intent(in) :: qmax(Psps%ntypat) 184 type(Pawrad_type),intent(in) :: Pawrad(Psps%ntypat) 185 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 186 type(pawpwff_t),intent(out) :: Paw_pwff(Psps%ntypat) 187 188 !Local variables------------------------------- 189 !scalars 190 integer :: dim1,dim2,itypat,nq 191 real(dp) :: dq 192 193 !************************************************************************ 194 195 !@pawpwff_t 196 197 ! === Evaluate form factors for the radial part of phi.phj-tphi.tphj === 198 do itypat=1,Psps%ntypat 199 Paw_pwff(itypat)%method = method 200 201 select case (method) 202 case (PWIJ_ARNAUD) 203 dim1 = Pawtab(itypat)%l_size-1 204 dim2 = Pawtab(itypat)%ij_size 205 case (PWIJ_SHISHKIN) 206 dim1 = Pawtab(itypat)%l_size**2 207 dim2 = Pawtab(itypat)%lmn2_size 208 case default 209 ABI_BUG(sjoin("Wrong method:", itoa(method))) 210 end select 211 212 Paw_pwff(itypat)%dim1 = dim1 213 Paw_pwff(itypat)%dim2 = dim2 214 Paw_pwff(itypat)%gmet = gmet 215 216 ! Setup of the q-mesh for spline. It can be type-dependent. 217 nq = nq_spl(itypat) 218 dq = qmax(itypat)/(one*(nq-1)) 219 !write(std_out,*)"nq,dq",nq,dq 220 221 Paw_pwff(itypat)%nq_spl = nq 222 ABI_MALLOC(Paw_pwff(itypat)%qgrid_spl,(nq)) 223 Paw_pwff(itypat)%qgrid_spl = arth(zero,dq,nq) 224 ! 225 ! === Calculate form factors depending on method === 226 ABI_MALLOC(Paw_pwff(itypat)%pwff_spl,(nq,2,0:dim1,dim2)) 227 228 call paw_mkrhox_spl(itypat,Psps%ntypat,method,dim1,dim2,nq, & 229 Paw_pwff(itypat)%qgrid_spl,Pawrad,Pawtab,Paw_pwff(itypat)%pwff_spl) 230 end do ! itypat 231 232 end subroutine pawpwff_init
m_pawpwij/pawpwff_t [ Types ]
[ Top ] [ m_pawpwij ] [ Types ]
NAME
pawpwff_t
FUNCTION
PAW form factors used to evaluate $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$
SOURCE
62 type,public :: pawpwff_t 63 64 integer :: method = -1 65 ! 1 For Arnaud-Alouani"s exact expression. 66 ! 2 For Shishkin-Kresse"s approximated expression. 67 68 integer :: dim1 = -1, dim2 = -1 69 ! Dimensions of pwff_spl, depending on method. 70 71 integer :: nq_spl = -1 72 ! Number of points in the reciprocal space grid on which 73 ! the radial integrals are evaluated. 74 75 real(dp) :: gmet(3,3) = zero 76 ! Reciprocal space metric tensor in Bohr**-2 77 78 real(dp),allocatable :: qgrid_spl(:) 79 ! qgrid_spl(nq_spl) 80 ! The coordinates of the points of the radial grid for the integrals used in the spline. 81 82 real(dp),allocatable :: pwff_spl(:,:,:,:) 83 ! pwff_spl(nq_spl,2,0:dim1,dim2) 84 ! The different integrals on the radial |q| grid, for a given atom type. 85 86 end type pawpwff_t 87 88 public :: pawpwff_init ! Initialize form factors for spline. 89 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
SOURCE
400 subroutine pawpwij_free_d1(Pwij) 401 402 !Arguments ------------------------------------ 403 type(pawpwij_t),intent(inout) :: Pwij(:) 404 405 !Local variables------------------------------- 406 integer :: ii 407 !************************************************************************ 408 409 do ii=1,SIZE(Pwij) 410 ABI_SFREE(Pwij(ii)%mqpgij) 411 end do 412 413 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
SOURCE
427 subroutine pawpwij_free_d2(Pwij) 428 429 !Arguments ------------------------------------ 430 !scalars 431 type(pawpwij_t),intent(inout) :: Pwij(:,:) 432 433 !Local variables------------------------------- 434 integer :: jj 435 436 !************************************************************************ 437 438 do jj=1,SIZE(Pwij,DIM=2) 439 call pawpwij_free_d1(Pwij(:,jj)) 440 end do 441 442 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
105 type,public :: pawpwij_t 106 107 integer :: istpw 108 ! Storage mode (similar to istwfk), not used at present 109 110 integer :: npw 111 ! The number of plane waves 112 113 integer :: lmn_size 114 ! Number of (l,m,n) elements for the paw basis 115 116 integer :: lmn2_size 117 ! lmn2_size=lmn_size*(lmn_size+1)/2 118 ! where lmn_size is the number of (l,m,n) elements for the paw basis 119 120 real(dp) :: qpt(3) 121 ! The q-point in e^{-i(q+G)}.r} 122 123 real(dp),allocatable :: mqpgij(:,:,:) 124 ! pwij(2,npw,lmn2_size) 125 ! $<phi|e^{-i(q+G).r}|phj> - <tphi|e^{-i(q+G).r}|tphj>$ 126 127 end type pawpwij_t 128 129 public :: pawpwij_init ! Calculate onsite matrix elements of a set of plane waves. 130 public :: pawpwij_free ! Deallocate dynamic memory in the structure. 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 ]
NAME
pawpwij_init
FUNCTION
Calculate 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 in gvec 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 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
OUTPUT
Pwij(%ntypat)<pawpwij_t>=Structure containing the onsite matrix elements of e^{-i(q+G).r}. Completely initialized in output.
SOURCE
295 subroutine pawpwij_init(Pwij, npw, qpt_in, gvec, rprimd, Psps, Pawtab, Paw_pwff) 296 297 !Arguments ------------------------------------ 298 !scalars 299 integer,intent(in) :: npw 300 type(Pseudopotential_type),intent(in) :: Psps 301 !arrays 302 integer,intent(in) :: gvec(3,npw) 303 real(dp),intent(in) :: qpt_in(3),rprimd(3,3) 304 type(Pawtab_type),intent(in) :: Pawtab(Psps%ntypat) 305 type(pawpwff_t),intent(in) :: Paw_pwff(Psps%ntypat) 306 type(pawpwij_t),intent(out) :: Pwij(Psps%ntypat) 307 308 !Local variables------------------------------- 309 !scalars 310 integer,parameter :: unkg0=0,unylm0=0 311 integer :: dim1,dim2,method, my_mqmem,my_nqpt,optder,two_lmaxp1,itypat 312 integer :: dummy_nsppol,lmn_size,lmn2_size,nq_spl,ierr 313 real(dp) :: ucvol 314 type(MPI_type) :: MPI_enreg_seq 315 !arrays 316 integer,allocatable :: npwarr(:),dummy_nband(:) 317 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 318 real(dp),allocatable :: my_qtmp(:,:), ylm_q(:,:),ylmgr_q(:,:,:) 319 320 ! ********************************************************************* 321 322 ! =============================================== 323 ! === Get real spherical harmonics in G space === 324 ! =============================================== 325 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 326 327 ! Set up of REAL Ylm(q+G) up to 2*l_max for this q-point. 328 my_mqmem=1; two_lmaxp1=2*Psps%mpsang-1; optder=0 329 330 ABI_MALLOC(ylm_q ,(npw*my_mqmem,two_lmaxp1**2)) 331 ABI_MALLOC(ylmgr_q,(npw*my_mqmem,3+6*(optder/2),two_lmaxp1**2)) 332 333 my_nqpt=1 334 ABI_MALLOC(my_qtmp,(3,my_nqpt)) 335 my_qtmp(:,1)=qpt_in(:) 336 337 ! dummy_nband and dummy_nsppol are not used in sequential mode. 338 dummy_nsppol=1 339 ABI_MALLOC(dummy_nband,(my_nqpt*dummy_nsppol)) 340 dummy_nband=0 341 ABI_MALLOC(npwarr,(my_nqpt)) 342 npwarr(:)=npw 343 344 ! Fake MPI_type for sequential part. 345 call initmpi_seq(MPI_enreg_seq) 346 347 call initylmg(gprimd,gvec,my_qtmp,my_mqmem,MPI_enreg_seq,two_lmaxp1,npw,dummy_nband,my_nqpt,& 348 npwarr,dummy_nsppol,optder,rprimd,ylm_q,ylmgr_q) 349 350 call destroy_mpi_enreg(MPI_enreg_seq) 351 352 ABI_FREE(my_qtmp) 353 ABI_FREE(dummy_nband) 354 ABI_FREE(npwarr) 355 ABI_FREE(ylmgr_q) 356 357 ! Construct the Pwij structure 358 do itypat=1,Psps%ntypat 359 360 Pwij(itypat)%istpw = 1 361 Pwij(itypat)%npw = npw 362 363 lmn_size = Pawtab(itypat)%lmn_size 364 lmn2_size = lmn_size*(lmn_size+1)/2 365 Pwij(itypat)%lmn_size = lmn_size 366 Pwij(itypat)%lmn2_size = lmn2_size 367 368 Pwij(itypat)%qpt(:) = qpt_in(:) 369 370 ! Prepare the call to paw_mkrhox 371 method = Paw_pwff(itypat)%method 372 dim1 = Paw_pwff(itypat)%dim1 373 dim2 = Paw_pwff(itypat)%dim2 374 nq_spl = Paw_pwff(itypat)%nq_spl 375 gmet = Paw_pwff(itypat)%gmet 376 377 ABI_MALLOC_OR_DIE(Pwij(itypat)%mqpgij,(2,npw,lmn2_size), ierr) 378 379 ! Evaluate oscillator matrix elements mqpgij 380 call paw_mkrhox(itypat,lmn2_size,method,dim1,dim2,nq_spl,Paw_pwff(itypat)%qgrid_spl,Paw_pwff(itypat)%pwff_spl,& 381 gmet,qpt_in,npw,gvec,ylm_q,Psps,Pawtab,Pwij(itypat)%mqpgij) 382 end do ! itypat 383 384 ABI_FREE(ylm_q) 385 386 end subroutine pawpwij_init