TABLE OF CONTENTS
ABINIT/pawsushat [ Functions ]
NAME
pawsushat
FUNCTION
PAW only, for susceptibility matrix: Compute contribution to the product of two wavefunctions (exchange charge density) from hat (compensation charge) density (in reciprocal space and eventually in real space): sushat_{ij,R}(g)=Sum_{L}[Q^L_ijR(g)]
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
atindx(natom)=index table for atoms, inverse of atindx
cprj_k(natom,nspinor*nband_k)= wave functions projected with non-local projectors:
cprj_k=<p_i|Cnk> where p_i is a non-local projector.
WARNING: cprj(iatom,:) ARE SORTED BY ATOM TYPE !!!
distribfft<type(distribfft_type)>=--optional-- contains infos related to FFT parallelism
gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for small FFT sphere.
gylmg_diel(npwdiel,lmax_diel**2,ntypat)= -PAW only- Fourier transform of g_l(r).Y_ml(r) shape functions
iband1,iband2= indices of the bands concerned with
ispinor1,ispinor2= indices of spinorial components concerned with
istwf_k=input option parameter that describes the storage of wfs
kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix.
lmax_diel=1+max. value of l angular momentum used for dielectric matrix
me_g0=--optional-- 1 if the current process treat the g=0 plane-wave (only needed when comm_fft is present)
mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix
mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
comm_atom=--optional-- MPI communicator over atoms
comm_fft=--optional-- MPI communicator over FT components
natom=number of atoms in cell
nband=number of bands at this k point for that spin polarization
ndiel4,ndiel5,ndiel6= FFT dimensions, modified to avoid cache trashing
nfftdiel=number of FFT grid points for the small (diel) grid
ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix
nspinor=number of spinorial components of the wavefunctions
ntypat=number of types of atoms in unit cell.
optreal=0 if WF product has to be output in reciprocal space
1 if WF product has to be output in real space
paral_kgb=--optional-- 1 if "band-FFT" parallelism is activated (only needed when comm_fft is present)
pawang <type(pawang_type)>=paw angular mesh and related data
pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
ph3d_diel(2,npwdiel,natom*usepaw)=3-dim structure factors, for each atom and plane wave, for dielectric matrix
typat(natom)=type (integer) for each atom
SIDE EFFECTS
=== if optreal=0
wfprod(2,npwdiel)=PAW contrib. to product of two wavefunctions (iband1,iband2):
is added (in reciprocal space)
=== if optreal=1
wfraug(2,ndiel4,ndiel5,ndiel6)=PAW contrib. to product of two wavefunctions (iband1,iband2)
is added (in real space)
PARENTS
susk,suskmm
CHILDREN
destroy_distribfft,fourwf,free_my_atmtab,get_my_atmtab init_distribfft_seq,initmpi_seq,set_mpi_enreg_fft,unset_mpi_enreg_fft xmpi_sum
SOURCE
70 #if defined HAVE_CONFIG_H 71 #include "config.h" 72 #endif 73 74 #include "abi_common.h" 75 76 subroutine pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband1,iband2,ispinor1,ispinor2,istwf_k,kg_diel,& 77 & lmax_diel,mgfftdiel,natom,nband,ndiel4,ndiel5,ndiel6,& 78 & ngfftdiel,npwdiel,nspinor,ntypat,optreal,& 79 & pawang,pawtab,ph3d_diel,typat,wfprod,wfraug, & 80 & mpi_atmtab,comm_atom,comm_fft,me_g0,paral_kgb,distribfft) ! optional arguments (parallelism) 81 82 use defs_basis 83 use defs_abitypes 84 use m_profiling_abi 85 use m_distribfft 86 use m_errors 87 use m_xmpi, only : xmpi_comm_self,xmpi_sum 88 89 use m_pawang, only : pawang_type 90 use m_pawtab, only : pawtab_type 91 use m_pawcprj, only : pawcprj_type 92 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 93 use m_mpinfo, only : set_mpi_enreg_fft, unset_mpi_enreg_fft 94 95 !This section has been created automatically by the script Abilint (TD). 96 !Do not modify the following lines by hand. 97 #undef ABI_FUNC 98 #define ABI_FUNC 'pawsushat' 99 use interfaces_51_manage_mpi 100 use interfaces_53_ffts 101 !End of the abilint section 102 103 implicit none 104 105 !Arguments --------------------------------------------- 106 !scalars 107 integer,intent(in) :: iband1,iband2,ispinor1,ispinor2,istwf_k,lmax_diel,mgfftdiel 108 integer,intent(in) :: natom,nband,ndiel4,ndiel5,ndiel6,npwdiel,nspinor 109 integer,intent(in) :: ntypat,optreal 110 integer,optional,intent(in) :: me_g0,comm_atom,comm_fft,paral_kgb 111 type(distribfft_type),optional,intent(in),target :: distribfft 112 type(pawang_type),intent(in) :: pawang 113 !arrays 114 integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2) 115 integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom) 116 integer,optional,target,intent(in) :: mpi_atmtab(:) 117 real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat) 118 real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom) 119 real(dp),intent(inout) :: wfprod(2,npwdiel*(1-optreal)) 120 real(dp),intent(inout) :: wfraug(2,ndiel4,ndiel5,ndiel6*optreal) 121 type(pawcprj_type),intent(in) :: cprj_k(natom,nspinor*nband) 122 type(pawtab_type),intent(in) :: pawtab(ntypat) 123 124 !Local variables --------------------------------------- 125 !scalars 126 integer :: cplex,iatm,iatom,iatom_tot,ibsp1,ibsp2,ierr,il,ilmn,ils,ilslm,ipw 127 integer :: itypat,j0lmn,jlmn,klm,klmn,lmax,lmin,mm,my_comm_atom,my_comm_fft,my_natom,paral_kgb_fft,tim_fourwf 128 real(dp) :: phil1,phil2,sgn,weight_dum,wf1,wf2 129 logical :: my_atmtab_allocated,parity,paral_atom 130 type(distribfft_type),pointer :: my_distribfft 131 type(mpi_type) :: mpi_enreg_fft 132 !arrays 133 integer,pointer :: my_atmtab(:) 134 real(dp) :: ro(2),ro_ql(2) 135 real(dp),allocatable :: dummy(:,:),wfprod_paw(:,:),wfraug_paw(:,:,:,:) 136 137 ! ************************************************************************* 138 139 DBG_ENTER("COLL") 140 141 if (present(comm_fft)) then 142 if ((.not.present(paral_kgb)).or.(.not.present(me_g0))) then 143 MSG_BUG('Need paral_kgb and me_g0 with comm_fft !') 144 end if 145 end if 146 147 !Set up parallelism over atoms 148 paral_atom=(present(comm_atom).and.(my_natom/=natom)) 149 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 150 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 151 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom) 152 my_natom=natom;if (paral_atom) my_natom=size(my_atmtab) 153 154 cplex=1;if (istwf_k>1) cplex=2 155 ABI_ALLOCATE(wfprod_paw,(2,npwdiel)) 156 wfprod_paw(:,:)=zero 157 ibsp1=(iband1-1)*nspinor+ispinor1 158 ibsp2=(iband2-1)*nspinor+ispinor2 159 160 !------------------------------------------------------------------------ 161 !----- Loop over atoms 162 !------------------------------------------------------------------------ 163 do iatom=1,my_natom 164 iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom) 165 iatm=atindx(iatom_tot) 166 itypat=typat(iatom_tot) 167 168 ! ------------------------------------------------------------------------ 169 ! ----- Loop over ij channels (basis components) 170 ! ------------------------------------------------------------------------ 171 do jlmn=1,pawtab(itypat)%lmn_size 172 j0lmn=jlmn*(jlmn-1)/2 173 do ilmn=1,jlmn 174 klmn=j0lmn+ilmn 175 klm =pawtab(itypat)%indklmn(1,klmn) 176 lmin=pawtab(itypat)%indklmn(3,klmn) 177 lmax=pawtab(itypat)%indklmn(4,klmn) 178 179 ro(1)=cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn) 180 if (cplex==2) then 181 ro(1)=ro(1)+cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn) 182 ro(2)=cprj_k(iatm,ibsp1)%cp(2,ilmn)*cprj_k(iatm,ibsp2)%cp(1,jlmn) & 183 & -cprj_k(iatm,ibsp1)%cp(1,ilmn)*cprj_k(iatm,ibsp2)%cp(2,jlmn) 184 end if 185 ro(1:cplex)=ro(1:cplex)*pawtab(itypat)%dltij(klmn) 186 187 do ils=lmin,lmax,2 188 il=mod(ils,4);parity=(mod(il,2)==0) 189 sgn=one;if (il>1) sgn=-one 190 191 do mm=-ils,ils 192 ilslm=ils*ils+ils+mm+1 193 if (pawang%gntselect(ilslm,klm)>0) then 194 195 ro_ql(1:cplex)=pawtab(itypat)%qijl(ilslm,klmn)*ro(1:cplex) 196 197 ! Compute: Sum_{ijR} [ cpi* cpj qij^l (-i)^l g_l(g) S_lm(g) ] 198 199 if (cplex==1) then 200 if (parity) then 201 do ipw=1,npwdiel 202 phil1= sgn*ph3d_diel(1,ipw,iatm) ! (i)^l.exp(i.g.R) 203 phil2= sgn*ph3d_diel(2,ipw,iatm) 204 wf1= phil1*ro_ql(1) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 205 wf2=-phil2*ro_ql(1) 206 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 207 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 208 end do 209 else 210 do ipw=1,npwdiel 211 phil1=-sgn*ph3d_diel(2,ipw,iatm) ! (i)^l.exp(i.g.R) 212 phil2= sgn*ph3d_diel(1,ipw,iatm) 213 wf1= phil1*ro_ql(1) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 214 wf2=-phil2*ro_ql(1) 215 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 216 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 217 end do 218 end if 219 220 else 221 222 if (parity) then 223 do ipw=1,npwdiel 224 phil1= sgn*ph3d_diel(1,ipw,iatm) ! (i)^l.exp(i.g.R) 225 phil2= sgn*ph3d_diel(2,ipw,iatm) 226 wf1=phil1*ro_ql(1)+phil2*ro_ql(2) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 227 wf2=phil1*ro_ql(2)-phil2*ro_ql(1) 228 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 229 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 230 end do 231 else 232 do ipw=1,npwdiel 233 phil1=-sgn*ph3d_diel(2,ipw,iatm) ! (i)^l.exp(i.g.R) 234 phil2= sgn*ph3d_diel(1,ipw,iatm) 235 wf1=phil1*ro_ql(1)+phil2*ro_ql(2) ! cpi* cpj qij^l (-i)^l.exp(-i.g.R) 236 wf2=phil1*ro_ql(2)-phil2*ro_ql(1) 237 wfprod_paw(1,ipw)=wfprod_paw(1,ipw)+wf1*gylmg_diel(ipw,ilslm,itypat) 238 wfprod_paw(2,ipw)=wfprod_paw(2,ipw)+wf2*gylmg_diel(ipw,ilslm,itypat) 239 end do 240 end if 241 242 end if 243 end if 244 end do 245 end do 246 247 ! ----- End loop over ij channels 248 end do 249 end do 250 251 ! ----- End loop over atoms 252 end do 253 254 !Reduction in case of parallelism over atoms 255 if (paral_atom) then 256 call xmpi_sum(wfprod_paw,my_comm_atom,ierr) 257 end if 258 259 if (optreal==0) then 260 261 ! === Output in reciprocal space 262 wfprod(:,:)=wfprod(:,:)+wfprod_paw(:,:) 263 264 else 265 ! === Output in reciprocal space 266 tim_fourwf=17;weight_dum=0 267 ! Create fake mpi_enreg to wrap fourdp 268 if (present(distribfft)) then 269 my_distribfft => distribfft 270 else 271 ABI_DATATYPE_ALLOCATE(my_distribfft,) 272 call init_distribfft_seq(my_distribfft,'c',ngfftdiel(2),ngfftdiel(3),'fourwf') 273 end if 274 call initmpi_seq(mpi_enreg_fft) 275 ABI_DATATYPE_DEALLOCATE(mpi_enreg_fft%distribfft) 276 if (present(comm_fft)) then 277 call set_mpi_enreg_fft(mpi_enreg_fft,comm_fft,my_distribfft,me_g0,paral_kgb) 278 my_comm_fft=comm_fft;paral_kgb_fft=paral_kgb 279 else 280 my_comm_fft=xmpi_comm_self;paral_kgb_fft=0; 281 mpi_enreg_fft%distribfft => my_distribfft 282 end if 283 ! do FFT 284 ABI_ALLOCATE(wfraug_paw,(2,ndiel4,ndiel5,ndiel6)) 285 call fourwf(1,dummy,wfprod_paw,dummy,wfraug_paw,gbound_diel,gbound_diel,& 286 & istwf_k,kg_diel,kg_diel,mgfftdiel,mpi_enreg_fft,1,ngfftdiel,1,npwdiel,& 287 & ndiel4,ndiel5,ndiel6,0,paral_kgb_fft,tim_fourwf,weight_dum,weight_dum) 288 wfraug(:,:,:,:)=wfraug(:,:,:,:)+wfraug_paw(:,:,:,:) 289 ABI_DEALLOCATE(wfraug_paw) 290 call unset_mpi_enreg_fft(mpi_enreg_fft) 291 if (.not.present(distribfft)) then 292 call destroy_distribfft(my_distribfft) 293 ABI_DATATYPE_DEALLOCATE(my_distribfft) 294 end if 295 end if 296 297 ABI_DEALLOCATE(wfprod_paw) 298 299 !Destroy atom table used for parallelism 300 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 301 302 DBG_EXIT("COLL") 303 304 end subroutine pawsushat