TABLE OF CONTENTS
ABINIT/susk [ Functions ]
NAME
susk
FUNCTION
Compute the contribution of one k point to the susceptibility matrix from input wavefunctions, band occupations, and k point wts. Include the usual sum-over-state terms, but also the corrections due to the change of the Fermi level in the metallic case, as well as implicit sum over higher lying conduction states, thanks to the closure relation (referred to as an extrapolation). Compared to the routine suskmm, there is no particular attention to the use of the memory, so the code is simpler.
COPYRIGHT
Copyright (C) 1998-2017 ABINIT group (XG). 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 bdtot_index=index for the number of the band cg(2,mcg)=wfs in G space 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. doccde(mband*nkpt*nsppol)=derivative of occupancies wrt the energy for each band and k point eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree) extrap: if==1, the closure relation (an extrapolation) must be used gbound(2*mgfftdiel+8,2)=G sphere boundary for going from WF sphere to medium size FFT grid gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for going from medium size FFT grid to small sphere. gylmg_diel(npwdiel,lmax_diel,ntypat*usepaw)= -PAW only- Fourier transform of g_l(r).Y_ml(r) shape functions for dielectric matrix icg=index for cg ikpt=number of the k point isp=number of the current spin istwf_k=input option parameter that describes the storage of wfs kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. kg_k(3,npw_k)=coordinates of planewaves in basis sphere. lmax_diel=1+max. value of l angular momentum used for dielectric matrix mband=maximum number of bands mcg=dimension of cg mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix mpi_enreg=information about MPI parallelization natom=number of atoms in cell nband_k=number of bands at this k point for that spin polarization ndiel4,ndiel5,ndiel6= FFT dimensions, modified to avoid cache trashing neglect_pawhat=1 if PAW contribution from hat density (compensation charge) has to be neglected (to be used when only an estimation of suscep. matrix has to be evaluated, i.e. for SCF precondictioning) ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix, see ~abinit/doc/input_variables/vargs.htm#ngfft nkpt=number of k points npwdiel=third and fifth dimension of the susmat array. npw_k=number of plane waves at this k point nspden=number of spin-density components nspden_eff=number of spin-density components actually computed in sussceptibility nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized ntypat=number of types of atoms in unit cell. occ(mband*nkpt*nsppol)= occupation numbers for each band (usually 2.0) at each k point occopt=option for occupancies occ_deavg(mband)=factor for extrapolation (occup. divided by an energy gap) 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 ucvol=unit cell volume (Bohr**3) usepaw=flag for PAW wtk(nkpt)=k point weights (they sum to 1.0)
OUTPUT
(see side effects)
SIDE EFFECTS
These quantities are accumulated in this routine: drhode(2,npwdiel,nspden_eff)=weighted density, needed to compute the effect of change of fermi energy rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)=density-like array, needed for the extrapolation procedure. sumdocc=sum of weighted occupation numbers, needed to compute the effect of change of fermi energy susmat(2,npwdiel,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space
NOTES
Band-fft parallel treatment: Each processor will treat his own band, but susmat will be known by all. This means that cg will not have the same meaning in sequential or parallel mode. In parallel mode, it will contain the set of all bands treated by the currrent processor. To achieve this, the argument cg has been replaced by cg_mpi, with the "target" attribute. In sequential mode, the pointer cg will point towards cg_mpi. In parallel mode, cg will point to a new array cg_local, containing the bands treated by the currrent processor. This allows to minimize the overhead incurred by the parallelization of the sequential version. A similar treatment is performed on kg_k, npw_k. A future version might have objects like kg_k_gather as arguments, instead of computing them. This is in slight violation of programming rules, but I think it is safe, since the pointers remain local GZ
PARENTS
suscep_stat
CHILDREN
destroy_mpi_enreg,fourwf,init_distribfft_seq,initmpi_seq,pawsushat sphereboundary,timab,xmpi_allgather,xmpi_allgatherv,xmpi_alltoallv xmpi_sum
SOURCE
115 #if defined HAVE_CONFIG_H 116 #include "config.h" 117 #endif 118 119 #include "abi_common.h" 120 121 122 subroutine susk(atindx,bdtot_index,cg_mpi,cprj_k,doccde,drhode,eigen,extrap,gbound,& 123 & gbound_diel,gylmg_diel,icg_mpi,ikpt,isp,istwf_k,kg_diel,kg_k_mpi,& 124 & lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 125 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 126 & npwdiel,npw_k_mpi,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,& 127 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 128 & susmat,typat,ucvol,usepaw,wtk) 129 130 use defs_basis 131 use defs_abitypes 132 use m_profiling_abi 133 use m_xmpi 134 135 use m_pawang, only : pawang_type 136 use m_pawtab, only : pawtab_type 137 use m_pawcprj, only : pawcprj_type 138 use m_mpinfo, only : destroy_mpi_enreg 139 140 !This section has been created automatically by the script Abilint (TD). 141 !Do not modify the following lines by hand. 142 #undef ABI_FUNC 143 #define ABI_FUNC 'susk' 144 use interfaces_18_timing 145 use interfaces_51_manage_mpi 146 use interfaces_52_fft_mpi_noabirule 147 use interfaces_53_ffts 148 use interfaces_65_paw 149 !End of the abilint section 150 151 implicit none 152 153 !Arguments ------------------------------------ 154 !This type is defined in defs_mpi 155 !scalars 156 integer,intent(in) :: bdtot_index,extrap,ikpt,isp,istwf_k,lmax_diel,mband,mcg 157 integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat 158 integer,intent(in) :: nkpt,npwdiel,nspden,nspden_eff,nspinor,nsppol 159 integer,intent(in) :: ntypat,occopt,usepaw 160 integer,intent(in),target :: icg_mpi,npw_k_mpi 161 real(dp),intent(in) :: ucvol 162 real(dp),intent(inout) :: sumdocc 163 type(MPI_type),intent(in) :: mpi_enreg 164 type(pawang_type),intent(in) :: pawang 165 !arrays 166 integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2) 167 integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom) 168 integer,intent(in),target :: kg_k_mpi(3,npw_k_mpi) 169 integer,intent(inout) :: gbound(2*mgfftdiel+8,2) 170 integer,pointer :: kg_k(:,:) 171 real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol) 172 real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw) 173 real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband) 174 real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt) 175 real(dp),intent(in),target :: cg_mpi(2,mcg) 176 real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff) 177 real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor) 178 real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 179 type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw) 180 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 181 182 !Local variables------------------------------- 183 ! real(dp), allocatable :: cg_disk(:,:) 184 !Local variables for MPI 185 !scalars 186 integer :: blocksize,i1,i2,i3,iband,iband_loc,ibd1,ibd2,ibdblock,ier 187 integer :: iproc,iproc_fft,ipw,ipw1,ipw2,isp1,isp2,ispinor,iwf,jsp,me_bandfft 188 integer :: nbdblock,ndatarecv,ndiel1,ndiel2,ndiel3 189 integer :: sizemax_per_proc,spaceComm,testocc,tim_fourwf 190 integer,pointer :: icg,npw_k 191 integer,target :: icg_loc=0,npw_k_loc,npw_tot 192 real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2 193 type(MPI_type) :: mpi_enreg_diel 194 !arrays 195 integer,allocatable :: band_loc(:),kg_k_gather(:,:),npw_per_proc(:),rdispls(:) 196 integer,allocatable :: rdispls_all(:),rdisplsloc(:),recvcounts(:) 197 integer,allocatable :: recvcountsloc(:),sdispls(:),sdisplsloc(:),sendcounts(:) 198 integer,allocatable :: sendcountsloc(:),susmat_mpi(:,:,:) 199 integer,allocatable,target :: kg_k_gather_all(:,:) 200 real(dp) :: tsec(2) 201 real(dp),allocatable :: cwavef(:,:),cwavef_alltoall(:,:) 202 real(dp),allocatable :: cwavef_alltoall_gather(:,:),dummy(:,:),rhoaug(:,:,:) 203 real(dp),allocatable :: wfprod(:,:),wfraug(:,:,:,:),wfrspa(:,:,:,:,:,:) 204 real(dp),allocatable,target :: cg_local(:,:) 205 real(dp),pointer :: cg(:,:) 206 logical,allocatable :: treat_band(:) 207 208 ! ************************************************************************* 209 210 !DEBUG 211 !write(std_out,*)' susk : enter ' 212 !if(.true.)stop 213 !ENDDEBUG 214 215 call timab(750,1,tsec) 216 call timab(751,1,tsec) 217 218 ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3) 219 220 !The dielectric stuff is performed in sequential mode. 221 !Set mpi_enreg_diel accordingly 222 call initmpi_seq(mpi_enreg_diel) 223 call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 224 me_bandfft=xmpi_comm_rank(mpi_enreg%comm_bandfft) 225 226 testocc=1 227 !DEBUG 228 !write(std_out,*)' susk : set testocc to 0 ' 229 !testocc=0 230 !write(std_out,*)' susk : set extrap to 0 ' 231 !extrap=0 232 !ENDDEBUG 233 234 !Allocations, initializations 235 ABI_ALLOCATE(rhoaug,(ndiel4,ndiel5,ndiel6)) 236 ABI_ALLOCATE(wfraug,(2,ndiel4,ndiel5,ndiel6)) 237 ABI_ALLOCATE(wfprod,(2,npwdiel)) 238 ABI_ALLOCATE(wfrspa,(2,ndiel4,ndiel5,ndiel6,nspinor,mband)) 239 ABI_ALLOCATE(dummy,(2,1)) 240 wfrspa(:,:,:,:,:,:)=zero 241 ABI_ALLOCATE(treat_band,(nband_k)) 242 treat_band(:)=.true. 243 isp1=isp;isp2=isp 244 if (nspden_eff==2.and.nspinor==2) isp2=isp+1 245 246 !BAND-FFT parallelism 247 if (mpi_enreg%paral_kgb==1) then 248 treat_band(:)=.false. 249 ! We gather the wavefunctions treated by this proc in cg_local 250 spaceComm=mpi_enreg%comm_band 251 blocksize=mpi_enreg%nproc_band 252 nbdblock=nband_k/blocksize 253 ABI_ALLOCATE(sdispls,(blocksize)) 254 ABI_ALLOCATE(sdisplsloc,(blocksize)) 255 ABI_ALLOCATE(sendcounts,(blocksize)) 256 ABI_ALLOCATE(sendcountsloc,(blocksize)) 257 ABI_ALLOCATE(rdispls,(blocksize)) 258 ABI_ALLOCATE(rdisplsloc,(blocksize)) 259 ABI_ALLOCATE(recvcounts,(blocksize)) 260 ABI_ALLOCATE(recvcountsloc,(blocksize)) 261 ! First gather the kg_k in kg_k_gather_all 262 npw_k_loc=npw_k_mpi 263 call xmpi_allgather(npw_k_loc,recvcounts,spaceComm,ier) 264 rdispls(1)=0 265 do iproc=2,blocksize 266 rdispls(iproc)=rdispls(iproc-1)+recvcounts(iproc-1) 267 end do 268 ndatarecv=rdispls(blocksize)+recvcounts(blocksize) 269 ABI_ALLOCATE(kg_k_gather,(3,ndatarecv)) 270 recvcountsloc(:)=recvcounts(:)*3 271 rdisplsloc(:)=rdispls(:)*3 272 call xmpi_allgatherv(kg_k_mpi,3*npw_k_loc,kg_k_gather,recvcountsloc(:),rdisplsloc,spaceComm,ier) 273 ABI_ALLOCATE(npw_per_proc,(mpi_enreg%nproc_fft)) 274 ABI_ALLOCATE(rdispls_all,(mpi_enreg%nproc_fft)) 275 spaceComm=mpi_enreg%comm_fft 276 call xmpi_allgather(ndatarecv,npw_per_proc,spaceComm,ier) 277 rdispls_all(1)=0 278 do iproc=2,mpi_enreg%nproc_fft 279 rdispls_all(iproc)=rdispls_all(iproc-1)+npw_per_proc(iproc-1) 280 end do 281 npw_tot=rdispls_all(mpi_enreg%nproc_fft)+npw_per_proc(mpi_enreg%nproc_fft) 282 ABI_ALLOCATE(kg_k_gather_all,(3,npw_tot)) 283 call xmpi_allgatherv(kg_k_gather,3*ndatarecv,kg_k_gather_all,3*npw_per_proc(:),3*rdispls_all,spaceComm,ier) 284 ! At this point kg_k_gather_all contains all the kg 285 if(allocated(cwavef)) then 286 ABI_DEALLOCATE(cwavef) 287 end if 288 ABI_ALLOCATE(cwavef,(2,npw_k_loc*nspinor*blocksize)) 289 sizemax_per_proc=nband_k/(mpi_enreg%nproc_band*mpi_enreg%nproc_fft)+1 290 ABI_ALLOCATE(band_loc,(sizemax_per_proc)) 291 ABI_ALLOCATE(cg_local,(2,sizemax_per_proc*npw_tot*nspinor)) 292 iband_loc=0 293 do ibdblock=1,nbdblock 294 cwavef(:,1:npw_k_loc*nspinor*blocksize)=& 295 & cg_mpi(:,1+(ibdblock-1)*npw_k_loc*nspinor*blocksize+icg_mpi:ibdblock*npw_k_loc*nspinor*blocksize+icg_mpi) 296 sendcounts(:)=npw_k_loc 297 do iproc=1,blocksize 298 sdispls(iproc)=(iproc-1)*npw_k_loc 299 end do 300 ABI_ALLOCATE(cwavef_alltoall,(2,ndatarecv*nspinor)) 301 recvcountsloc(:)=recvcounts(:)*2*nspinor 302 rdisplsloc(:)=rdispls(:)*2*nspinor 303 sendcountsloc(:)=sendcounts(:)*2*nspinor 304 sdisplsloc(:)=sdispls(:)*2*nspinor 305 call timab(547,1,tsec) 306 spaceComm=mpi_enreg%comm_band 307 call xmpi_alltoallv(cwavef,sendcountsloc,sdisplsloc,cwavef_alltoall,recvcountsloc,rdisplsloc,spaceComm,ier) 308 call timab(547,2,tsec) 309 ABI_ALLOCATE(cwavef_alltoall_gather,(2,npw_tot*nspinor)) 310 blocksize=mpi_enreg%nproc_band 311 spaceComm=mpi_enreg%comm_fft 312 call xmpi_allgatherv(cwavef_alltoall,2*nspinor*ndatarecv,cwavef_alltoall_gather,& 313 & 2*nspinor*npw_per_proc,2*nspinor*rdispls_all,spaceComm,ier) 314 iproc_fft=modulo(ibdblock-1,mpi_enreg%nproc_fft) 315 if(mpi_enreg%me_fft==iproc_fft) then !All nproc_band procs of index me_fft will treat these bands 316 iband_loc=iband_loc+1 317 iband=1+mpi_enreg%me_band+mpi_enreg%nproc_band*mpi_enreg%me_fft+(iband_loc-1)*mpi_enreg%nproc_fft*mpi_enreg%nproc_band 318 treat_band(iband)=.true. 319 band_loc(iband_loc)=iband 320 cg_local(:,1+(iband_loc-1)*npw_tot*nspinor:iband_loc*npw_tot*nspinor)=cwavef_alltoall_gather(:,1:npw_tot*nspinor) 321 end if 322 ABI_DEALLOCATE(cwavef_alltoall_gather) 323 ABI_DEALLOCATE(cwavef_alltoall) 324 end do 325 ! On exit: 326 ! npw_tot will be npw 327 ! kg_k_gather_all will be kg_k 328 ! cg_local will be cg 329 ! icg will be zero 330 npw_k=>npw_tot 331 kg_k=>kg_k_gather_all(:,:) 332 cg=>cg_local(:,:) 333 icg=>icg_loc 334 call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k) 335 ABI_DEALLOCATE(npw_per_proc) 336 ABI_DEALLOCATE(rdispls_all) 337 ABI_DEALLOCATE(sendcounts) 338 ABI_DEALLOCATE(recvcounts) 339 ABI_DEALLOCATE(sdispls) 340 ABI_DEALLOCATE(rdispls) 341 ABI_DEALLOCATE(sendcountsloc) 342 ABI_DEALLOCATE(sdisplsloc) 343 ABI_DEALLOCATE(recvcountsloc) 344 ABI_DEALLOCATE(rdisplsloc) 345 ABI_DEALLOCATE(kg_k_gather) 346 ABI_DEALLOCATE(cwavef) 347 ! Because they will be summed over all procs, and arrive on input, rescale drhode and rhoextrap 348 if(occopt>=3)drhode(:,:,isp1:isp2)=drhode(:,:,isp1:isp2)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 349 if(extrap==1)rhoextrap(:,:,:,:)=rhoextrap(:,:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 350 do i1=isp1,isp2 351 susmat(:,:,i1,:,i1)=susmat(:,:,i1,:,i1)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 352 end do 353 354 ! No BAND-FFT parallelism 355 else ! use argument variables 356 cg=>cg_mpi 357 kg_k=>kg_k_mpi 358 npw_k=>npw_k_mpi 359 icg=>icg_mpi 360 end if 361 iband_loc=0 362 363 call timab(751,2,tsec) 364 call timab(752,1,tsec) 365 366 !Loop over bands to fft and store Fourier transform of wavefunction 367 ABI_ALLOCATE(cwavef,(2,npw_k)) 368 do iband=1,nband_k 369 if(.not. treat_band(iband)) cycle ! I am not treating this band (only for the parallel case) 370 iband_loc=iband_loc+1 371 372 ! Loop on spinorial components 373 do ispinor=1,nspinor 374 iwf=(ispinor-1)*npw_k+(iband_loc-1)*npw_k*nspinor+icg 375 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 376 377 ! Obtain Fourier transform in fft box 378 tim_fourwf=8 379 cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf) 380 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 381 & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,& 382 & 0,mpi_enreg_diel%paral_kgb,tim_fourwf,weight,weight) 383 384 wfrspa(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:) 385 386 if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then 387 ! In the case of metallic occupation, or if the extrapolation 388 ! over higher bands is included, must compute the 389 ! Fourier transform of the density of each band, then 390 ! generate the part of the susceptibility matrix due 391 ! varying occupation numbers. 392 393 weight=-two*occ_deavg(iband)*wtk(ikpt)/ucvol 394 do i3=1,ndiel3 395 do i2=1,ndiel2 396 do i1=1,ndiel1 397 wfraug(1,i1,i2,i3)=wfraug(1,i1,i2,i3)**2+wfraug(2,i1,i2,i3)**2 398 wfraug(2,i1,i2,i3)=zero 399 end do 400 end do 401 ! If extrapolation, accumulate density in real space 402 if(extrap==1.and.usepaw==0)then 403 do i2=1,ndiel2 404 do i1=1,ndiel1 405 rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3) 406 end do 407 end do 408 end if 409 end do 410 411 ! In case of PAW, add compensation charge contribution 412 if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then 413 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,& 414 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 415 & ngfftdiel,npwdiel,nspinor,ntypat,1,& 416 & pawang,pawtab,ph3d_diel,typat,dummy,wfraug,& 417 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 418 rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:) 419 end if 420 421 ! Performs the Fourier Transform of the density of the band, 422 ! and store it in wfprod 423 tim_fourwf=9 424 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 425 & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 426 & ndiel4,ndiel5,ndiel6,3,mpi_enreg_diel%paral_kgb,tim_fourwf,weight,weight) 427 ! In case of PAW, add compensation charge contribution if not already done 428 if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then 429 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,& 430 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 431 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 432 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 433 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 434 end if 435 436 ! Perform now the summation of terms related to direct change of eigenvalues 437 ! or extrapolation over higher bands 438 wght1=zero ; wght2=zero 439 if(occopt>=3 .and. testocc==1) wght1=doccde(iband+bdtot_index)*wtk(ikpt)/ucvol 440 if(extrap==1) wght2=two*occ_deavg(iband)*wtk(ikpt)/ucvol 441 weight=wght1+wght2 442 443 if (abs(weight)>tol12) then 444 do ipw2=1,npwdiel 445 ! Only fills lower half of the matrix (here, the susceptibility matrix) 446 ! Note that wfprod of the first index must behave like a density, 447 ! so that it is used as generated by fourwf, while wfprod of the 448 ! second index will be implicitely used to make a scalar product 449 ! with a potential change, meaning that its complex conjugate must be 450 ! used. This explains the following signs... 451 do ipw1=ipw2,npwdiel 452 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 453 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 454 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 455 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 456 end do 457 end do 458 end if 459 460 if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then 461 ! Accumulate product of band densities by their doccde, for the 462 ! computation of the effect of change of Fermi level. 463 do ipw=1,npwdiel 464 drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1 465 drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1 466 end do 467 ! Also accumulate weighted sum of doccde 468 sumdocc=sumdocc+wght1 469 end if 470 471 ! End condition of metallic occupancies or extrapolation 472 end if 473 474 ! End loop on spinorial components 475 end do 476 ! End loop on iband 477 end do 478 479 call timab(752,2,tsec) 480 call timab(753,1,tsec) 481 482 ABI_DEALLOCATE(cwavef) 483 484 !Stuff for parallelism (bands-FFT) 485 if(mpi_enreg%paral_kgb==1) then 486 call xmpi_sum(wfrspa,mpi_enreg%comm_bandfft,ier) 487 if(occopt>=3) then 488 call xmpi_sum(drhode(:,:,isp1:isp2),mpi_enreg%comm_bandfft,ier) 489 end if 490 if(extrap==1) then 491 call xmpi_sum(rhoextrap,mpi_enreg%comm_bandfft,ier) 492 end if 493 if(occopt>=3) then 494 call xmpi_sum(sumdocc,mpi_enreg%comm_bandfft,ier) 495 end if 496 ABI_ALLOCATE(susmat_mpi,(2,npwdiel,npwdiel)) 497 do i1=isp1,isp2 498 susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1) 499 call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier) 500 susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 501 end do 502 ABI_DEALLOCATE(susmat_mpi) 503 end if 504 call timab(753,2,tsec) 505 506 !-- Wavefunctions have been generated in real space ------------------------ 507 !-- Now, compute product of wavefunctions for different bands -------------- 508 call timab(754,1,tsec) 509 !if (occopt<3) then 510 tolocc=1.0d-3 511 !else 512 !tolocc=1.0d-8 513 !end if 514 iproc=-1 515 516 if(nband_k>1)then 517 do ibd1=1,nband_k-1 518 do ibd2=ibd1+1,nband_k 519 iproc=iproc+1 520 if(modulo(iproc,mpi_enreg%nproc_fft*mpi_enreg%nproc_band) /= me_bandfft) cycle 521 ! If the occupation numbers are sufficiently different, or 522 ! if extrapolation is used and the corresponding factor is not zero, 523 ! then there is a contribution 524 occdiff=occ(ibd1+bdtot_index)-occ(ibd2+bdtot_index) 525 if( abs(occdiff)>tolocc .or. & 526 & ( extrap==1 .and. & 527 & ( abs(occ_deavg(ibd1)) + abs(occ_deavg(ibd2)) ) >tolocc ) & 528 & ) then 529 530 eigdiff=eigen(ibd1+bdtot_index)-eigen(ibd2+bdtot_index) 531 ! DEBUG 532 ! write(std_out,*)' susk : contribution from bands',ibd1,ibd2 533 ! write(std_out,*)' occ diff =',occdiff 534 ! write(std_out,*)' eig diff =',eigdiff 535 ! ENDDEBUG 536 537 ! Loop on spinorial components 538 do ispinor=1,nspinor 539 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 540 541 ! Store the contribution in wfraug 542 do i3=1,ndiel3 543 do i2=1,ndiel2 544 do i1=1,ndiel1 545 wfraug(1,i1,i2,i3)=wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)& 546 & +wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2) 547 wfraug(2,i1,i2,i3)=wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)& 548 & -wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2) 549 end do 550 end do 551 end do 552 553 ! Performs the Fourier Transform of the product, and store it in wfprod 554 tim_fourwf=19 555 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 556 & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 557 & ndiel4,ndiel5,ndiel6,3,mpi_enreg_diel%paral_kgb,tim_fourwf,weight,weight) 558 559 ! In case of PAW, add compensation charge contribution 560 if (usepaw==1.and.neglect_pawhat==0) then 561 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,& 562 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 563 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 564 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 565 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 566 end if 567 568 ! Perform now the summation 569 wght1=zero ; wght2=zero 570 if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol 571 if(extrap==1) wght2=(occ_deavg(ibd1)+occ_deavg(ibd2)) * two*wtk(ikpt)/ucvol 572 weight=wght1+wght2 573 574 ! DEBUG 575 ! write(std_out,*)' weight =',weight 576 ! norm=zero 577 ! do ipw=1,npwdiel 578 ! norm=norm+wfprod(1,ipw)**2+wfprod(2,ipw)**2 579 ! end do 580 ! write(std_out,*)' norm in reciprocal space =',norm 581 ! ENDDEBUG 582 583 if (abs(weight)>tol12) then 584 do ipw2=1,npwdiel 585 ! Only fills lower half of the matrix (here, the susceptibility matrix) 586 ! Note that wfprod of the first index must behave like a density, 587 ! so that it is used as generated by fourwf, while wfprod of the 588 ! second index will be implicitely used to make a scalar product 589 ! with a potential change, meaning that its complex conjugate must be 590 ! used. This explains the following signs... 591 do ipw1=ipw2,npwdiel 592 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 593 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 594 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 595 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 596 end do 597 end do 598 end if 599 600 ! End loop on spinorial components 601 end do 602 ! End condition of different occupation numbers or extrapolation 603 end if 604 ! End internal loop over bands 605 end do 606 ! End external loop over bands 607 end do 608 ! End condition of having more than one band 609 end if 610 611 call timab(754,2,tsec) 612 call timab(755,1,tsec) 613 614 if(mpi_enreg%paral_kgb==1) then 615 ABI_ALLOCATE(susmat_mpi,(2,npwdiel,npwdiel)) 616 do i1=isp1,isp2 617 susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1) 618 call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier) 619 susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:) 620 end do 621 ABI_DEALLOCATE(susmat_mpi) 622 ABI_DEALLOCATE(band_loc) 623 ABI_DEALLOCATE(treat_band) 624 ABI_DEALLOCATE(cg_local) 625 ABI_DEALLOCATE(kg_k_gather_all) 626 end if 627 628 call destroy_mpi_enreg(mpi_enreg_diel) 629 ABI_DEALLOCATE(dummy) 630 ABI_DEALLOCATE(rhoaug) 631 ABI_DEALLOCATE(wfprod) 632 ABI_DEALLOCATE(wfraug) 633 ABI_DEALLOCATE(wfrspa) 634 635 call timab(755,2,tsec) 636 call timab(750,2,tsec) 637 638 end subroutine susk