TABLE OF CONTENTS
ABINIT/m_suscep_stat [ Modules ]
NAME
m_suscep_stat
FUNCTION
Compute the susceptibility matrix
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (XG, AR, MB) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_suscep_stat 23 24 use defs_basis 25 use m_xmpi 26 use m_errors 27 use m_abicore 28 use m_distribfft 29 30 use defs_abitypes, only : MPI_type 31 use m_time, only : timab 32 use m_pawang, only : pawang_type 33 use m_pawtab, only : pawtab_type 34 use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, pawcprj_free 35 use m_mpinfo, only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle 36 use m_kg, only : ph1d3d 37 use m_gsphere, only : symg 38 use m_fftcore, only : sphereboundary 39 use m_fft, only : fftpac, fourwf 40 use m_spacepar, only : symrhg 41 use m_paw_finegrid, only : pawgylmg 42 use m_paw_nhat, only : pawsushat 43 44 implicit none 45 46 private
m_suscep_stat/suscep_stat [ Functions ]
[ Top ] [ m_suscep_stat ] [ Functions ]
NAME
suscep_stat
FUNCTION
Compute 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).
INPUTS
atindx(natom)=index table for atoms (see scfcv.f) atindx1(natom)=index table for atoms, inverse of atindx (see scfcv.f) cg(2,mcg)=wf in G space cprj(natom,mcg*usecprj)= wave functions projected with non-local projectors: cprj_nk(i)=<p_i|Cnk> where p_i is a non-local projector. dielar(7)=input parameters for dielectric matrix and susceptibility: diecut,dielng,diemac,diemix,diegap,dielam,diemixmag. dimcprj(natom*usepaw)=array of dimensions of array cprj (ordered by atom-type) 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) gbound_diel(2*mgfftdiel+8,2)=G sphere boundary for the dielectric matrix gprimd(3,3)=dimensional reciprocal space primitive translations irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data istwfk(nkpt)=input option parameter that describes the storage of wfs kg(3,mpw*mkmem)=reduced planewave coordinates. kg_diel(3,npwdiel)=reduced planewave coordinates for the dielectric matrix. lmax_diel=1+max. value of l angular momentum used for dielectric matrix mband=maximum number of bands mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol mcprj=size of projected wave-functions array (cprj) =nspinor*mband*mkmem*nsppol mgfftdiel=maximum size of 1D FFTs, for the computation of the dielectric matrix mkmem=number of k points treated by this node mpi_enreg=information about MPI parallelization mpw=maximum allowed value for npw natom=number of atoms in cell nband(nkpt*nsppol)=number of bands to be included in summation at each k point for each spin channel 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) nfftdiel=number of fft grid points for the computation of the diel matrix ngfftdiel(18)=contain all needed information about 3D FFT, for dielectric matrix, see ~abinit/doc/variables/vargs.htm#ngfft nkpt=number of k points npwarr(nkpt)=number of planewaves and boundary planewaves at each k, for going from the WF sphere to the medium size FFT grid. npwdiel=third and fifth dimension of the susmat array. nspden=number of spin-density components nspinor=number of spinorial components of the wavefunctions nsppol=1 for unpolarized, 2 for spin-polarized nsym=number of symmetry elements in group (at least 1 for identity) 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 pawang <type(pawang_type)>=paw angular mesh and related data pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4))=nonsymmorphic translation phases ph1ddiel(2,3*(2*mgfftdiel+1)*natom*usepaw)=one-dimensional structure factor information for the dielectric matrix rprimd(3,3)=dimensional real space primitive translations symafm(nsym)=(anti)ferromagnetic part of symmetry operations symrel(3,3,nsym)=symmetry matrices in real space (integers) tnons(3,nsym)=reduced nonsymmorphic translations (symrel and tnons are in terms of real space primitive translations) typat(natom)=type (integer) for each atom ucvol=unit cell volume (Bohr**3) unpaw=unit number for cprj PAW data (if used) usecprj= 1 if cprj array is stored in memory usepaw=flag for PAW usetimerev=1 if Time-Reversal symmetry has to be used when symmetrizing susceptibility wtk(nkpt)=k point weights (they sum to 1.0) ylmdiel(npwdiel,lmax_diel**2)= real spherical harmonics for each G and k point for the dielectric matrix
OUTPUT
susmat(2,npwdiel,nspden,npwdiel,nspden)= the susceptibility (or density-density response) matrix in reciprocal space
NOTES
Case of non-collinear magnetism: In principle, we should compute 16 susceptibility matrix: chi0-(s1,s2),(s3,s4) (where s1, s2, s3,and s4 are spin indexes)... But, for the time being, the susceptibility is only used to compute the dielectric matrix within RPA approximation; in this approximation, only four susceptibilities are non-zero: chi0-(s1,s1),(s3,s3). They are stored in susmat(:,ipw1,1:2,ipw2,1:2)
SOURCE
148 subroutine suscep_stat(atindx,atindx1,cg,cprj,dielar,dimcprj,doccde,& 149 & eigen,gbound_diel,gprimd,irrzondiel,istwfk,kg,& 150 & kg_diel,lmax_diel,& 151 & mband,mcg,mcprj,mgfftdiel,mkmem,mpi_enreg,mpw,natom,nband,& 152 & neglect_pawhat,nfftdiel,ngfftdiel,nkpt,npwarr,& 153 & npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occ,occopt,& 154 & pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,& 155 & susmat,symafm,symrel,tnons,typat,ucvol,unpaw,usecprj,usepaw,usetimerev,& 156 & wtk,ylmdiel) 157 158 !Arguments ------------------------------------ 159 !scalars 160 integer,intent(in) :: lmax_diel,mband,mcg,mcprj,mgfftdiel,mkmem,mpw,natom,neglect_pawhat 161 integer,intent(in) :: nfftdiel,nkpt,npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occopt 162 integer,intent(in) :: unpaw,usecprj,usepaw,usetimerev 163 real(dp),intent(in) :: ucvol 164 type(MPI_type),intent(in) :: mpi_enreg 165 type(pawang_type),intent(in) :: pawang 166 !arrays 167 integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom*usepaw) 168 integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2) 169 !no_abirules 170 !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise 171 integer,intent(in) :: irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),& 172 & istwfk(nkpt) 173 integer,intent(in) :: kg(3,mpw*mkmem),kg_diel(3,npwdiel),& 174 & nband(nkpt*nsppol),ngfftdiel(18) 175 integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym),typat(ntypat) 176 real(dp),intent(in) :: cg(2,mcg),dielar(7) 177 real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol) 178 real(dp),intent(in) :: gprimd(3,3),occ(mband*nkpt*nsppol) 179 !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise 180 real(dp),intent(in) :: phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),& 181 & tnons(3,nsym),wtk(nkpt) 182 real(dp),intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*usepaw),rprimd(3,3) 183 real(dp),intent(in) :: ylmdiel(npwdiel,lmax_diel**2) 184 real(dp),intent(out) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 185 type(pawcprj_type) :: cprj(natom,mcprj*usecprj) 186 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 187 188 !Local variables------------------------------- 189 !scalars 190 integer :: bdtot_index,diag,extrap,i1,i2,i3,iband,ibg,icg,ier,ierr 191 integer :: ifft,ii,ikg,ikpt,indx,iorder_cprj,ipw1,ipw2,isp,isp1,isp2 192 integer :: ispinor,istwf_k,isym,j1,j2,j3,jj,jsp,k1,k2,k3 193 integer :: my_nspinor,nband_k,nband_loc,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6 194 integer :: nkpg_diel,npw_k,npwsp,nspden_eff,nspden_tmp,nsym1,nsym2 195 integer :: spaceComm,t1,t2,testocc 196 real(dp) :: ai,ai2,ar,ar2,diegap,dielam,emax,invnsym 197 real(dp) :: invnsym1,invnsym2,phi1,phi12,phi2,phr1,phr12 198 real(dp) :: phr2,sumdocc,weight 199 logical :: antiferro 200 character(len=500) :: message 201 type(MPI_type) :: mpi_enreg_diel 202 203 !arrays 204 integer,allocatable :: gbound(:,:),kg_k(:,:),sym_g(:,:) 205 integer,allocatable :: tmrev_g(:) 206 real(dp) :: kpt_diel(3,1),tsec(2) 207 real(dp),allocatable :: drhode(:,:,:),drhode_wk(:,:,:) 208 real(dp),allocatable :: eig_diel(:),gylmg_diel(:,:,:),kpg_dum(:,:) 209 real(dp),allocatable :: occ_deavg(:),ph3d_diel(:,:,:),phdiel(:,:,:) 210 real(dp),allocatable :: phkxred_diel(:,:),rhoextrap(:,:,:,:),rhoextrg(:,:) 211 real(dp),allocatable :: rhoextrr(:,:),sush(:),sussum(:),susvec(:,:,:) 212 real(dp),allocatable :: suswk(:,:,:,:),zhpev1(:,:),zhpev2(:) 213 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_loc(:,:) 214 215 ! ************************************************************************* 216 217 call timab(740,1,tsec) 218 call timab(741,1,tsec) 219 220 ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!") 221 222 223 !----- Initialisations ----------------------------------------------------------- 224 !--------------------------------------------------------------------------------- 225 226 if (usecprj==0.and.usepaw==1) then 227 write (message,'(3a)')& 228 & ' cprj datastructure must be allocated !',ch10,& 229 & ' Action: change pawusecp input keyword.' 230 ABI_ERROR(message) 231 end if 232 233 if (mpi_enreg%paral_spinor==1) then 234 message = ' not yet allowed for parallelization over spinors !' 235 ABI_ERROR(message) 236 end if 237 238 !Init mpicomm 239 if(mpi_enreg%paral_kgb==1) then 240 spaceComm=mpi_enreg%comm_kpt 241 else 242 spaceComm=mpi_enreg%comm_cell 243 end if 244 245 !The dielectric stuff is performed in sequential mode. 246 !Set mpi_enreg_diel accordingly 247 call initmpi_seq(mpi_enreg_diel) 248 call init_distribfft_seq(MPI_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 249 250 !testocc to be taken away 251 testocc=1 252 253 my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor) 254 255 iorder_cprj=0 ! order for the cprj reading... 256 257 !Initialize some scalar quantities 258 antiferro=(nsppol==1.and.nspden==2) 259 nspden_eff=min(max(nsppol,nspden),2) ! Size for the computed part of susmat 260 bdtot_index=0 ; icg=0 ; ibg=0 261 ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3) 262 263 !ndiel4,ndiel5,ndiel6 are FFT dimensions, modified to avoid cache trashing 264 ndiel4=ngfftdiel(4) ; ndiel5=ngfftdiel(5) ; ndiel6=ngfftdiel(6) 265 diegap=dielar(5) ; dielam=dielar(6) 266 extrap=0 267 268 !If dielam is too small, there is no extrapolation. 269 if(dielam>1.0d-6)extrap=1 270 271 !Some stuff for symmetries 272 nsym1=sum(symafm,mask=symafm==1) 273 nsym2=nsym-nsym1 274 invnsym =one/dble(nsym) 275 invnsym1=one/dble(nsym1) 276 invnsym2=one 277 !FIXME: make sure this is consistent with following code 278 !div by 0 for several v5 tests 279 if (nsym2 > 0) invnsym2=one/dble(nsym2) 280 281 !Allocations 282 ABI_MALLOC(occ_deavg,(mband)) 283 if(occopt>=3) then 284 ABI_MALLOC(drhode,(2,npwdiel,nspden_eff)) 285 else 286 ABI_MALLOC(drhode,(0,0,0)) 287 end if 288 if(extrap==1) then 289 ABI_MALLOC(rhoextrap,(ndiel4,ndiel5,ndiel6,nspinor)) 290 else 291 ABI_MALLOC(rhoextrap,(0,0,0,0)) 292 end if 293 294 !zero the susceptibility matrix and other needed quantities 295 susmat(:,:,:,:,:)=zero 296 if(occopt>=3)then 297 drhode(:,:,:)=zero 298 sumdocc=zero 299 end if 300 301 !PAW additional initializations 302 if (usepaw==1) then 303 ABI_MALLOC(gylmg_diel,(npwdiel,lmax_diel**2,ntypat)) 304 ABI_MALLOC(ph3d_diel,(2,npwdiel,natom)) 305 if (neglect_pawhat==0) then 306 ABI_MALLOC(phkxred_diel,(2,natom)) 307 ABI_MALLOC(kpg_dum,(0,0)) 308 kpt_diel(1:3,1)=zero;phkxred_diel(1,:)=one;phkxred_diel(2,:)=zero;nkpg_diel=0 309 ! write(std_out,*) ' lmax_diel ', lmax_diel 310 call pawgylmg(gprimd,gylmg_diel,kg_diel,kpg_dum,kpt_diel,lmax_diel,nkpg_diel,npwdiel,ntypat,pawtab,ylmdiel) 311 call ph1d3d(1,natom,kg_diel,natom,natom,npwdiel,ndiel1,ndiel2,ndiel3,phkxred_diel,ph1ddiel,ph3d_diel) 312 ABI_FREE(phkxred_diel) 313 ABI_FREE(kpg_dum) 314 else 315 gylmg_diel=zero;ph3d_diel=one 316 end if 317 else 318 ABI_MALLOC(gylmg_diel,(0,0,0)) 319 ABI_MALLOC(ph3d_diel,(0,0,0)) 320 end if 321 322 call timab(741,2,tsec) 323 324 325 326 !--BIG loop over spins ------------------------------------------------------------ 327 !--------------------------------------------------------------------------------- 328 329 do isp=1,nsppol 330 ikg=0 331 332 if(extrap==1)rhoextrap(:,:,:,:)=zero 333 334 ! --BIG loop over k-points -------------------------------------------------------- 335 ! --------------------------------------------------------------------------------- 336 337 do ikpt=1,nkpt 338 339 nband_k=nband(ikpt+(isp-1)*nkpt) 340 istwf_k=istwfk(ikpt) 341 npw_k=npwarr(ikpt) 342 343 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isp,mpi_enreg%me_kpt)) then 344 bdtot_index=bdtot_index+nband_k 345 cycle 346 end if 347 348 call timab(742,1,tsec) 349 350 351 ABI_MALLOC(gbound,(2*mgfftdiel+8,2)) 352 ABI_MALLOC(kg_k,(3,npw_k)) 353 354 if (usepaw==1) then 355 ABI_MALLOC(cprj_k,(natom,my_nspinor*nband_k)) 356 if (neglect_pawhat==0) then 357 call pawcprj_alloc(cprj_k,0,dimcprj) 358 if (mpi_enreg%nproc_band==1) then 359 call pawcprj_get(atindx1,cprj_k,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,& 360 & mband,mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,unpaw,& 361 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 362 else 363 nband_loc=nband_k/mpi_enreg%nproc_band 364 ABI_MALLOC(cprj_loc,(natom,my_nspinor*nband_loc)) 365 call pawcprj_alloc(cprj_loc,0,dimcprj) 366 call pawcprj_get(atindx1,cprj_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,& 367 & mband/mpi_enreg%nproc_band,mkmem,natom,nband_loc,nband_loc,my_nspinor,nsppol,unpaw,& 368 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 369 call pawcprj_mpi_allgather(cprj_loc,cprj_k,natom,my_nspinor*nband_loc,mpi_enreg%bandpp,& 370 & dimcprj,0,mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.true.) 371 call pawcprj_free(cprj_loc) 372 ABI_FREE(cprj_loc) 373 end if 374 else 375 !call pawcprj_nullify(cprj_k) 376 end if 377 else 378 ABI_MALLOC(cprj_k,(0,0)) 379 end if 380 381 kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg) 382 call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k) 383 384 if(extrap==1)then 385 ! Compute inverse of average dielectric gap for each band 386 ! and multiply by occupation factor 387 emax=maxval(eigen(1+bdtot_index:nband_k+bdtot_index)) 388 do iband=1,nband_k 389 occ_deavg(iband)= occ(iband+bdtot_index)*dielam & 390 & / ( emax-eigen(iband+bdtot_index) + diegap ) 391 end do 392 else 393 occ_deavg(:)=zero 394 end if 395 396 call timab(742,2,tsec) 397 call timab(743,1,tsec) 398 399 ! Compute the contribution of each k-point to susmat, rhoextrap, drhode and sumdocc. 400 if(mpi_enreg%paral_kgb==1)then !Only this version is in parallel 401 ! Use either the simpler implementation 402 ! Should provide a test !!! 403 call susk(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,& 404 & gbound_diel,gylmg_diel,icg,ikpt,& 405 & isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 406 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 407 & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,& 408 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 409 & susmat,typat,ucvol,usepaw,wtk) 410 else 411 ! Or the more sophisticated one, needed to save memory. 412 call suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,& 413 & gbound_diel,gylmg_diel,icg,ikpt,& 414 & isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 415 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 416 & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,& 417 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 418 & susmat,typat,ucvol,usepaw,wtk) 419 end if 420 421 call timab(743,2,tsec) 422 423 ABI_FREE(gbound) 424 ABI_FREE(kg_k) 425 426 bdtot_index=bdtot_index+nband_k 427 428 if (mkmem/=0) then 429 ibg=ibg+my_nspinor*nband_k 430 icg=icg+my_nspinor*npw_k*nband_k 431 ikg=ikg+npw_k 432 end if 433 if (usepaw==1) then 434 if (neglect_pawhat==0) then 435 call pawcprj_free(cprj_k) 436 end if 437 end if 438 ABI_FREE(cprj_k) 439 440 ! End loop on ikpt: -------------------------------------------------------- 441 end do 442 443 ! Here include the contribution from the extrapolation to susmat, 444 ! diagonal part 445 if(extrap==1)then 446 447 call timab(744,1,tsec) 448 449 ! Transfer extrapolating density on augmented fft grid to 450 ! normal fft grid in real space. 451 ! Warning1 : if collinear magnetism, must treat only one spin at a time 452 ! Warning2 : if non-collinear magnetism, treat both spins 453 ! Warning3 : this is subtle for antiferro magnetism 454 nspden_tmp=1;if (antiferro) nspden_tmp=2 455 ABI_MALLOC(rhoextrr,(nfftdiel,nspden_tmp)) 456 ABI_MALLOC(rhoextrg,(2,nfftdiel)) 457 if (nspden==1.and.nspinor==2) rhoextrap(:,:,:,1)=rhoextrap(:,:,:,1)+rhoextrap(:,:,:,2) 458 459 do ispinor=1,min(nspinor,nspden) 460 jsp=isp+ispinor-1 461 462 call fftpac(1,mpi_enreg_diel,1,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6,& 463 & ngfftdiel,rhoextrr(:,1),rhoextrap(:,:,:,ispinor),1) 464 465 ! Generate the density in reciprocal space, and symmetrize it 466 ! (note symrhg also make the reverse FFT, to get symmetrized density; 467 ! this is useless here, and should be made an option) 468 call symrhg(1,gprimd,irrzondiel,mpi_enreg_diel,nfftdiel,nfftdiel,ngfftdiel,& 469 & nspden_tmp,1,nsym,phnonsdiel,rhoextrg,rhoextrr,rprimd,symafm,symrel,tnons) 470 471 do ipw2=1,npwdiel 472 j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2) 473 ! static: Only fills lower half of the matrix (here, the susceptibility matrix) 474 ! dynamical: fill all, will not affect susopt==2 for which extrap==0 475 do ipw1=1,npwdiel 476 i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1) 477 ! NOTE that there is a FFT folding (superposition) bias here 478 ! Should use kgindex, in the same spirit as in prcref 479 k1=i1-j1; k1=modulo(k1,ndiel1) 480 k2=i2-j2; k2=modulo(k2,ndiel2) 481 k3=i3-j3; k3=modulo(k3,ndiel3) 482 ifft=k1+1+ndiel1*(k2+ndiel2*k3) 483 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+rhoextrg(1,ifft) 484 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+rhoextrg(2,ifft) 485 end do 486 end do 487 488 end do 489 ABI_FREE(rhoextrg) 490 ABI_FREE(rhoextrr) 491 492 call timab(744,2,tsec) 493 494 end if 495 496 ! End loop over spins --------------------------------------------------------- 497 end do 498 499 ABI_FREE(occ_deavg) 500 ABI_FREE(rhoextrap) 501 ABI_FREE(gylmg_diel) 502 ABI_FREE(ph3d_diel) 503 !end if 504 505 call destroy_mpi_enreg(mpi_enreg_diel) 506 507 !-- Stuff for parallelism -------------------------------------------------------- 508 !--------------------------------------------------------------------------------- 509 510 if(xmpi_paral==1)then 511 call timab(746,1,tsec) 512 ABI_MALLOC(sussum,(2*npwdiel*nspden*npwdiel*nspden)) 513 ! Recreate full susmat on all proc. 514 ! This should be coded more efficiently, 515 ! since half of the matrix is still empty, and 516 ! it is spin-diagonal. 517 sussum(:)=reshape(susmat(:,:,:,:,:),(/2*npwdiel*nspden*npwdiel*nspden/)) 518 call xmpi_sum(sussum,spaceComm,ierr) 519 susmat(:,:,:,:,:)=reshape(sussum(:),(/2,npwdiel,nspden,npwdiel,nspden/)) 520 ABI_FREE(sussum) 521 ! Recreate full drhode on all proc. 522 if(occopt>=3 .and. testocc==1)then 523 call xmpi_sum(drhode,spaceComm,ierr) 524 ! Should use only one mpi-allreduce call instead of the three 525 call xmpi_sum(sumdocc,spaceComm,ierr) 526 end if 527 call timab(746,2,tsec) 528 end if 529 530 !-- Apply spatial hermitian/symmetries on spin-diagonal susceptibility matrix ---- 531 !--------------------------------------------------------------------------------- 532 533 call timab(747,1,tsec) 534 535 !If antiferro magnetism, has to divide (spin-diagonal) susceptibility by 2 (due to dble occupations) 536 if (antiferro) then 537 do ipw2=1,npwdiel 538 do ipw1=ipw2,npwdiel 539 susmat(:,ipw1,1,ipw2,1)=half*susmat(:,ipw1,1,ipw2,1) 540 end do 541 end do 542 end if 543 544 !Generate upper half of the spin-diagonal matrix (still the susceptibility matrix) 545 do isp=1,nspden_eff 546 do ipw2=2,npwdiel 547 do ipw1=1,ipw2-1 548 susmat(1,ipw1,isp,ipw2,isp)= susmat(1,ipw2,isp,ipw1,isp) 549 susmat(2,ipw1,isp,ipw2,isp)=-susmat(2,ipw2,isp,ipw1,isp) 550 end do 551 end do 552 end do 553 554 !Compute symmetric of G-vectors and eventual phases 555 !(either time-reversal or spatial symmetries) 556 ABI_MALLOC(tmrev_g,(npwdiel)) 557 ABI_MALLOC(sym_g,(npwdiel,nsym)) 558 ABI_MALLOC(phdiel,(2,npwdiel,nsym)) 559 call symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons) 560 561 !Impose spatial symmetries to the spin-diagonal susceptibility matrix 562 ABI_MALLOC(suswk,(2,npwdiel,npwdiel,nspden_eff)) 563 do isp=1,nspden_eff 564 suswk(:,:,:,isp)=susmat(:,:,isp,:,isp) ! Temporary storage 565 end do 566 567 do isp=1,nspden_eff 568 jsp=min(3-isp,nsppol) 569 do ipw2=1,npwdiel 570 do ipw1=1,npwdiel 571 ar=suswk(1,ipw1,ipw2,isp) 572 ai=suswk(2,ipw1,ipw2,isp) 573 ar2=zero;ai2=zero 574 if(nsym>1)then 575 do isym=2,nsym 576 t1=sym_g(ipw1,isym) ; t2=sym_g(ipw2,isym) 577 ! Not all symmetries are non-symmorphic. Should save time here ... 578 phr1=phdiel(1,ipw1,isym) ; phi1=phdiel(2,ipw1,isym) 579 phr2=phdiel(1,ipw2,isym) ; phi2=phdiel(2,ipw2,isym) 580 phr12= phr1*phr2+phi1*phi2 ; phi12=phi1*phr2-phr1*phi2 581 if (symafm(isym)==1) then 582 ar=ar+suswk(1,t1,t2,isp)*phr12-suswk(2,t1,t2,isp)*phi12 583 ai=ai+suswk(2,t1,t2,isp)*phr12+suswk(1,t1,t2,isp)*phi12 584 else 585 ar2=ar2+suswk(1,t1,t2,jsp)*phr12-suswk(2,t1,t2,jsp)*phi12 586 ai2=ai2+suswk(2,t1,t2,jsp)*phr12+suswk(1,t1,t2,jsp)*phi12 587 end if 588 end do 589 end if 590 if (antiferro) then 591 susmat(1,ipw1,1,ipw2,1)=ar*invnsym1 592 susmat(2,ipw1,1,ipw2,1)=ai*invnsym1 593 susmat(1,ipw1,2,ipw2,2)=ar2*invnsym2 594 susmat(2,ipw1,2,ipw2,2)=ai2*invnsym2 595 else 596 susmat(1,ipw1,isp,ipw2,isp)=(ar+ar2)*invnsym 597 susmat(2,ipw1,isp,ipw2,isp)=(ai+ai2)*invnsym 598 end if 599 end do 600 end do 601 end do 602 ABI_FREE(suswk) 603 604 605 !-- Add contribibution to susceptibility due to change of Fermi level ----------- 606 !--------------------------------------------------------------------------------- 607 608 if (occopt>=3.and.testocc==1) then 609 610 ! Impose spatial symmetries to drhode 611 ABI_MALLOC(drhode_wk,(2,npwdiel,nspden_eff)) 612 do isp=1,nspden_eff 613 jsp=min(3-isp,nsppol) 614 do ipw1=1,npwdiel 615 ar=drhode(1,ipw1,isp) 616 ai=drhode(2,ipw1,isp) 617 ar2=zero;ai2=zero 618 if (nsym>1) then 619 do isym=2,nsym 620 t1=sym_g(ipw1,isym) 621 ! Not all symmetries are non-symmorphic. Should save time here ... 622 phr1=phdiel(1,ipw1,isym);phi1=phdiel(2,ipw1,isym) 623 if (symafm(isym)==1) then 624 ar=ar+drhode(1,t1,isp)*phr1-drhode(2,t1,isp)*phi1 625 ai=ai+drhode(2,t1,isp)*phr1+drhode(1,t1,isp)*phi1 626 else 627 ar2=ar2+drhode(1,t1,jsp)*phr1-drhode(2,t1,jsp)*phi1 628 ai2=ai2+drhode(2,t1,jsp)*phr1+drhode(1,t1,jsp)*phi1 629 end if 630 end do 631 end if 632 if (antiferro) then ! 1/2 factor due to (dble) occupations 633 drhode_wk(1,ipw1,1)=half*ar*invnsym1 634 drhode_wk(2,ipw1,1)=half*ai*invnsym1 635 drhode_wk(1,ipw1,2)=half*ar2*invnsym2 636 drhode_wk(2,ipw1,2)=half*ai2*invnsym2 637 else 638 drhode_wk(1,ipw1,isp)=(ar+ar2)*invnsym 639 drhode_wk(2,ipw1,isp)=(ai+ai2)*invnsym 640 end if 641 end do 642 end do 643 644 ! Add contribution to non-diagonal susceptibility 645 ! Presently fills complete susceptibility matrix, not only lower half 646 weight=one/sumdocc 647 do isp2=1,nspden_eff 648 do ipw2=1,npwdiel 649 do isp1=1,nspden_eff 650 do ipw1=1,npwdiel 651 susmat(1,ipw1,isp1,ipw2,isp2)=susmat(1,ipw1,isp1,ipw2,isp2)- & 652 & weight*( drhode_wk(1,ipw1,isp1)*drhode_wk(1,ipw2,isp2) & 653 & +drhode_wk(2,ipw1,isp1)*drhode_wk(2,ipw2,isp2) ) 654 susmat(2,ipw1,isp1,ipw2,isp2)=susmat(2,ipw1,isp1,ipw2,isp2)- & 655 & weight*( drhode_wk(2,ipw1,isp1)*drhode_wk(1,ipw2,isp2) & 656 & -drhode_wk(1,ipw1,isp1)*drhode_wk(2,ipw2,isp2) ) 657 end do 658 end do 659 end do 660 end do 661 ABI_FREE(drhode_wk) 662 663 end if 664 !if (occopt>=3) then 665 ABI_FREE(drhode) 666 !end if 667 668 669 !--- Impose the time-reversal symmetry to the susceptibility matrix -------------- 670 !--------------------------------------------------------------------------------- 671 672 if (usetimerev==1) then 673 ABI_MALLOC(suswk,(2,npwdiel,npwdiel,1)) 674 675 ! Impose the time-reversal symmetry to the spin-diagonal susceptibility matrix 676 do isp=1,nspden_eff 677 suswk(:,:,:,1)=susmat(:,:,isp,:,isp) ! Temporary storage 678 do ipw2=1,npwdiel 679 t2=tmrev_g(ipw2) 680 do ipw1=1,npwdiel 681 t1=tmrev_g(ipw1) 682 susmat(1,ipw1,isp,ipw2,isp)=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1)) 683 susmat(2,ipw1,isp,ipw2,isp)=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1)) 684 end do 685 end do 686 end do 687 688 ! Impose the time-reversal symmetry to the off-diagonal susceptibility matrix 689 if (nspden_eff/=1.and.occopt>=3.and.testocc==1) then 690 suswk(:,:,:,1)=susmat(:,:,1,:,2) ! Temporary storage 691 do ipw2=1,npwdiel 692 t2=tmrev_g(ipw2) 693 do ipw1=1,npwdiel 694 t1=tmrev_g(ipw1) 695 ar=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1)) 696 ai=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1)) 697 susmat(1,ipw1,1,ipw2,2)= ar 698 susmat(2,ipw1,1,ipw2,2)= ai 699 susmat(1,ipw1,2,ipw2,1)= ar 700 susmat(2,ipw1,2,ipw2,1)=-ai 701 end do 702 end do 703 end if 704 ABI_FREE(suswk) 705 end if 706 707 ABI_FREE(phdiel) 708 ABI_FREE(sym_g) 709 ABI_FREE(tmrev_g) 710 711 712 !-- The full susceptibility matrix is computed ----------------------------------- 713 !-- Now, eventually diagonalize it and stop -------------------------------------- 714 !--------------------------------------------------------------------------------- 715 716 !Must turn on this flag to make the diagonalisation 717 diag=0 718 if(diag==1)then 719 720 npwsp=npwdiel*nspden_eff 721 ABI_MALLOC(sush,(npwsp*(npwsp+1))) 722 ABI_MALLOC(susvec,(2,npwsp,npwsp)) 723 ABI_MALLOC(eig_diel,(npwsp)) 724 ABI_MALLOC(zhpev1,(2,2*npwsp-1)) 725 ABI_MALLOC(zhpev2,(3*npwsp-2)) 726 ier=0 727 728 ! Store the susceptibility matrix in proper mode before calling zhpev 729 indx=1 730 do ii=1,npwdiel 731 do jj=1,ii 732 sush(indx )=susmat(1,jj,1,ii,1) 733 sush(indx+1)=susmat(2,jj,1,ii,1) 734 indx=indx+2 735 end do 736 end do 737 738 ! If spin-polarized, need to store other parts of the matrix 739 if(nspden_eff/=1)then 740 do ii=1,npwdiel 741 ! Here, spin-flip contribution 742 do jj=1,npwdiel 743 sush(indx )=susmat(1,jj,1,ii,2) 744 sush(indx+1)=susmat(2,jj,1,ii,2) 745 indx=indx+2 746 end do 747 ! Here spin down-spin down upper matrix 748 do jj=1,ii 749 sush(indx )=susmat(1,jj,2,ii,2) 750 sush(indx+1)=susmat(2,jj,2,ii,2) 751 indx=indx+2 752 end do 753 end do 754 end if 755 756 call ZHPEV ('V','U',npwsp,sush,eig_diel,susvec,npwsp,zhpev1,& 757 & zhpev2,ier) 758 759 write(std_out,*)' suscep_stat : print eigenvalues of the susceptibility matrix' 760 do ii=1,npwsp 761 write(std_out,'(i5,es16.6)' )ii,eig_diel(ii) 762 end do 763 764 ABI_FREE(sush) 765 ABI_FREE(susvec) 766 ABI_FREE(eig_diel) 767 ABI_FREE(zhpev1) 768 ABI_FREE(zhpev2) 769 ABI_ERROR("Stopping here!") 770 end if 771 772 call timab(747,2,tsec) 773 call timab(740,2,tsec) 774 775 end subroutine suscep_stat
m_suscep_stat/susk [ Functions ]
[ Top ] [ m_suscep_stat ] [ 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.
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/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
SOURCE
877 subroutine susk(atindx,bdtot_index,cg_mpi,cprj_k,doccde,drhode,eigen,extrap,gbound,& 878 & gbound_diel,gylmg_diel,icg_mpi,ikpt,isp,istwf_k,kg_diel,kg_k_mpi,& 879 & lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 880 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 881 & npwdiel,npw_k_mpi,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,& 882 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 883 & susmat,typat,ucvol,usepaw,wtk) 884 885 !Arguments ------------------------------------ 886 !This type is defined in defs_mpi 887 !scalars 888 integer,intent(in) :: bdtot_index,extrap,ikpt,isp,istwf_k,lmax_diel,mband,mcg 889 integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat 890 integer,intent(in) :: nkpt,npwdiel,nspden,nspden_eff,nspinor,nsppol 891 integer,intent(in) :: ntypat,occopt,usepaw 892 integer,intent(in),target :: icg_mpi,npw_k_mpi 893 real(dp),intent(in) :: ucvol 894 real(dp),intent(inout) :: sumdocc 895 type(MPI_type),intent(in) :: mpi_enreg 896 type(pawang_type),intent(in) :: pawang 897 !arrays 898 integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2) 899 integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom) 900 integer,intent(in),target :: kg_k_mpi(3,npw_k_mpi) 901 integer,intent(inout) :: gbound(2*mgfftdiel+8,2) 902 integer,pointer :: kg_k(:,:) 903 real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol) 904 real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw) 905 real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband) 906 real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt) 907 real(dp),intent(in),target :: cg_mpi(2,mcg) 908 real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff) 909 real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor) 910 real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 911 type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw) 912 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 913 914 !Local variables------------------------------- 915 ! real(dp), allocatable :: cg_disk(:,:) 916 !Local variables for MPI 917 !scalars 918 integer :: blocksize,i1,i2,i3,iband,iband_loc,ibd1,ibd2,ibdblock,ier 919 integer :: iproc,iproc_fft,ipw,ipw1,ipw2,isp1,isp2,ispinor,iwf,jsp,me_bandfft 920 integer :: nbdblock,ndatarecv,ndiel1,ndiel2,ndiel3 921 integer :: sizemax_per_proc,spaceComm,testocc,tim_fourwf 922 integer,pointer :: icg,npw_k 923 integer,target :: icg_loc=0,npw_k_loc,npw_tot 924 real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2 925 type(MPI_type) :: mpi_enreg_diel 926 !arrays 927 integer,allocatable :: band_loc(:),kg_k_gather(:,:),npw_per_proc(:),rdispls(:) 928 integer,allocatable :: rdispls_all(:),rdisplsloc(:),recvcounts(:) 929 integer,allocatable :: recvcountsloc(:),sdispls(:),sdisplsloc(:),sendcounts(:) 930 integer,allocatable :: sendcountsloc(:) 931 integer,allocatable,target :: kg_k_gather_all(:,:) 932 real(dp) :: tsec(2) 933 real(dp),allocatable :: cwavef(:,:),cwavef_alltoall(:,:) 934 real(dp),allocatable :: cwavef_alltoall_gather(:,:),dummy(:,:),rhoaug(:,:,:) 935 real(dp),allocatable :: susmat_mpi(:,:,:) 936 real(dp),allocatable :: wfprod(:,:),wfraug(:,:,:,:),wfrspa(:,:,:,:,:,:) 937 real(dp),allocatable,target :: cg_local(:,:) 938 real(dp),pointer :: cg(:,:) 939 logical,allocatable :: treat_band(:) 940 941 ! ************************************************************************* 942 943 !DEBUG 944 !write(std_out,*)' susk : enter ' 945 !if(.true.)stop 946 !ENDDEBUG 947 948 call timab(750,1,tsec) 949 call timab(751,1,tsec) 950 951 ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3) 952 953 !The dielectric stuff is performed in sequential mode. 954 !Set mpi_enreg_diel accordingly 955 call initmpi_seq(mpi_enreg_diel) 956 call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 957 me_bandfft=xmpi_comm_rank(mpi_enreg%comm_bandfft) 958 959 testocc=1 960 !DEBUG 961 !write(std_out,*)' susk : set testocc to 0 ' 962 !testocc=0 963 !write(std_out,*)' susk : set extrap to 0 ' 964 !extrap=0 965 !ENDDEBUG 966 967 !Allocations, initializations 968 ABI_MALLOC(rhoaug,(ndiel4,ndiel5,ndiel6)) 969 ABI_MALLOC(wfraug,(2,ndiel4,ndiel5,ndiel6)) 970 ABI_MALLOC(wfprod,(2,npwdiel)) 971 ABI_MALLOC(wfrspa,(2,ndiel4,ndiel5,ndiel6,nspinor,mband)) 972 ABI_MALLOC(dummy,(2,1)) 973 wfrspa(:,:,:,:,:,:)=zero 974 ABI_MALLOC(treat_band,(nband_k)) 975 treat_band(:)=.true. 976 isp1=isp;isp2=isp 977 if (nspden_eff==2.and.nspinor==2) isp2=isp+1 978 979 !BAND-FFT parallelism 980 if (mpi_enreg%paral_kgb==1) then 981 treat_band(:)=.false. 982 ! We gather the wavefunctions treated by this proc in cg_local 983 spaceComm=mpi_enreg%comm_band 984 blocksize=mpi_enreg%nproc_band 985 nbdblock=nband_k/blocksize 986 ABI_MALLOC(sdispls,(blocksize)) 987 ABI_MALLOC(sdisplsloc,(blocksize)) 988 ABI_MALLOC(sendcounts,(blocksize)) 989 ABI_MALLOC(sendcountsloc,(blocksize)) 990 ABI_MALLOC(rdispls,(blocksize)) 991 ABI_MALLOC(rdisplsloc,(blocksize)) 992 ABI_MALLOC(recvcounts,(blocksize)) 993 ABI_MALLOC(recvcountsloc,(blocksize)) 994 ! First gather the kg_k in kg_k_gather_all 995 npw_k_loc=npw_k_mpi 996 call xmpi_allgather(npw_k_loc,recvcounts,spaceComm,ier) 997 rdispls(1)=0 998 do iproc=2,blocksize 999 rdispls(iproc)=rdispls(iproc-1)+recvcounts(iproc-1) 1000 end do 1001 ndatarecv=rdispls(blocksize)+recvcounts(blocksize) 1002 ABI_MALLOC(kg_k_gather,(3,ndatarecv)) 1003 recvcountsloc(:)=recvcounts(:)*3 1004 rdisplsloc(:)=rdispls(:)*3 1005 call xmpi_allgatherv(kg_k_mpi,3*npw_k_loc,kg_k_gather,recvcountsloc(:),rdisplsloc,spaceComm,ier) 1006 ABI_MALLOC(npw_per_proc,(mpi_enreg%nproc_fft)) 1007 ABI_MALLOC(rdispls_all,(mpi_enreg%nproc_fft)) 1008 spaceComm=mpi_enreg%comm_fft 1009 call xmpi_allgather(ndatarecv,npw_per_proc,spaceComm,ier) 1010 rdispls_all(1)=0 1011 do iproc=2,mpi_enreg%nproc_fft 1012 rdispls_all(iproc)=rdispls_all(iproc-1)+npw_per_proc(iproc-1) 1013 end do 1014 npw_tot=rdispls_all(mpi_enreg%nproc_fft)+npw_per_proc(mpi_enreg%nproc_fft) 1015 ABI_MALLOC(kg_k_gather_all,(3,npw_tot)) 1016 call xmpi_allgatherv(kg_k_gather,3*ndatarecv,kg_k_gather_all,3*npw_per_proc(:),3*rdispls_all,spaceComm,ier) 1017 ! At this point kg_k_gather_all contains all the kg 1018 if(allocated(cwavef)) then 1019 ABI_FREE(cwavef) 1020 end if 1021 ABI_MALLOC(cwavef,(2,npw_k_loc*nspinor*blocksize)) 1022 sizemax_per_proc=nband_k/(mpi_enreg%nproc_band*mpi_enreg%nproc_fft)+1 1023 ABI_MALLOC(band_loc,(sizemax_per_proc)) 1024 ABI_MALLOC(cg_local,(2,sizemax_per_proc*npw_tot*nspinor)) 1025 iband_loc=0 1026 do ibdblock=1,nbdblock 1027 cwavef(:,1:npw_k_loc*nspinor*blocksize)=& 1028 & cg_mpi(:,1+(ibdblock-1)*npw_k_loc*nspinor*blocksize+icg_mpi:ibdblock*npw_k_loc*nspinor*blocksize+icg_mpi) 1029 sendcounts(:)=npw_k_loc 1030 do iproc=1,blocksize 1031 sdispls(iproc)=(iproc-1)*npw_k_loc 1032 end do 1033 ABI_MALLOC(cwavef_alltoall,(2,ndatarecv*nspinor)) 1034 recvcountsloc(:)=recvcounts(:)*2*nspinor 1035 rdisplsloc(:)=rdispls(:)*2*nspinor 1036 sendcountsloc(:)=sendcounts(:)*2*nspinor 1037 sdisplsloc(:)=sdispls(:)*2*nspinor 1038 call timab(547,1,tsec) 1039 spaceComm=mpi_enreg%comm_band 1040 call xmpi_alltoallv(cwavef,sendcountsloc,sdisplsloc,cwavef_alltoall,recvcountsloc,rdisplsloc,spaceComm,ier) 1041 call timab(547,2,tsec) 1042 ABI_MALLOC(cwavef_alltoall_gather,(2,npw_tot*nspinor)) 1043 blocksize=mpi_enreg%nproc_band 1044 spaceComm=mpi_enreg%comm_fft 1045 call xmpi_allgatherv(cwavef_alltoall,2*nspinor*ndatarecv,cwavef_alltoall_gather,& 1046 & 2*nspinor*npw_per_proc,2*nspinor*rdispls_all,spaceComm,ier) 1047 iproc_fft=modulo(ibdblock-1,mpi_enreg%nproc_fft) 1048 if(mpi_enreg%me_fft==iproc_fft) then !All nproc_band procs of index me_fft will treat these bands 1049 iband_loc=iband_loc+1 1050 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 1051 treat_band(iband)=.true. 1052 band_loc(iband_loc)=iband 1053 cg_local(:,1+(iband_loc-1)*npw_tot*nspinor:iband_loc*npw_tot*nspinor)=cwavef_alltoall_gather(:,1:npw_tot*nspinor) 1054 end if 1055 ABI_FREE(cwavef_alltoall_gather) 1056 ABI_FREE(cwavef_alltoall) 1057 end do 1058 ! On exit: 1059 ! npw_tot will be npw 1060 ! kg_k_gather_all will be kg_k 1061 ! cg_local will be cg 1062 ! icg will be zero 1063 npw_k=>npw_tot 1064 kg_k=>kg_k_gather_all(:,:) 1065 cg=>cg_local(:,:) 1066 icg=>icg_loc 1067 call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k) 1068 ABI_FREE(npw_per_proc) 1069 ABI_FREE(rdispls_all) 1070 ABI_FREE(sendcounts) 1071 ABI_FREE(recvcounts) 1072 ABI_FREE(sdispls) 1073 ABI_FREE(rdispls) 1074 ABI_FREE(sendcountsloc) 1075 ABI_FREE(sdisplsloc) 1076 ABI_FREE(recvcountsloc) 1077 ABI_FREE(rdisplsloc) 1078 ABI_FREE(kg_k_gather) 1079 ABI_FREE(cwavef) 1080 ! Because they will be summed over all procs, and arrive on input, rescale drhode and rhoextrap 1081 if(occopt>=3)drhode(:,:,isp1:isp2)=drhode(:,:,isp1:isp2)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 1082 if(extrap==1)rhoextrap(:,:,:,:)=rhoextrap(:,:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 1083 do i1=isp1,isp2 1084 susmat(:,:,i1,:,i1)=susmat(:,:,i1,:,i1)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 1085 end do 1086 1087 ! No BAND-FFT parallelism 1088 else ! use argument variables 1089 cg=>cg_mpi 1090 kg_k=>kg_k_mpi 1091 npw_k=>npw_k_mpi 1092 icg=>icg_mpi 1093 end if 1094 iband_loc=0 1095 1096 call timab(751,2,tsec) 1097 call timab(752,1,tsec) 1098 1099 !Loop over bands to fft and store Fourier transform of wavefunction 1100 ABI_MALLOC(cwavef,(2,npw_k)) 1101 do iband=1,nband_k 1102 if(.not. treat_band(iband)) cycle ! I am not treating this band (only for the parallel case) 1103 iband_loc=iband_loc+1 1104 1105 ! Loop on spinorial components 1106 do ispinor=1,nspinor 1107 iwf=(ispinor-1)*npw_k+(iband_loc-1)*npw_k*nspinor+icg 1108 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 1109 1110 ! Obtain Fourier transform in fft box 1111 tim_fourwf=8 1112 cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf) 1113 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 1114 & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,& 1115 & 0,tim_fourwf,weight,weight) 1116 1117 wfrspa(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:) 1118 1119 if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then 1120 ! In the case of metallic occupation, or if the extrapolation 1121 ! over higher bands is included, must compute the 1122 ! Fourier transform of the density of each band, then 1123 ! generate the part of the susceptibility matrix due 1124 ! varying occupation numbers. 1125 1126 weight=-two*occ_deavg(iband)*wtk(ikpt)/ucvol 1127 do i3=1,ndiel3 1128 do i2=1,ndiel2 1129 do i1=1,ndiel1 1130 wfraug(1,i1,i2,i3)=wfraug(1,i1,i2,i3)**2+wfraug(2,i1,i2,i3)**2 1131 wfraug(2,i1,i2,i3)=zero 1132 end do 1133 end do 1134 ! If extrapolation, accumulate density in real space 1135 if(extrap==1.and.usepaw==0)then 1136 do i2=1,ndiel2 1137 do i1=1,ndiel1 1138 rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3) 1139 end do 1140 end do 1141 end if 1142 end do 1143 1144 ! In case of PAW, add compensation charge contribution 1145 if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then 1146 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,& 1147 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1148 & ngfftdiel,npwdiel,nspinor,ntypat,1,& 1149 & pawang,pawtab,ph3d_diel,typat,dummy,wfraug,& 1150 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1151 rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:) 1152 end if 1153 1154 ! Performs the Fourier Transform of the density of the band, 1155 ! and store it in wfprod 1156 tim_fourwf=9 1157 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 1158 & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 1159 & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight) 1160 ! In case of PAW, add compensation charge contribution if not already done 1161 if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then 1162 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,& 1163 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1164 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 1165 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 1166 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1167 end if 1168 1169 ! Perform now the summation of terms related to direct change of eigenvalues 1170 ! or extrapolation over higher bands 1171 wght1=zero ; wght2=zero 1172 if(occopt>=3 .and. testocc==1) wght1=doccde(iband+bdtot_index)*wtk(ikpt)/ucvol 1173 if(extrap==1) wght2=two*occ_deavg(iband)*wtk(ikpt)/ucvol 1174 weight=wght1+wght2 1175 1176 if (abs(weight)>tol12) then 1177 do ipw2=1,npwdiel 1178 ! Only fills lower half of the matrix (here, the susceptibility matrix) 1179 ! Note that wfprod of the first index must behave like a density, 1180 ! so that it is used as generated by fourwf, while wfprod of the 1181 ! second index will be implicitely used to make a scalar product 1182 ! with a potential change, meaning that its complex conjugate must be 1183 ! used. This explains the following signs... 1184 do ipw1=ipw2,npwdiel 1185 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 1186 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 1187 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 1188 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 1189 end do 1190 end do 1191 end if 1192 1193 if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then 1194 ! Accumulate product of band densities by their doccde, for the 1195 ! computation of the effect of change of Fermi level. 1196 do ipw=1,npwdiel 1197 drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1 1198 drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1 1199 end do 1200 ! Also accumulate weighted sum of doccde 1201 sumdocc=sumdocc+wght1 1202 end if 1203 1204 ! End condition of metallic occupancies or extrapolation 1205 end if 1206 1207 ! End loop on spinorial components 1208 end do 1209 ! End loop on iband 1210 end do 1211 1212 call timab(752,2,tsec) 1213 call timab(753,1,tsec) 1214 1215 ABI_FREE(cwavef) 1216 1217 !Stuff for parallelism (bands-FFT) 1218 if(mpi_enreg%paral_kgb==1) then 1219 call xmpi_sum(wfrspa,mpi_enreg%comm_bandfft,ier) 1220 if(occopt>=3) then 1221 call xmpi_sum(drhode(:,:,isp1:isp2),mpi_enreg%comm_bandfft,ier) 1222 end if 1223 if(extrap==1) then 1224 call xmpi_sum(rhoextrap,mpi_enreg%comm_bandfft,ier) 1225 end if 1226 if(occopt>=3) then 1227 call xmpi_sum(sumdocc,mpi_enreg%comm_bandfft,ier) 1228 end if 1229 ABI_MALLOC(susmat_mpi,(2,npwdiel,npwdiel)) 1230 do i1=isp1,isp2 1231 susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1) 1232 call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier) 1233 susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp) 1234 end do 1235 ABI_FREE(susmat_mpi) 1236 end if 1237 call timab(753,2,tsec) 1238 1239 !-- Wavefunctions have been generated in real space ------------------------ 1240 !-- Now, compute product of wavefunctions for different bands -------------- 1241 call timab(754,1,tsec) 1242 !if (occopt<3) then 1243 tolocc=1.0d-3 1244 !else 1245 !tolocc=1.0d-8 1246 !end if 1247 iproc=-1 1248 1249 if(nband_k>1)then 1250 do ibd1=1,nband_k-1 1251 do ibd2=ibd1+1,nband_k 1252 iproc=iproc+1 1253 if(modulo(iproc,mpi_enreg%nproc_fft*mpi_enreg%nproc_band) /= me_bandfft) cycle 1254 ! If the occupation numbers are sufficiently different, or 1255 ! if extrapolation is used and the corresponding factor is not zero, 1256 ! then there is a contribution 1257 occdiff=occ(ibd1+bdtot_index)-occ(ibd2+bdtot_index) 1258 if( abs(occdiff)>tolocc .or. & 1259 & ( extrap==1 .and. & 1260 & ( abs(occ_deavg(ibd1)) + abs(occ_deavg(ibd2)) ) >tolocc ) & 1261 & ) then 1262 1263 eigdiff=eigen(ibd1+bdtot_index)-eigen(ibd2+bdtot_index) 1264 ! DEBUG 1265 ! write(std_out,*)' susk : contribution from bands',ibd1,ibd2 1266 ! write(std_out,*)' occ diff =',occdiff 1267 ! write(std_out,*)' eig diff =',eigdiff 1268 ! ENDDEBUG 1269 1270 ! Loop on spinorial components 1271 do ispinor=1,nspinor 1272 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 1273 1274 ! Store the contribution in wfraug 1275 do i3=1,ndiel3 1276 do i2=1,ndiel2 1277 do i1=1,ndiel1 1278 wfraug(1,i1,i2,i3)=wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)& 1279 & +wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2) 1280 wfraug(2,i1,i2,i3)=wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)& 1281 & -wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2) 1282 end do 1283 end do 1284 end do 1285 1286 ! Performs the Fourier Transform of the product, and store it in wfprod 1287 tim_fourwf=19 1288 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 1289 & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 1290 & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight) 1291 1292 ! In case of PAW, add compensation charge contribution 1293 if (usepaw==1.and.neglect_pawhat==0) then 1294 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,& 1295 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1296 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 1297 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 1298 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1299 end if 1300 1301 ! Perform now the summation 1302 wght1=zero ; wght2=zero 1303 if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol 1304 if(extrap==1) wght2=(occ_deavg(ibd1)+occ_deavg(ibd2)) * two*wtk(ikpt)/ucvol 1305 weight=wght1+wght2 1306 1307 ! DEBUG 1308 ! write(std_out,*)' weight =',weight 1309 ! norm=zero 1310 ! do ipw=1,npwdiel 1311 ! norm=norm+wfprod(1,ipw)**2+wfprod(2,ipw)**2 1312 ! end do 1313 ! write(std_out,*)' norm in reciprocal space =',norm 1314 ! ENDDEBUG 1315 1316 if (abs(weight)>tol12) then 1317 do ipw2=1,npwdiel 1318 ! Only fills lower half of the matrix (here, the susceptibility matrix) 1319 ! Note that wfprod of the first index must behave like a density, 1320 ! so that it is used as generated by fourwf, while wfprod of the 1321 ! second index will be implicitely used to make a scalar product 1322 ! with a potential change, meaning that its complex conjugate must be 1323 ! used. This explains the following signs... 1324 do ipw1=ipw2,npwdiel 1325 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 1326 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 1327 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 1328 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 1329 end do 1330 end do 1331 end if 1332 1333 ! End loop on spinorial components 1334 end do 1335 ! End condition of different occupation numbers or extrapolation 1336 end if 1337 ! End internal loop over bands 1338 end do 1339 ! End external loop over bands 1340 end do 1341 ! End condition of having more than one band 1342 end if 1343 1344 call timab(754,2,tsec) 1345 call timab(755,1,tsec) 1346 1347 if(mpi_enreg%paral_kgb==1) then 1348 ABI_MALLOC(susmat_mpi,(2,npwdiel,npwdiel)) 1349 do i1=isp1,isp2 1350 susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1) 1351 call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier) 1352 susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:) 1353 end do 1354 ABI_FREE(susmat_mpi) 1355 ABI_FREE(band_loc) 1356 ABI_FREE(treat_band) 1357 ABI_FREE(cg_local) 1358 ABI_FREE(kg_k_gather_all) 1359 end if 1360 1361 call destroy_mpi_enreg(mpi_enreg_diel) 1362 ABI_FREE(dummy) 1363 ABI_FREE(rhoaug) 1364 ABI_FREE(wfprod) 1365 ABI_FREE(wfraug) 1366 ABI_FREE(wfrspa) 1367 1368 call timab(755,2,tsec) 1369 call timab(750,2,tsec) 1370 1371 end subroutine susk
m_suscep_stat/suskmm [ Functions ]
[ Top ] [ m_suscep_stat ] [ Functions ]
NAME
suskmm
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). This routine is similar to susk, but use blocking on wavefunctions to decrease memory requirements, at the expense of CPU time.
NOTES
There is still room for optimization !!
INPUTS
atindx(natom)=index table for atoms bdtot_index=index for the number of the band cg(2,mcg)=wf 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**2,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)=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/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
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
SOURCE
1463 subroutine suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,& 1464 & gbound_diel,gylmg_diel,icg,ikpt,isp,istwf_k,kg_diel,kg_k,& 1465 & lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 1466 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 1467 & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,& 1468 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 1469 & susmat,typat,ucvol,usepaw,wtk) 1470 1471 !Arguments ------------------------------------ 1472 !scalars 1473 integer,intent(in) :: bdtot_index,extrap,icg,ikpt,isp,istwf_k,lmax_diel,mband,mcg 1474 integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat 1475 integer,intent(in) :: nkpt,npw_k,npwdiel,nspden,nspden_eff,nspinor 1476 integer,intent(in) :: nsppol,ntypat,occopt,usepaw 1477 real(dp),intent(in) :: ucvol 1478 real(dp),intent(inout) :: sumdocc 1479 type(MPI_type),intent(in) :: mpi_enreg 1480 type(pawang_type),intent(in) :: pawang 1481 !arrays 1482 integer,intent(in) :: atindx(natom),gbound(2*mgfftdiel+8,2) 1483 integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2) 1484 integer,intent(in) :: kg_diel(3,npwdiel),kg_k(3,npw_k),ngfftdiel(18) 1485 integer,intent(in) :: typat(natom) 1486 real(dp),intent(in) :: cg(2,mcg),doccde(mband*nkpt*nsppol) 1487 real(dp),intent(in) :: eigen(mband*nkpt*nsppol) 1488 real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw) 1489 real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband) 1490 real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt) 1491 real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff) 1492 real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor) 1493 real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 1494 type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw) 1495 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 1496 1497 !Local variables------------------------------- 1498 !scalars 1499 integer :: comm_fft,i1,i2,i3,iband,iband_shift,iband_shift2,ibd1,ibd2,ibdshft1,ibdshft2 1500 integer :: iblk1,iblk2,ipw,ipw1,ipw2,ispinor,iwf,jsp,mblk 1501 integer :: nblk,nbnd_current,nbnd_in_blk,nbnd_in_blk1,ndiel1,ndiel2,ndiel3 1502 integer :: testocc,tim_fourwf 1503 real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2 1504 character(len=500) :: message 1505 type(MPI_type) :: mpi_enreg_diel 1506 !arrays 1507 real(dp) :: tsec(2) 1508 real(dp),allocatable :: cwavef(:,:),dummy(:,:),rhoaug(:,:,:),wfprod(:,:) 1509 real(dp),allocatable :: wfraug(:,:,:,:),wfrspa1(:,:,:,:,:,:) 1510 real(dp),allocatable :: wfrspa2(:,:,:,:,:,:) 1511 1512 ! ************************************************************************* 1513 1514 call timab(760,1,tsec) 1515 call timab(761,1,tsec) 1516 1517 !Allocations, initializations 1518 ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3) 1519 testocc=1 1520 ABI_MALLOC(rhoaug,(ndiel4,ndiel5,ndiel6)) 1521 ABI_MALLOC(wfraug,(2,ndiel4,ndiel5,ndiel6)) 1522 ABI_MALLOC(wfprod,(2,npwdiel)) 1523 ABI_MALLOC(dummy,(2,1)) 1524 1525 !The dielectric stuff is performed in sequential mode. 1526 !Set mpi_enreg_diel accordingly 1527 call initmpi_seq(mpi_enreg_diel) 1528 call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 1529 1530 comm_fft=mpi_enreg%comm_fft 1531 1532 !Prepare the blocking : compute the number of blocks, 1533 !the number of bands in each normal block, 1534 !and the number in the first one, usually smaller. 1535 1536 !Consider that if the number of bands is large, there are at most 8 blocks 1537 nbnd_in_blk=0 1538 if(nband_k>=48)then 1539 mblk=8 1540 nbnd_in_blk=(nband_k-1)/mblk+1 1541 ! If the number of bands is medium, place 6 bands per block 1542 else if(nband_k>=12)then 1543 nbnd_in_blk=6 1544 ! Otherwise, must have at least 2 blocks 1545 else if(nband_k>=2)then 1546 mblk=2 1547 nbnd_in_blk=(nband_k-1)/mblk+1 1548 else 1549 write(message, '(a,a,a,i2,a,a,a)')& 1550 & ' The number of bands must be larger or equal to 2, in suskmm.',ch10,& 1551 & ' It is equal to ',nband_k,'.',ch10,& 1552 & ' Action : choose another preconditioner.' 1553 ABI_ERROR(message) 1554 end if 1555 1556 !Compute the effective number of blocks, and the number of bands in 1557 !the first block. 1558 nblk=(nband_k-1)/nbnd_in_blk+1 1559 nbnd_in_blk1=nband_k-(nblk-1)*nbnd_in_blk 1560 1561 !DEBUG 1562 !write(std_out,*)' suskmm : nband_k,nblk,nbnd_in_blk,nbnd_in_blk1 ' 1563 !write(std_out,*)nband_k,nblk,nbnd_in_blk,nbnd_in_blk1 1564 !stop 1565 !ENDDEBUG 1566 1567 !wfrspa1 will contain the wavefunctions of the slow sampling (iblk1) 1568 ABI_MALLOC(wfrspa1,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk)) 1569 !wfrspa2 will contain the wavefunctions of the rapid sampling (iblk2) 1570 ABI_MALLOC(wfrspa2,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk)) 1571 1572 ABI_MALLOC(cwavef,(2,npw_k)) 1573 1574 call timab(761,2,tsec) 1575 1576 !First loop over blocks 1577 do iblk1=1,nblk 1578 1579 call timab(762,1,tsec) 1580 1581 ! Initialisation 1582 if(iblk1==1)then 1583 1584 nbnd_current=nbnd_in_blk1 1585 iband_shift=0 1586 ! Loop over bands to fft and store Fourier transform of wavefunction 1587 do iband=1,nbnd_current 1588 ! Loop on spinorial components 1589 do ispinor=1,nspinor 1590 iwf=(ispinor-1)*npw_k+(iband-1)*npw_k*nspinor+icg 1591 ! Obtain Fourier transform in fft box 1592 tim_fourwf=21 1593 cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf) 1594 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 1595 & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,& 1596 & 0,tim_fourwf,weight,weight) 1597 wfrspa1(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:) 1598 end do 1599 end do 1600 1601 else 1602 1603 ! The Fourier transform of wavefunctions have already been obtained 1604 nbnd_current=nbnd_in_blk 1605 iband_shift=nbnd_in_blk1+(iblk1-2)*nbnd_in_blk 1606 1607 end if 1608 1609 ! Loop over bands of this block, to generate band-diagonal 1610 do iband=1,nbnd_current 1611 1612 ! Loop on spinorial components 1613 do ispinor=1,nspinor 1614 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 1615 1616 if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then 1617 ! In the case of metallic occupation, or if the extrapolation 1618 ! over higher bands is included, must compute the 1619 ! Fourier transform of the density of each band, then 1620 ! generate the part of the susceptibility matrix due 1621 ! varying occupation numbers. 1622 weight=-two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol 1623 do i3=1,ndiel3 1624 do i2=1,ndiel2 1625 do i1=1,ndiel1 1626 wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,iband)**2& 1627 & +wfrspa1(2,i1,i2,i3,ispinor,iband)**2 1628 wfraug(2,i1,i2,i3)=zero 1629 end do 1630 end do 1631 ! If extrapolation, accumulate density in real space 1632 if(extrap==1.and.usepaw==0)then 1633 do i2=1,ndiel2 1634 do i1=1,ndiel1 1635 rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3) 1636 end do 1637 end do 1638 end if 1639 end do 1640 1641 ! In case of PAW, add compensation charge contribution 1642 if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then 1643 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,& 1644 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1645 & ngfftdiel,npwdiel,nspinor,ntypat,1,& 1646 & pawang,pawtab,ph3d_diel,typat,dummy,wfraug,& 1647 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1648 rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:) 1649 end if 1650 1651 ! Performs the Fourier Transform of the density of the band, 1652 ! and store it in wfprod 1653 tim_fourwf=31 1654 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 1655 & 1,kg_diel,kg_diel,& 1656 & mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight) 1657 ! In case of PAW, add compensation charge contribution if not already done 1658 if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then 1659 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,1,1,1,kg_diel,& 1660 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1661 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 1662 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 1663 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1664 end if 1665 1666 ! Perform now the summation of terms related to direct change of eigenvalues 1667 ! or extrapolation over higher bands 1668 wght1=zero ; wght2=zero 1669 if(occopt>=3 .and. testocc==1) wght1=doccde(iband+iband_shift+bdtot_index)*wtk(ikpt)/ucvol 1670 if(extrap==1) wght2=two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol 1671 weight=wght1+wght2 1672 1673 if (abs(weight)>tol12) then 1674 do ipw2=1,npwdiel 1675 ! Only fills lower half of the matrix (here, the susceptibility matrix) 1676 ! Note that wfprod of the first index must behave like a density, 1677 ! so that it is used as generated by fourwf, while wfprod of the 1678 ! second index will be implicitely used to make a scalar product 1679 ! with a potential change, meaning that its complex conjugate must be 1680 ! used. This explains the following signs... 1681 do ipw1=ipw2,npwdiel 1682 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 1683 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 1684 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 1685 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 1686 end do 1687 end do 1688 end if 1689 1690 if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then 1691 ! Accumulate product of band densities by their doccde, for the 1692 ! computation of the effect of change of Fermi level. 1693 do ipw=1,npwdiel 1694 drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1 1695 drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1 1696 end do 1697 ! Also accumulate weighted sum of doccde 1698 sumdocc=sumdocc+wght1 1699 end if 1700 1701 ! End condition of metallic occupancies or extrapolation 1702 end if 1703 1704 ! End loop on spinorial components 1705 end do 1706 ! End loop on iband 1707 end do 1708 1709 call timab(762,2,tsec) 1710 1711 ! -- Compute now off-band-diagonal terms ------------------------------------ 1712 ! -- Compute product of wavefunctions for different bands, inside the block - 1713 1714 call timab(763,1,tsec) 1715 1716 ! if (occopt<3) then 1717 tolocc=1.0d-3 1718 ! else 1719 ! tolocc=1.0d-8 1720 ! end if 1721 1722 if(nbnd_current>1)then 1723 do ibd1=1,nbnd_current-1 1724 ibdshft1=ibd1+iband_shift 1725 do ibd2=ibd1+1,nbnd_current 1726 ibdshft2=ibd2+iband_shift 1727 1728 ! If the occupation numbers are sufficiently different, or 1729 ! if extrapolation is used and the corresponding factor is not zero, 1730 ! then there is a contribution 1731 occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index) 1732 if( abs(occdiff)>tolocc .or. & 1733 & ( extrap==1 .and. & 1734 & ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) & 1735 & ) then 1736 1737 eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index) 1738 1739 ! Loop on spinorial components 1740 do ispinor=1,nspinor 1741 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 1742 1743 ! Store the contribution in wfraug 1744 do i3=1,ndiel3 1745 do i2=1,ndiel2 1746 do i1=1,ndiel1 1747 wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)& 1748 & +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2) 1749 wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)& 1750 & -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2) 1751 end do 1752 end do 1753 end do 1754 1755 ! Performs the Fourier Transform of the product, and store it in wfprod 1756 tim_fourwf=32 1757 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 1758 & 1,kg_diel,kg_diel, mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 1759 & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight) 1760 1761 ! In case of PAW, add compensation charge contribution 1762 if (usepaw==1.and.neglect_pawhat==0) then 1763 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,& 1764 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1765 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 1766 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 1767 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1768 end if 1769 1770 ! Perform now the summation 1771 wght1=zero ; wght2=zero 1772 if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol 1773 if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol 1774 weight=wght1+wght2 1775 1776 if (abs(weight)>tol12) then 1777 do ipw2=1,npwdiel 1778 ! Only fills lower half of the matrix (here, the susceptibility matrix) 1779 ! Note that wfprod of the first index must behave like a density, 1780 ! so that it is used as generated by fourwf, while wfprod of the 1781 ! second index will be implicitely used to make a scalar product 1782 ! with a potential change, meaning that its complex conjugate must be 1783 ! used. This explains the following signs... 1784 do ipw1=ipw2,npwdiel 1785 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 1786 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 1787 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 1788 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 1789 end do 1790 end do 1791 end if 1792 1793 ! End loop on spinorial components 1794 end do 1795 ! End condition of different occupation numbers or extrapolation 1796 end if 1797 ! End internal loop over bands 1798 end do 1799 ! End external loop over bands 1800 end do 1801 ! End condition of having more than one band 1802 end if 1803 1804 ! Loop on secondary block, with fast varying index, in decreasing order. 1805 if(iblk1/=nblk)then 1806 do iblk2=nblk,iblk1+1,-1 1807 iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk 1808 1809 ! Loop over bands to fft and store Fourier transform of wavefunction 1810 iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk 1811 do iband=1,nbnd_in_blk 1812 ! Loop on spinorial components 1813 do ispinor=1,nspinor 1814 iwf=(ispinor-1)*npw_k+(iband+iband_shift2-1)*npw_k*nspinor+icg 1815 1816 ! Obtain Fourier transform in fft box 1817 tim_fourwf=22 1818 cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf) 1819 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 1820 & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,& 1821 & ndiel4,ndiel5,ndiel6,0,tim_fourwf,weight,weight) 1822 wfrspa2(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:) 1823 end do 1824 end do 1825 1826 do ibd1=1,nbnd_current 1827 ibdshft1=ibd1+iband_shift 1828 do ibd2=1,nbnd_in_blk 1829 ibdshft2=ibd2+iband_shift2 1830 1831 ! If the occupation numbers are sufficiently different, or 1832 ! if extrapolation is used and the corresponding factor is not zero, 1833 ! then there is a contribution 1834 occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index) 1835 if( abs(occdiff)>tolocc .or. & 1836 & ( extrap==1 .and. & 1837 & ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) & 1838 & ) then 1839 1840 eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index) 1841 1842 ! Loop on spinorial components 1843 do ispinor=1,nspinor 1844 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 1845 1846 ! Store the contribution in wfraug 1847 do i3=1,ndiel3 1848 do i2=1,ndiel2 1849 do i1=1,ndiel1 1850 wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)& 1851 & +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2) 1852 wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)& 1853 & -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2) 1854 end do 1855 end do 1856 end do 1857 1858 ! Performs the Fourier Transform of the product, and store it in wfprod 1859 tim_fourwf=32 1860 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 1861 & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 1862 & ndiel4,ndiel5,ndiel6,3,tim_fourwf,weight,weight) 1863 1864 ! In case of PAW, add compensation charge contribution 1865 if (usepaw==1.and.neglect_pawhat==0) then 1866 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibdshft2,ispinor,ispinor,1,kg_diel,& 1867 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 1868 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 1869 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 1870 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 1871 end if 1872 1873 ! Perform now the summation 1874 wght1=zero ; wght2=zero 1875 if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol 1876 if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol 1877 weight=wght1+wght2 1878 1879 if (abs(weight)>tol12) then 1880 do ipw2=1,npwdiel 1881 ! Only fills lower half of the matrix (here, the susceptibility matrix) 1882 ! Note that wfprod of the first index must behave like a density, 1883 ! so that it is used as generated by fourwf, while wfprod of the 1884 ! second index will be implicitely used to make a scalar product 1885 ! with a potential change, meaning that its complex conjugate must be 1886 ! used. This explains the following signs... 1887 do ipw1=ipw2,npwdiel 1888 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 1889 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 1890 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 1891 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 1892 end do 1893 end do 1894 end if 1895 1896 ! End loop on spinorial components 1897 end do 1898 ! End condition of different occupation numbers or extrapolation 1899 end if 1900 ! End internal loop over bands 1901 end do 1902 ! End external loop over bands 1903 end do 1904 ! End loop on bloks 1905 end do 1906 1907 ! Finish the loop on blok with iblk2=iblk1+1, so can use the 1908 ! FFTd wavefunctions for the next iblk1. 1909 do iband=1,nbnd_in_blk 1910 wfrspa1(:,:,:,:,1:nspinor,iband)=wfrspa2(:,:,:,:,1:nspinor,iband) 1911 end do 1912 1913 ! End condition of iblk1/=nblk 1914 end if 1915 1916 call timab(763,2,tsec) 1917 1918 ! End loop on iblk1 1919 end do 1920 1921 !DEBUG 1922 !write(std_out,*)' suskmm : exit ' 1923 !do ipw1=1,npwdiel 1924 !write(std_out,*)ipw1,susmat(1,ipw1,1,ipw1,1),susmat(2,ipw1,1,ipw1,1) 1925 !end do 1926 !write(std_out,*)' suskmm : end of susmat ' 1927 !stop 1928 !ENDDEBUG 1929 1930 call destroy_mpi_enreg(mpi_enreg_diel) 1931 ABI_FREE(cwavef) 1932 ABI_FREE(dummy) 1933 ABI_FREE(rhoaug) 1934 ABI_FREE(wfprod) 1935 ABI_FREE(wfraug) 1936 ABI_FREE(wfrspa1) 1937 ABI_FREE(wfrspa2) 1938 1939 call timab(760,2,tsec) 1940 1941 end subroutine suskmm