TABLE OF CONTENTS
ABINIT/suskmm [ 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 !!
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)=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/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
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
PARENTS
suscep_stat
CHILDREN
destroy_mpi_enreg,fourwf,init_distribfft_seq,initmpi_seq,pawsushat timab
SOURCE
105 #if defined HAVE_CONFIG_H 106 #include "config.h" 107 #endif 108 109 #include "abi_common.h" 110 111 112 subroutine suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,& 113 & gbound_diel,gylmg_diel,icg,ikpt,isp,istwf_k,kg_diel,kg_k,& 114 & lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,& 115 & natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,& 116 & npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,paral_kgb,& 117 & pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,& 118 & susmat,typat,ucvol,usepaw,wtk) 119 120 use defs_basis 121 use defs_abitypes 122 use m_errors 123 use m_profiling_abi 124 use m_mpinfo, only : destroy_mpi_enreg 125 126 use m_pawang, only : pawang_type 127 use m_pawtab, only : pawtab_type 128 use m_pawcprj, only : pawcprj_type 129 130 !This section has been created automatically by the script Abilint (TD). 131 !Do not modify the following lines by hand. 132 #undef ABI_FUNC 133 #define ABI_FUNC 'suskmm' 134 use interfaces_18_timing 135 use interfaces_51_manage_mpi 136 use interfaces_53_ffts 137 use interfaces_65_paw 138 !End of the abilint section 139 140 implicit none 141 142 !Arguments ------------------------------------ 143 !This type is defined in defs_mpi 144 !scalars 145 integer,intent(in) :: bdtot_index,extrap,icg,ikpt,isp,istwf_k,lmax_diel,mband,mcg 146 integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat 147 integer,intent(in) :: nkpt,npw_k,npwdiel,nspden,nspden_eff,nspinor 148 integer,intent(in) :: nsppol,ntypat,occopt,paral_kgb,usepaw 149 real(dp),intent(in) :: ucvol 150 real(dp),intent(inout) :: sumdocc 151 type(MPI_type),intent(in) :: mpi_enreg 152 type(pawang_type),intent(in) :: pawang 153 !arrays 154 integer,intent(in) :: atindx(natom),gbound(2*mgfftdiel+8,2) 155 integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2) 156 integer,intent(in) :: kg_diel(3,npwdiel),kg_k(3,npw_k),ngfftdiel(18) 157 integer,intent(in) :: typat(natom) 158 real(dp),intent(in) :: cg(2,mcg),doccde(mband*nkpt*nsppol) 159 real(dp),intent(in) :: eigen(mband*nkpt*nsppol) 160 real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw) 161 real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband) 162 real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt) 163 real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff) 164 real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor) 165 real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden) 166 type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw) 167 type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw) 168 169 !Local variables------------------------------- 170 !scalars 171 integer :: comm_fft,i1,i2,i3,iband,iband_shift,iband_shift2,ibd1,ibd2,ibdshft1,ibdshft2 172 integer :: iblk1,iblk2,ipw,ipw1,ipw2,ispinor,iwf,jsp,mblk 173 integer :: nblk,nbnd_current,nbnd_in_blk,nbnd_in_blk1,ndiel1,ndiel2,ndiel3 174 integer :: testocc,tim_fourwf 175 real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2 176 character(len=500) :: message 177 type(MPI_type) :: mpi_enreg_diel 178 !arrays 179 real(dp) :: tsec(2) 180 real(dp),allocatable :: cwavef(:,:),dummy(:,:),rhoaug(:,:,:),wfprod(:,:) 181 real(dp),allocatable :: wfraug(:,:,:,:),wfrspa1(:,:,:,:,:,:) 182 real(dp),allocatable :: wfrspa2(:,:,:,:,:,:) 183 184 ! ************************************************************************* 185 186 !DEBUG 187 !write(std_out,*)' suskmm : enter ' 188 !if(.true.)stop 189 !ENDDEBUG 190 191 call timab(760,1,tsec) 192 call timab(761,1,tsec) 193 194 !Allocations, initializations 195 ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3) 196 testocc=1 197 ABI_ALLOCATE(rhoaug,(ndiel4,ndiel5,ndiel6)) 198 ABI_ALLOCATE(wfraug,(2,ndiel4,ndiel5,ndiel6)) 199 ABI_ALLOCATE(wfprod,(2,npwdiel)) 200 ABI_ALLOCATE(dummy,(2,1)) 201 202 !The dielectric stuff is performed in sequential mode. 203 !Set mpi_enreg_diel accordingly 204 call initmpi_seq(mpi_enreg_diel) 205 call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all') 206 207 comm_fft=mpi_enreg%comm_fft 208 209 !Prepare the blocking : compute the number of blocks, 210 !the number of bands in each normal block, 211 !and the number in the first one, usually smaller. 212 213 !Consider that if the number of bands is large, there are at most 8 blocks 214 nbnd_in_blk=0 215 if(nband_k>=48)then 216 mblk=8 217 nbnd_in_blk=(nband_k-1)/mblk+1 218 ! If the number of bands is medium, place 6 bands per block 219 else if(nband_k>=12)then 220 nbnd_in_blk=6 221 ! Otherwise, must have at least 2 blocks 222 else if(nband_k>=2)then 223 mblk=2 224 nbnd_in_blk=(nband_k-1)/mblk+1 225 else 226 write(message, '(a,a,a,i2,a,a,a)')& 227 & ' The number of bands must be larger or equal to 2, in suskmm.',ch10,& 228 & ' It is equal to ',nband_k,'.',ch10,& 229 & ' Action : choose another preconditioner.' 230 MSG_ERROR(message) 231 end if 232 233 !Compute the effective number of blocks, and the number of bands in 234 !the first block. 235 nblk=(nband_k-1)/nbnd_in_blk+1 236 nbnd_in_blk1=nband_k-(nblk-1)*nbnd_in_blk 237 238 !DEBUG 239 !write(std_out,*)' suskmm : nband_k,nblk,nbnd_in_blk,nbnd_in_blk1 ' 240 !write(std_out,*)nband_k,nblk,nbnd_in_blk,nbnd_in_blk1 241 !stop 242 !ENDDEBUG 243 244 !wfrspa1 will contain the wavefunctions of the slow sampling (iblk1) 245 ABI_ALLOCATE(wfrspa1,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk)) 246 !wfrspa2 will contain the wavefunctions of the rapid sampling (iblk2) 247 ABI_ALLOCATE(wfrspa2,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk)) 248 249 ABI_ALLOCATE(cwavef,(2,npw_k)) 250 251 call timab(761,2,tsec) 252 253 !First loop over blocks 254 do iblk1=1,nblk 255 256 call timab(762,1,tsec) 257 258 ! Initialisation 259 if(iblk1==1)then 260 261 nbnd_current=nbnd_in_blk1 262 iband_shift=0 263 ! Loop over bands to fft and store Fourier transform of wavefunction 264 do iband=1,nbnd_current 265 ! Loop on spinorial components 266 do ispinor=1,nspinor 267 iwf=(ispinor-1)*npw_k+(iband-1)*npw_k*nspinor+icg 268 ! Obtain Fourier transform in fft box 269 tim_fourwf=21 270 cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf) 271 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 272 & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,& 273 & 0,paral_kgb,tim_fourwf,weight,weight) 274 wfrspa1(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:) 275 end do 276 end do 277 278 else 279 280 ! The Fourier transform of wavefunctions have already been obtained 281 nbnd_current=nbnd_in_blk 282 iband_shift=nbnd_in_blk1+(iblk1-2)*nbnd_in_blk 283 284 end if 285 286 ! Loop over bands of this block, to generate band-diagonal 287 do iband=1,nbnd_current 288 289 ! Loop on spinorial components 290 do ispinor=1,nspinor 291 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 292 293 if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then 294 ! In the case of metallic occupation, or if the extrapolation 295 ! over higher bands is included, must compute the 296 ! Fourier transform of the density of each band, then 297 ! generate the part of the susceptibility matrix due 298 ! varying occupation numbers. 299 weight=-two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol 300 do i3=1,ndiel3 301 do i2=1,ndiel2 302 do i1=1,ndiel1 303 wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,iband)**2& 304 & +wfrspa1(2,i1,i2,i3,ispinor,iband)**2 305 wfraug(2,i1,i2,i3)=zero 306 end do 307 end do 308 ! If extrapolation, accumulate density in real space 309 if(extrap==1.and.usepaw==0)then 310 do i2=1,ndiel2 311 do i1=1,ndiel1 312 rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3) 313 end do 314 end do 315 end if 316 end do 317 318 ! In case of PAW, add compensation charge contribution 319 if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then 320 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,& 321 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 322 & ngfftdiel,npwdiel,nspinor,ntypat,1,& 323 & pawang,pawtab,ph3d_diel,typat,dummy,wfraug,& 324 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 325 rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:) 326 end if 327 328 ! Performs the Fourier Transform of the density of the band, 329 ! and store it in wfprod 330 tim_fourwf=31 331 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 332 & 1,kg_diel,kg_diel,& 333 & mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,ndiel4,ndiel5,ndiel6,3,paral_kgb,tim_fourwf,weight,weight) 334 ! In case of PAW, add compensation charge contribution if not already done 335 if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then 336 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,1,1,1,kg_diel,& 337 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 338 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 339 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 340 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 341 end if 342 343 ! Perform now the summation of terms related to direct change of eigenvalues 344 ! or extrapolation over higher bands 345 wght1=zero ; wght2=zero 346 if(occopt>=3 .and. testocc==1) wght1=doccde(iband+iband_shift+bdtot_index)*wtk(ikpt)/ucvol 347 if(extrap==1) wght2=two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol 348 weight=wght1+wght2 349 350 if (abs(weight)>tol12) then 351 do ipw2=1,npwdiel 352 ! Only fills lower half of the matrix (here, the susceptibility matrix) 353 ! Note that wfprod of the first index must behave like a density, 354 ! so that it is used as generated by fourwf, while wfprod of the 355 ! second index will be implicitely used to make a scalar product 356 ! with a potential change, meaning that its complex conjugate must be 357 ! used. This explains the following signs... 358 do ipw1=ipw2,npwdiel 359 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 360 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 361 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 362 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 363 end do 364 end do 365 end if 366 367 if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then 368 ! Accumulate product of band densities by their doccde, for the 369 ! computation of the effect of change of Fermi level. 370 do ipw=1,npwdiel 371 drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1 372 drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1 373 end do 374 ! Also accumulate weighted sum of doccde 375 sumdocc=sumdocc+wght1 376 end if 377 378 ! End condition of metallic occupancies or extrapolation 379 end if 380 381 ! End loop on spinorial components 382 end do 383 ! End loop on iband 384 end do 385 386 call timab(762,2,tsec) 387 388 ! -- Compute now off-band-diagonal terms ------------------------------------ 389 ! -- Compute product of wavefunctions for different bands, inside the block - 390 391 call timab(763,1,tsec) 392 393 ! if (occopt<3) then 394 tolocc=1.0d-3 395 ! else 396 ! tolocc=1.0d-8 397 ! end if 398 399 if(nbnd_current>1)then 400 do ibd1=1,nbnd_current-1 401 ibdshft1=ibd1+iband_shift 402 do ibd2=ibd1+1,nbnd_current 403 ibdshft2=ibd2+iband_shift 404 405 ! If the occupation numbers are sufficiently different, or 406 ! if extrapolation is used and the corresponding factor is not zero, 407 ! then there is a contribution 408 occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index) 409 if( abs(occdiff)>tolocc .or. & 410 & ( extrap==1 .and. & 411 & ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) & 412 & ) then 413 414 eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index) 415 416 ! Loop on spinorial components 417 do ispinor=1,nspinor 418 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 419 420 ! Store the contribution in wfraug 421 do i3=1,ndiel3 422 do i2=1,ndiel2 423 do i1=1,ndiel1 424 wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)& 425 & +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2) 426 wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)& 427 & -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2) 428 end do 429 end do 430 end do 431 432 ! Performs the Fourier Transform of the product, and store it in wfprod 433 tim_fourwf=32 434 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 435 & 1,kg_diel,kg_diel, mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 436 & ndiel4,ndiel5,ndiel6,3,paral_kgb,tim_fourwf,weight,weight) 437 438 ! In case of PAW, add compensation charge contribution 439 if (usepaw==1.and.neglect_pawhat==0) then 440 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,& 441 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 442 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 443 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 444 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 445 end if 446 447 ! Perform now the summation 448 wght1=zero ; wght2=zero 449 if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol 450 if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol 451 weight=wght1+wght2 452 453 if (abs(weight)>tol12) then 454 do ipw2=1,npwdiel 455 ! Only fills lower half of the matrix (here, the susceptibility matrix) 456 ! Note that wfprod of the first index must behave like a density, 457 ! so that it is used as generated by fourwf, while wfprod of the 458 ! second index will be implicitely used to make a scalar product 459 ! with a potential change, meaning that its complex conjugate must be 460 ! used. This explains the following signs... 461 do ipw1=ipw2,npwdiel 462 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 463 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 464 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 465 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 466 end do 467 end do 468 end if 469 470 ! End loop on spinorial components 471 end do 472 ! End condition of different occupation numbers or extrapolation 473 end if 474 ! End internal loop over bands 475 end do 476 ! End external loop over bands 477 end do 478 ! End condition of having more than one band 479 end if 480 481 ! Loop on secondary block, with fast varying index, in decreasing order. 482 if(iblk1/=nblk)then 483 do iblk2=nblk,iblk1+1,-1 484 iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk 485 486 ! Loop over bands to fft and store Fourier transform of wavefunction 487 iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk 488 do iband=1,nbnd_in_blk 489 ! Loop on spinorial components 490 do ispinor=1,nspinor 491 iwf=(ispinor-1)*npw_k+(iband+iband_shift2-1)*npw_k*nspinor+icg 492 493 ! Obtain Fourier transform in fft box 494 tim_fourwf=22 495 cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf) 496 call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,& 497 & istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,& 498 & ndiel4,ndiel5,ndiel6,0,paral_kgb,tim_fourwf,weight,weight) 499 wfrspa2(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:) 500 end do 501 end do 502 503 do ibd1=1,nbnd_current 504 ibdshft1=ibd1+iband_shift 505 do ibd2=1,nbnd_in_blk 506 ibdshft2=ibd2+iband_shift2 507 508 ! If the occupation numbers are sufficiently different, or 509 ! if extrapolation is used and the corresponding factor is not zero, 510 ! then there is a contribution 511 occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index) 512 if( abs(occdiff)>tolocc .or. & 513 & ( extrap==1 .and. & 514 & ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) & 515 & ) then 516 517 eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index) 518 519 ! Loop on spinorial components 520 do ispinor=1,nspinor 521 jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp 522 523 ! Store the contribution in wfraug 524 do i3=1,ndiel3 525 do i2=1,ndiel2 526 do i1=1,ndiel1 527 wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)& 528 & +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2) 529 wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)& 530 & -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2) 531 end do 532 end do 533 end do 534 535 ! Performs the Fourier Transform of the product, and store it in wfprod 536 tim_fourwf=32 537 call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,& 538 & 1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,& 539 & ndiel4,ndiel5,ndiel6,3,paral_kgb,tim_fourwf,weight,weight) 540 541 ! In case of PAW, add compensation charge contribution 542 if (usepaw==1.and.neglect_pawhat==0) then 543 call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibdshft2,ispinor,ispinor,1,kg_diel,& 544 & lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,& 545 & ngfftdiel,npwdiel,nspinor,ntypat,0,& 546 & pawang,pawtab,ph3d_diel,typat,wfprod,dummy,& 547 & mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom) 548 end if 549 550 ! Perform now the summation 551 wght1=zero ; wght2=zero 552 if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol 553 if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol 554 weight=wght1+wght2 555 556 if (abs(weight)>tol12) then 557 do ipw2=1,npwdiel 558 ! Only fills lower half of the matrix (here, the susceptibility matrix) 559 ! Note that wfprod of the first index must behave like a density, 560 ! so that it is used as generated by fourwf, while wfprod of the 561 ! second index will be implicitely used to make a scalar product 562 ! with a potential change, meaning that its complex conjugate must be 563 ! used. This explains the following signs... 564 do ipw1=ipw2,npwdiel 565 susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+& 566 & weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2)) 567 susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+& 568 & weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2)) 569 end do 570 end do 571 end if 572 573 ! End loop on spinorial components 574 end do 575 ! End condition of different occupation numbers or extrapolation 576 end if 577 ! End internal loop over bands 578 end do 579 ! End external loop over bands 580 end do 581 ! End loop on bloks 582 end do 583 584 ! Finish the loop on blok with iblk2=iblk1+1, so can use the 585 ! FFTd wavefunctions for the next iblk1. 586 do iband=1,nbnd_in_blk 587 wfrspa1(:,:,:,:,1:nspinor,iband)=wfrspa2(:,:,:,:,1:nspinor,iband) 588 end do 589 590 ! End condition of iblk1/=nblk 591 end if 592 593 call timab(763,2,tsec) 594 595 ! End loop on iblk1 596 end do 597 598 !DEBUG 599 !write(std_out,*)' suskmm : exit ' 600 !do ipw1=1,npwdiel 601 !write(std_out,*)ipw1,susmat(1,ipw1,1,ipw1,1),susmat(2,ipw1,1,ipw1,1) 602 !end do 603 !write(std_out,*)' suskmm : end of susmat ' 604 !stop 605 !ENDDEBUG 606 607 call destroy_mpi_enreg(mpi_enreg_diel) 608 ABI_DEALLOCATE(cwavef) 609 ABI_DEALLOCATE(dummy) 610 ABI_DEALLOCATE(rhoaug) 611 ABI_DEALLOCATE(wfprod) 612 ABI_DEALLOCATE(wfraug) 613 ABI_DEALLOCATE(wfrspa1) 614 ABI_DEALLOCATE(wfrspa2) 615 616 call timab(760,2,tsec) 617 618 end subroutine suskmm