TABLE OF CONTENTS


ABINIT/m_suscep_stat [ Modules ]

[ Top ] [ Modules ]

NAME

 m_suscep_stat

FUNCTION

 Compute the susceptibility matrix

COPYRIGHT

  Copyright (C) 2008-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 MODULE m_suscep_stat
28 
29  use defs_basis
30  use defs_abitypes
31  use m_xmpi
32  use m_errors
33  use m_abicore
34 
35  use m_time,    only : timab
36  use m_pawang,  only : pawang_type
37  use m_pawtab,  only : pawtab_type
38  use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_mpi_allgather, pawcprj_free
39  use m_mpinfo,  only : destroy_mpi_enreg, initmpi_seq, proc_distrb_cycle
40  use m_kg,      only : ph1d3d
41  use m_gsphere, only : symg
42  use m_fftcore, only : sphereboundary
43  use m_fft,     only : fftpac, fourwf
44  use m_spacepar,     only : symrhg
45  use m_paw_finegrid, only : pawgylmg
46  use m_paw_nhat,     only : pawsushat
47 
48  implicit none
49 
50  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)

PARENTS

      vtorho

CHILDREN

      destroy_mpi_enreg,fftpac,init_distribfft_seq,initmpi_seq,pawcprj_alloc
      pawcprj_free,pawcprj_get,pawcprj_mpi_allgather,pawgylmg,ph1d3d
      sphereboundary,susk,suskmm,symg,symrhg,timab,xmpi_sum,zhpev

SOURCE

160 subroutine suscep_stat(atindx,atindx1,cg,cprj,dielar,dimcprj,doccde,&
161 &  eigen,gbound_diel,gprimd,irrzondiel,istwfk,kg,&
162 &  kg_diel,lmax_diel,&
163 &  mband,mcg,mcprj,mgfftdiel,mkmem,mpi_enreg,mpw,natom,nband,&
164 &  neglect_pawhat,nfftdiel,ngfftdiel,nkpt,npwarr,&
165 &  npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occ,occopt,&
166 &  pawang,pawtab,phnonsdiel,ph1ddiel,rprimd,&
167 &  susmat,symafm,symrel,tnons,typat,ucvol,unpaw,usecprj,usepaw,usetimerev,&
168 &  wtk,ylmdiel)
169 
170 
171 !This section has been created automatically by the script Abilint (TD).
172 !Do not modify the following lines by hand.
173 #undef ABI_FUNC
174 #define ABI_FUNC 'suscep_stat'
175 !End of the abilint section
176 
177  implicit none
178 
179 !Arguments ------------------------------------
180 !scalars
181  integer,intent(in) :: lmax_diel,mband,mcg,mcprj,mgfftdiel,mkmem,mpw,natom,neglect_pawhat
182  integer,intent(in) :: nfftdiel,nkpt,npwdiel,nspden,nspinor,nsppol,nsym,ntypat,occopt
183  integer,intent(in) :: unpaw,usecprj,usepaw,usetimerev
184  real(dp),intent(in) :: ucvol
185  type(MPI_type),intent(in) :: mpi_enreg
186  type(pawang_type),intent(in) :: pawang
187 !arrays
188  integer,intent(in) :: atindx(natom),atindx1(natom),dimcprj(natom*usepaw)
189  integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2)
190 !no_abirules
191 !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise
192  integer,intent(in) :: irrzondiel(nfftdiel**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4)),&
193  & istwfk(nkpt)
194  integer,intent(in) :: kg(3,mpw*mkmem),kg_diel(3,npwdiel),&
195  & nband(nkpt*nsppol),ngfftdiel(18)
196  integer,intent(in) :: npwarr(nkpt),symafm(nsym),symrel(3,3,nsym),typat(ntypat)
197  real(dp),intent(in) :: cg(2,mcg),dielar(7)
198  real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol)
199  real(dp),intent(in) :: gprimd(3,3),occ(mband*nkpt*nsppol)
200 !nfftdiel**(1-1/nsym) is 1 if nsym==1, and nfftdiel otherwise
201  real(dp),intent(in) :: phnonsdiel(2,nfftdiel**(1-1/nsym),(nspden/nsppol)-3*(nspden/4)),&
202  &                                 tnons(3,nsym),wtk(nkpt)
203  real(dp),intent(in) :: ph1ddiel(2,(3*(2*mgfftdiel+1)*natom)*usepaw),rprimd(3,3)
204  real(dp),intent(in) :: ylmdiel(npwdiel,lmax_diel**2)
205  real(dp),intent(out) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
206  type(pawcprj_type) :: cprj(natom,mcprj*usecprj)
207  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
208 
209 !Local variables-------------------------------
210 !scalars
211  integer :: bdtot_index,diag,extrap,i1,i2,i3,iband,ibg,icg,ier,ierr
212  integer :: ifft,ii,ikg,ikpt,indx,iorder_cprj,ipw1,ipw2,isp,isp1,isp2
213  integer :: ispinor,istwf_k,isym,j1,j2,j3,jj,jsp,k1,k2,k3
214  integer :: my_nspinor,nband_k,nband_loc,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6
215  integer :: nkpg_diel,npw_k,npwsp,nspden_eff,nspden_tmp,nsym1,nsym2
216  integer :: paral_kgb_diel,spaceComm,t1,t2,testocc
217  real(dp) :: ai,ai2,ar,ar2,diegap,dielam,emax,invnsym
218  real(dp) :: invnsym1,invnsym2,phi1,phi12,phi2,phr1,phr12
219  real(dp) :: phr2,sumdocc,weight
220  logical :: antiferro
221  character(len=500) :: message
222  type(MPI_type) :: mpi_enreg_diel
223 
224 !arrays
225  integer,allocatable :: gbound(:,:),kg_k(:,:),sym_g(:,:)
226  integer,allocatable :: tmrev_g(:)
227  real(dp) :: kpt_diel(3,1),tsec(2)
228  real(dp),allocatable :: drhode(:,:,:),drhode_wk(:,:,:)
229  real(dp),allocatable :: eig_diel(:),gylmg_diel(:,:,:),kpg_dum(:,:)
230  real(dp),allocatable :: occ_deavg(:),ph3d_diel(:,:,:),phdiel(:,:,:)
231  real(dp),allocatable :: phkxred_diel(:,:),rhoextrap(:,:,:,:),rhoextrg(:,:)
232  real(dp),allocatable :: rhoextrr(:,:),sush(:),sussum(:),susvec(:,:,:)
233  real(dp),allocatable :: suswk(:,:,:,:),zhpev1(:,:),zhpev2(:)
234  type(pawcprj_type),allocatable :: cprj_k(:,:),cprj_loc(:,:)
235 
236 ! *************************************************************************
237 
238  call timab(740,1,tsec)
239  call timab(741,1,tsec)
240 
241  ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
242 
243 
244 !----- Initialisations -----------------------------------------------------------
245 !---------------------------------------------------------------------------------
246 
247  if (usecprj==0.and.usepaw==1) then
248    write (message,'(3a)')&
249 &   ' cprj datastructure must be allocated !',ch10,&
250 &   ' Action: change pawusecp input keyword.'
251    MSG_ERROR(message)
252  end if
253 
254  if (mpi_enreg%paral_spinor==1) then
255    message = ' not yet allowed for parallelization over spinors !'
256    MSG_ERROR(message)
257  end if
258 
259 !Init mpicomm
260  if(mpi_enreg%paral_kgb==1) then
261    spaceComm=mpi_enreg%comm_kpt
262  else
263    spaceComm=mpi_enreg%comm_cell
264  end if
265 
266 !The dielectric stuff is performed in sequential mode.
267 !Set mpi_enreg_diel accordingly
268  call initmpi_seq(mpi_enreg_diel)
269  call init_distribfft_seq(MPI_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all')
270 
271 !testocc to be taken away
272  testocc=1
273 
274  my_nspinor=max(1,nspinor/mpi_enreg%nproc_spinor)
275 
276  iorder_cprj=0 ! order for the cprj reading...
277 
278 !Initialize some scalar quantities
279  antiferro=(nsppol==1.and.nspden==2)
280  nspden_eff=min(max(nsppol,nspden),2) ! Size for the computed part of susmat
281  bdtot_index=0 ; icg=0 ; ibg=0
282  ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
283 
284 !ndiel4,ndiel5,ndiel6 are FFT dimensions, modified to avoid cache trashing
285  ndiel4=ngfftdiel(4) ; ndiel5=ngfftdiel(5) ; ndiel6=ngfftdiel(6)
286  diegap=dielar(5) ; dielam=dielar(6)
287  extrap=0
288 
289 !If dielam is too small, there is no extrapolation.
290  if(dielam>1.0d-6)extrap=1
291 
292 !Some stuff for symmetries
293  nsym1=sum(symafm,mask=symafm==1)
294  nsym2=nsym-nsym1
295  invnsym =one/dble(nsym)
296  invnsym1=one/dble(nsym1)
297  invnsym2=one
298 !FIXME: make sure this is consistent with following code
299 !div by 0 for several v5 tests
300  if (nsym2 > 0) invnsym2=one/dble(nsym2)
301 
302 !Allocations
303  ABI_ALLOCATE(occ_deavg,(mband))
304  if(occopt>=3) then
305    ABI_ALLOCATE(drhode,(2,npwdiel,nspden_eff))
306  else
307    ABI_ALLOCATE(drhode,(0,0,0))
308  end if
309  if(extrap==1) then
310    ABI_ALLOCATE(rhoextrap,(ndiel4,ndiel5,ndiel6,nspinor))
311  else
312    ABI_ALLOCATE(rhoextrap,(0,0,0,0))
313  end if
314 
315 !zero the susceptibility matrix and other needed quantities
316  susmat(:,:,:,:,:)=zero
317  if(occopt>=3)then
318    drhode(:,:,:)=zero
319    sumdocc=zero
320  end if
321 
322 !PAW additional initializations
323  if (usepaw==1) then
324    ABI_ALLOCATE(gylmg_diel,(npwdiel,lmax_diel**2,ntypat))
325    ABI_ALLOCATE(ph3d_diel,(2,npwdiel,natom))
326    if (neglect_pawhat==0) then
327      ABI_ALLOCATE(phkxred_diel,(2,natom))
328      ABI_ALLOCATE(kpg_dum,(0,0))
329      kpt_diel(1:3,1)=zero;phkxred_diel(1,:)=one;phkxred_diel(2,:)=zero;nkpg_diel=0
330 !    write(std_out,*) ' lmax_diel ', lmax_diel
331      call pawgylmg(gprimd,gylmg_diel,kg_diel,kpg_dum,kpt_diel,lmax_diel,nkpg_diel,npwdiel,ntypat,pawtab,ylmdiel)
332      call ph1d3d(1,natom,kg_diel,natom,natom,npwdiel,ndiel1,ndiel2,ndiel3,phkxred_diel,ph1ddiel,ph3d_diel)
333      ABI_DEALLOCATE(phkxred_diel)
334      ABI_DEALLOCATE(kpg_dum)
335    else
336      gylmg_diel=zero;ph3d_diel=one
337    end if
338  else
339    ABI_ALLOCATE(gylmg_diel,(0,0,0))
340    ABI_ALLOCATE(ph3d_diel,(0,0,0))
341  end if
342 
343  call timab(741,2,tsec)
344 
345 
346 
347 !--BIG loop over spins ------------------------------------------------------------
348 !---------------------------------------------------------------------------------
349 
350  do isp=1,nsppol
351    ikg=0
352 
353    if(extrap==1)rhoextrap(:,:,:,:)=zero
354 
355 !  --BIG loop over k-points --------------------------------------------------------
356 !  ---------------------------------------------------------------------------------
357 
358    do ikpt=1,nkpt
359 
360      nband_k=nband(ikpt+(isp-1)*nkpt)
361      istwf_k=istwfk(ikpt)
362      npw_k=npwarr(ikpt)
363 
364      if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isp,mpi_enreg%me_kpt)) then
365        bdtot_index=bdtot_index+nband_k
366        cycle
367      end if
368 
369      call timab(742,1,tsec)
370 
371 
372      ABI_ALLOCATE(gbound,(2*mgfftdiel+8,2))
373      ABI_ALLOCATE(kg_k,(3,npw_k))
374 
375      if (usepaw==1) then
376        ABI_DATATYPE_ALLOCATE(cprj_k,(natom,my_nspinor*nband_k))
377        if (neglect_pawhat==0) then
378          call pawcprj_alloc(cprj_k,0,dimcprj)
379          if (mpi_enreg%nproc_band==1) then
380            call pawcprj_get(atindx1,cprj_k,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,&
381 &           mband,mkmem,natom,nband_k,nband_k,my_nspinor,nsppol,unpaw,&
382 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
383          else
384            nband_loc=nband_k/mpi_enreg%nproc_band
385            ABI_DATATYPE_ALLOCATE(cprj_loc,(natom,my_nspinor*nband_loc))
386            call pawcprj_alloc(cprj_loc,0,dimcprj)
387            call pawcprj_get(atindx1,cprj_loc,cprj,natom,1,ibg,ikpt,iorder_cprj,isp,&
388 &           mband/mpi_enreg%nproc_band,mkmem,natom,nband_loc,nband_loc,my_nspinor,nsppol,unpaw,&
389 &           mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb)
390            call pawcprj_mpi_allgather(cprj_loc,cprj_k,natom,my_nspinor*nband_loc,dimcprj,0,&
391 &           mpi_enreg%nproc_band,mpi_enreg%comm_band,ierr,rank_ordered=.true.)
392            call pawcprj_free(cprj_loc)
393            ABI_DATATYPE_DEALLOCATE(cprj_loc)
394          end if
395        else
396          !call pawcprj_nullify(cprj_k)
397        end if
398      else
399        ABI_DATATYPE_ALLOCATE(cprj_k,(0,0))
400      end if
401 
402      kg_k(:,1:npw_k)=kg(:,1+ikg:npw_k+ikg)
403      call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k)
404 
405      if(extrap==1)then
406 !      Compute inverse of average dielectric gap for each band
407 !      and multiply by occupation factor
408        emax=maxval(eigen(1+bdtot_index:nband_k+bdtot_index))
409        do iband=1,nband_k
410          occ_deavg(iband)= occ(iband+bdtot_index)*dielam &
411 &         / ( emax-eigen(iband+bdtot_index)  + diegap )
412        end do
413      else
414        occ_deavg(:)=zero
415      end if
416 
417      call timab(742,2,tsec)
418      call timab(743,1,tsec)
419 
420 !    Compute the contribution of each k-point to susmat, rhoextrap, drhode and sumdocc.
421      if(mpi_enreg%paral_kgb==1)then !Only this version is in parallel
422 !      Use either the simpler implementation
423 !      Should provide a test !!!
424        call susk(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,&
425 &       gbound_diel,gylmg_diel,icg,ikpt,&
426 &       isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
427 &       natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
428 &       npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,&
429 &       pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
430 &       susmat,typat,ucvol,usepaw,wtk)
431      else
432 !      Or the more sophisticated one, needed to save memory.
433        call suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,&
434 &       gbound_diel,gylmg_diel,icg,ikpt,&
435 &       isp,istwf_k,kg_diel,kg_k,lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
436 &       natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
437 &       npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,paral_kgb_diel,&
438 &       pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
439 &       susmat,typat,ucvol,usepaw,wtk)
440      end if
441 
442      call timab(743,2,tsec)
443 
444      ABI_DEALLOCATE(gbound)
445      ABI_DEALLOCATE(kg_k)
446 
447      bdtot_index=bdtot_index+nband_k
448 
449      if (mkmem/=0) then
450        ibg=ibg+my_nspinor*nband_k
451        icg=icg+my_nspinor*npw_k*nband_k
452        ikg=ikg+npw_k
453      end if
454      if (usepaw==1) then
455        if (neglect_pawhat==0) then
456          call pawcprj_free(cprj_k)
457        end if
458      end if
459      ABI_DATATYPE_DEALLOCATE(cprj_k)
460 
461 !    End loop on ikpt:  --------------------------------------------------------
462    end do
463 
464 !  Here include the contribution from the extrapolation to susmat,
465 !  diagonal part
466    if(extrap==1)then
467 
468      call timab(744,1,tsec)
469 
470 !    Transfer extrapolating density on augmented fft grid to
471 !    normal fft grid in real space.
472 !    Warning1 : if collinear magnetism, must treat only one spin at a time
473 !    Warning2 : if non-collinear magnetism, treat both spins
474 !    Warning3 : this is subtle for antiferro magnetism
475      nspden_tmp=1;if (antiferro) nspden_tmp=2
476      ABI_ALLOCATE(rhoextrr,(nfftdiel,nspden_tmp))
477      ABI_ALLOCATE(rhoextrg,(2,nfftdiel))
478      if (nspden==1.and.nspinor==2) rhoextrap(:,:,:,1)=rhoextrap(:,:,:,1)+rhoextrap(:,:,:,2)
479 
480      do ispinor=1,min(nspinor,nspden)
481        jsp=isp+ispinor-1
482 
483        call fftpac(1,mpi_enreg_diel,1,ndiel1,ndiel2,ndiel3,ndiel4,ndiel5,ndiel6,&
484 &       ngfftdiel,rhoextrr(:,1),rhoextrap(:,:,:,ispinor),1)
485 
486 !      Generate the density in reciprocal space, and symmetrize it
487 !      (note symrhg also make the reverse FFT, to get symmetrized density;
488 !      this is useless here, and should be made an option)
489        call symrhg(1,gprimd,irrzondiel,mpi_enreg_diel,nfftdiel,nfftdiel,ngfftdiel,&
490 &       nspden_tmp,1,nsym,paral_kgb_diel,phnonsdiel,rhoextrg,rhoextrr,rprimd,symafm,symrel)
491 
492        do ipw2=1,npwdiel
493          j1=kg_diel(1,ipw2) ; j2=kg_diel(2,ipw2) ; j3=kg_diel(3,ipw2)
494 !        static:    Only fills lower half of the matrix (here, the susceptibility matrix)
495 !        dynamical: fill all, will not affect susopt==2 for which extrap==0
496          do ipw1=1,npwdiel
497            i1=kg_diel(1,ipw1) ; i2=kg_diel(2,ipw1) ; i3=kg_diel(3,ipw1)
498 !          NOTE that there is a FFT folding (superposition) bias here
499 !          Should use kgindex, in the same spirit as in prcref
500            k1=i1-j1; k1=modulo(k1,ndiel1)
501            k2=i2-j2; k2=modulo(k2,ndiel2)
502            k3=i3-j3; k3=modulo(k3,ndiel3)
503            ifft=k1+1+ndiel1*(k2+ndiel2*k3)
504            susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+rhoextrg(1,ifft)
505            susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+rhoextrg(2,ifft)
506          end do
507        end do
508 
509      end do
510      ABI_DEALLOCATE(rhoextrg)
511      ABI_DEALLOCATE(rhoextrr)
512 
513      call timab(744,2,tsec)
514 
515    end if
516 
517 !  End loop over spins ---------------------------------------------------------
518  end do
519 
520  ABI_DEALLOCATE(occ_deavg)
521  ABI_DEALLOCATE(rhoextrap)
522  ABI_DEALLOCATE(gylmg_diel)
523  ABI_DEALLOCATE(ph3d_diel)
524  !end if
525 
526  call destroy_mpi_enreg(mpi_enreg_diel)
527 
528 !-- Stuff for parallelism --------------------------------------------------------
529 !---------------------------------------------------------------------------------
530 
531  if(xmpi_paral==1)then
532    call timab(746,1,tsec)
533    ABI_ALLOCATE(sussum,(2*npwdiel*nspden*npwdiel*nspden))
534 !  Recreate full susmat on all proc.
535 !  This should be coded more efficiently,
536 !  since half of the matrix is still empty, and
537 !  it is spin-diagonal.
538    sussum(:)=reshape(susmat(:,:,:,:,:),(/2*npwdiel*nspden*npwdiel*nspden/))
539    call xmpi_sum(sussum,spaceComm,ierr)
540    susmat(:,:,:,:,:)=reshape(sussum(:),(/2,npwdiel,nspden,npwdiel,nspden/))
541    ABI_DEALLOCATE(sussum)
542 !  Recreate full drhode on all proc.
543    if(occopt>=3 .and. testocc==1)then
544      call xmpi_sum(drhode,spaceComm,ierr)
545 !    Should use only one mpi-allreduce call instead of the three
546      call xmpi_sum(sumdocc,spaceComm,ierr)
547    end if
548    call timab(746,2,tsec)
549  end if
550 
551 !-- Apply spatial hermitian/symmetries on spin-diagonal susceptibility matrix ----
552 !---------------------------------------------------------------------------------
553 
554  call timab(747,1,tsec)
555 
556 !If antiferro magnetism, has to divide (spin-diagonal) susceptibility by 2 (due to dble occupations)
557  if (antiferro) then
558    do ipw2=1,npwdiel
559      do ipw1=ipw2,npwdiel
560        susmat(:,ipw1,1,ipw2,1)=half*susmat(:,ipw1,1,ipw2,1)
561      end do
562    end do
563  end if
564 
565 !Generate upper half of the spin-diagonal matrix (still the susceptibility matrix)
566  do isp=1,nspden_eff
567    do ipw2=2,npwdiel
568      do ipw1=1,ipw2-1
569        susmat(1,ipw1,isp,ipw2,isp)= susmat(1,ipw2,isp,ipw1,isp)
570        susmat(2,ipw1,isp,ipw2,isp)=-susmat(2,ipw2,isp,ipw1,isp)
571      end do
572    end do
573  end do
574 
575 !Compute symmetric of G-vectors and eventual phases
576 !(either time-reversal or spatial symmetries)
577  ABI_ALLOCATE(tmrev_g,(npwdiel))
578  ABI_ALLOCATE(sym_g,(npwdiel,nsym))
579  ABI_ALLOCATE(phdiel,(2,npwdiel,nsym))
580  call symg(kg_diel,npwdiel,nsym,phdiel,sym_g,symrel,tmrev_g,tnons)
581 
582 !Impose spatial symmetries to the spin-diagonal susceptibility matrix
583  ABI_ALLOCATE(suswk,(2,npwdiel,npwdiel,nspden_eff))
584  do isp=1,nspden_eff
585    suswk(:,:,:,isp)=susmat(:,:,isp,:,isp) ! Temporary storage
586  end do
587 
588  do isp=1,nspden_eff
589    jsp=min(3-isp,nsppol)
590    do ipw2=1,npwdiel
591      do ipw1=1,npwdiel
592        ar=suswk(1,ipw1,ipw2,isp)
593        ai=suswk(2,ipw1,ipw2,isp)
594        ar2=zero;ai2=zero
595        if(nsym>1)then
596          do isym=2,nsym
597            t1=sym_g(ipw1,isym) ; t2=sym_g(ipw2,isym)
598 !          Not all symmetries are non-symmorphic. Should save time here ...
599            phr1=phdiel(1,ipw1,isym) ; phi1=phdiel(2,ipw1,isym)
600            phr2=phdiel(1,ipw2,isym) ; phi2=phdiel(2,ipw2,isym)
601            phr12= phr1*phr2+phi1*phi2 ; phi12=phi1*phr2-phr1*phi2
602            if (symafm(isym)==1) then
603              ar=ar+suswk(1,t1,t2,isp)*phr12-suswk(2,t1,t2,isp)*phi12
604              ai=ai+suswk(2,t1,t2,isp)*phr12+suswk(1,t1,t2,isp)*phi12
605            else
606              ar2=ar2+suswk(1,t1,t2,jsp)*phr12-suswk(2,t1,t2,jsp)*phi12
607              ai2=ai2+suswk(2,t1,t2,jsp)*phr12+suswk(1,t1,t2,jsp)*phi12
608            end if
609          end do
610        end if
611        if (antiferro) then
612          susmat(1,ipw1,1,ipw2,1)=ar*invnsym1
613          susmat(2,ipw1,1,ipw2,1)=ai*invnsym1
614          susmat(1,ipw1,2,ipw2,2)=ar2*invnsym2
615          susmat(2,ipw1,2,ipw2,2)=ai2*invnsym2
616        else
617          susmat(1,ipw1,isp,ipw2,isp)=(ar+ar2)*invnsym
618          susmat(2,ipw1,isp,ipw2,isp)=(ai+ai2)*invnsym
619        end if
620      end do
621    end do
622  end do
623  ABI_DEALLOCATE(suswk)
624 
625 
626 !--  Add contribibution to susceptibility due to change of Fermi level -----------
627 !---------------------------------------------------------------------------------
628 
629  if (occopt>=3.and.testocc==1) then
630 
631 !  Impose spatial symmetries to drhode
632    ABI_ALLOCATE(drhode_wk,(2,npwdiel,nspden_eff))
633    do isp=1,nspden_eff
634      jsp=min(3-isp,nsppol)
635      do ipw1=1,npwdiel
636        ar=drhode(1,ipw1,isp)
637        ai=drhode(2,ipw1,isp)
638        ar2=zero;ai2=zero
639        if (nsym>1) then
640          do isym=2,nsym
641            t1=sym_g(ipw1,isym)
642 !          Not all symmetries are non-symmorphic. Should save time here ...
643            phr1=phdiel(1,ipw1,isym);phi1=phdiel(2,ipw1,isym)
644            if (symafm(isym)==1) then
645              ar=ar+drhode(1,t1,isp)*phr1-drhode(2,t1,isp)*phi1
646              ai=ai+drhode(2,t1,isp)*phr1+drhode(1,t1,isp)*phi1
647            else
648              ar2=ar2+drhode(1,t1,jsp)*phr1-drhode(2,t1,jsp)*phi1
649              ai2=ai2+drhode(2,t1,jsp)*phr1+drhode(1,t1,jsp)*phi1
650            end if
651          end do
652        end if
653        if (antiferro) then  ! 1/2 factor due to (dble) occupations
654          drhode_wk(1,ipw1,1)=half*ar*invnsym1
655          drhode_wk(2,ipw1,1)=half*ai*invnsym1
656          drhode_wk(1,ipw1,2)=half*ar2*invnsym2
657          drhode_wk(2,ipw1,2)=half*ai2*invnsym2
658        else
659          drhode_wk(1,ipw1,isp)=(ar+ar2)*invnsym
660          drhode_wk(2,ipw1,isp)=(ai+ai2)*invnsym
661        end if
662      end do
663    end do
664 
665 !  Add contribution to non-diagonal susceptibility
666 !  Presently fills complete susceptibility matrix, not only lower half
667    weight=one/sumdocc
668    do isp2=1,nspden_eff
669      do ipw2=1,npwdiel
670        do isp1=1,nspden_eff
671          do ipw1=1,npwdiel
672            susmat(1,ipw1,isp1,ipw2,isp2)=susmat(1,ipw1,isp1,ipw2,isp2)- &
673 &           weight*( drhode_wk(1,ipw1,isp1)*drhode_wk(1,ipw2,isp2)  &
674 &           +drhode_wk(2,ipw1,isp1)*drhode_wk(2,ipw2,isp2) )
675            susmat(2,ipw1,isp1,ipw2,isp2)=susmat(2,ipw1,isp1,ipw2,isp2)- &
676 &           weight*( drhode_wk(2,ipw1,isp1)*drhode_wk(1,ipw2,isp2)  &
677 &           -drhode_wk(1,ipw1,isp1)*drhode_wk(2,ipw2,isp2) )
678          end do
679        end do
680      end do
681    end do
682    ABI_DEALLOCATE(drhode_wk)
683 
684  end if
685  !if (occopt>=3)  then
686  ABI_DEALLOCATE(drhode)
687  !end if
688 
689 
690 !--- Impose the time-reversal symmetry to the susceptibility matrix --------------
691 !---------------------------------------------------------------------------------
692 
693  if (usetimerev==1) then
694    ABI_ALLOCATE(suswk,(2,npwdiel,npwdiel,1))
695 
696 !  Impose the time-reversal symmetry to the spin-diagonal susceptibility matrix
697    do isp=1,nspden_eff
698      suswk(:,:,:,1)=susmat(:,:,isp,:,isp) ! Temporary storage
699      do ipw2=1,npwdiel
700        t2=tmrev_g(ipw2)
701        do ipw1=1,npwdiel
702          t1=tmrev_g(ipw1)
703          susmat(1,ipw1,isp,ipw2,isp)=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1))
704          susmat(2,ipw1,isp,ipw2,isp)=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1))
705        end do
706      end do
707    end do
708 
709 !  Impose the time-reversal symmetry to the off-diagonal susceptibility matrix
710    if (nspden_eff/=1.and.occopt>=3.and.testocc==1) then
711      suswk(:,:,:,1)=susmat(:,:,1,:,2) ! Temporary storage
712      do ipw2=1,npwdiel
713        t2=tmrev_g(ipw2)
714        do ipw1=1,npwdiel
715          t1=tmrev_g(ipw1)
716          ar=half*(suswk(1,ipw1,ipw2,1)+suswk(1,t1,t2,1))
717          ai=half*(suswk(2,ipw1,ipw2,1)-suswk(2,t1,t2,1))
718          susmat(1,ipw1,1,ipw2,2)= ar
719          susmat(2,ipw1,1,ipw2,2)= ai
720          susmat(1,ipw1,2,ipw2,1)= ar
721          susmat(2,ipw1,2,ipw2,1)=-ai
722        end do
723      end do
724    end if
725    ABI_DEALLOCATE(suswk)
726  end if
727 
728  ABI_DEALLOCATE(phdiel)
729  ABI_DEALLOCATE(sym_g)
730  ABI_DEALLOCATE(tmrev_g)
731 
732 
733 !-- The full susceptibility matrix is computed -----------------------------------
734 !-- Now, eventually diagonalize it and stop --------------------------------------
735 !---------------------------------------------------------------------------------
736 
737 !Must turn on this flag to make the diagonalisation
738  diag=0
739  if(diag==1)then
740 
741    npwsp=npwdiel*nspden_eff
742    ABI_ALLOCATE(sush,(npwsp*(npwsp+1)))
743    ABI_ALLOCATE(susvec,(2,npwsp,npwsp))
744    ABI_ALLOCATE(eig_diel,(npwsp))
745    ABI_ALLOCATE(zhpev1,(2,2*npwsp-1))
746    ABI_ALLOCATE(zhpev2,(3*npwsp-2))
747    ier=0
748 
749 !  Store the susceptibility matrix in proper mode before calling zhpev
750    indx=1
751    do ii=1,npwdiel
752      do jj=1,ii
753        sush(indx  )=susmat(1,jj,1,ii,1)
754        sush(indx+1)=susmat(2,jj,1,ii,1)
755        indx=indx+2
756      end do
757    end do
758 
759 !  If spin-polarized, need to store other parts of the matrix
760    if(nspden_eff/=1)then
761      do ii=1,npwdiel
762 !      Here, spin-flip contribution
763        do jj=1,npwdiel
764          sush(indx  )=susmat(1,jj,1,ii,2)
765          sush(indx+1)=susmat(2,jj,1,ii,2)
766          indx=indx+2
767        end do
768 !      Here spin down-spin down upper matrix
769        do jj=1,ii
770          sush(indx  )=susmat(1,jj,2,ii,2)
771          sush(indx+1)=susmat(2,jj,2,ii,2)
772          indx=indx+2
773        end do
774      end do
775    end if
776 
777    call ZHPEV ('V','U',npwsp,sush,eig_diel,susvec,npwsp,zhpev1,&
778 &   zhpev2,ier)
779 
780    write(std_out,*)' suscep_stat : print eigenvalues of the susceptibility matrix'
781    do ii=1,npwsp
782      write(std_out,'(i5,es16.6)' )ii,eig_diel(ii)
783    end do
784 
785    ABI_DEALLOCATE(sush)
786    ABI_DEALLOCATE(susvec)
787    ABI_DEALLOCATE(eig_diel)
788    ABI_DEALLOCATE(zhpev1)
789    ABI_DEALLOCATE(zhpev2)
790    MSG_ERROR("Stopping here!")
791  end if
792 
793  call timab(747,2,tsec)
794  call timab(740,2,tsec)
795 
796 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

PARENTS

      suscep_stat

CHILDREN

      destroy_mpi_enreg,fourwf,init_distribfft_seq,initmpi_seq,pawsushat
      sphereboundary,timab,xmpi_allgather,xmpi_allgatherv,xmpi_alltoallv
      xmpi_sum

SOURCE

 906 subroutine susk(atindx,bdtot_index,cg_mpi,cprj_k,doccde,drhode,eigen,extrap,gbound,&
 907 &  gbound_diel,gylmg_diel,icg_mpi,ikpt,isp,istwf_k,kg_diel,kg_k_mpi,&
 908 &  lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
 909 &  natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
 910 &  npwdiel,npw_k_mpi,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,&
 911 &  pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
 912 &  susmat,typat,ucvol,usepaw,wtk)
 913 
 914 
 915 !This section has been created automatically by the script Abilint (TD).
 916 !Do not modify the following lines by hand.
 917 #undef ABI_FUNC
 918 #define ABI_FUNC 'susk'
 919 !End of the abilint section
 920 
 921  implicit none
 922 
 923 !Arguments ------------------------------------
 924 !This type is defined in defs_mpi
 925 !scalars
 926  integer,intent(in) :: bdtot_index,extrap,ikpt,isp,istwf_k,lmax_diel,mband,mcg
 927  integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat
 928  integer,intent(in) :: nkpt,npwdiel,nspden,nspden_eff,nspinor,nsppol
 929  integer,intent(in) :: ntypat,occopt,usepaw
 930  integer,intent(in),target :: icg_mpi,npw_k_mpi
 931  real(dp),intent(in) :: ucvol
 932  real(dp),intent(inout) :: sumdocc
 933  type(MPI_type),intent(in) :: mpi_enreg
 934  type(pawang_type),intent(in) :: pawang
 935 !arrays
 936  integer,intent(in) :: atindx(natom),gbound_diel(2*mgfftdiel+8,2)
 937  integer,intent(in) :: kg_diel(3,npwdiel),ngfftdiel(18),typat(natom)
 938  integer,intent(in),target :: kg_k_mpi(3,npw_k_mpi)
 939  integer,intent(inout) :: gbound(2*mgfftdiel+8,2)
 940  integer,pointer :: kg_k(:,:)
 941  real(dp),intent(in) :: doccde(mband*nkpt*nsppol),eigen(mband*nkpt*nsppol)
 942  real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw)
 943  real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband)
 944  real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt)
 945  real(dp),intent(in),target :: cg_mpi(2,mcg)
 946  real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff)
 947  real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)
 948  real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
 949  type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw)
 950  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
 951 
 952 !Local variables-------------------------------
 953 ! real(dp), allocatable :: cg_disk(:,:)
 954 !Local variables for MPI
 955 !scalars
 956  integer :: blocksize,i1,i2,i3,iband,iband_loc,ibd1,ibd2,ibdblock,ier
 957  integer :: iproc,iproc_fft,ipw,ipw1,ipw2,isp1,isp2,ispinor,iwf,jsp,me_bandfft
 958  integer :: nbdblock,ndatarecv,ndiel1,ndiel2,ndiel3
 959  integer :: sizemax_per_proc,spaceComm,testocc,tim_fourwf
 960  integer,pointer :: icg,npw_k
 961  integer,target :: icg_loc=0,npw_k_loc,npw_tot
 962  real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2
 963  type(MPI_type) :: mpi_enreg_diel
 964 !arrays
 965  integer,allocatable :: band_loc(:),kg_k_gather(:,:),npw_per_proc(:),rdispls(:)
 966  integer,allocatable :: rdispls_all(:),rdisplsloc(:),recvcounts(:)
 967  integer,allocatable :: recvcountsloc(:),sdispls(:),sdisplsloc(:),sendcounts(:)
 968  integer,allocatable :: sendcountsloc(:),susmat_mpi(:,:,:)
 969  integer,allocatable,target :: kg_k_gather_all(:,:)
 970  real(dp) :: tsec(2)
 971  real(dp),allocatable :: cwavef(:,:),cwavef_alltoall(:,:)
 972  real(dp),allocatable :: cwavef_alltoall_gather(:,:),dummy(:,:),rhoaug(:,:,:)
 973  real(dp),allocatable :: wfprod(:,:),wfraug(:,:,:,:),wfrspa(:,:,:,:,:,:)
 974  real(dp),allocatable,target :: cg_local(:,:)
 975  real(dp),pointer :: cg(:,:)
 976  logical,allocatable :: treat_band(:)
 977 
 978 ! *************************************************************************
 979 
 980 !DEBUG
 981 !write(std_out,*)' susk : enter '
 982 !if(.true.)stop
 983 !ENDDEBUG
 984 
 985  call timab(750,1,tsec)
 986  call timab(751,1,tsec)
 987 
 988  ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
 989 
 990 !The dielectric stuff is performed in sequential mode.
 991 !Set mpi_enreg_diel accordingly
 992  call initmpi_seq(mpi_enreg_diel)
 993  call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all')
 994  me_bandfft=xmpi_comm_rank(mpi_enreg%comm_bandfft)
 995 
 996  testocc=1
 997 !DEBUG
 998 !write(std_out,*)' susk : set testocc to 0 '
 999 !testocc=0
1000 !write(std_out,*)' susk : set extrap to 0 '
1001 !extrap=0
1002 !ENDDEBUG
1003 
1004 !Allocations, initializations
1005  ABI_ALLOCATE(rhoaug,(ndiel4,ndiel5,ndiel6))
1006  ABI_ALLOCATE(wfraug,(2,ndiel4,ndiel5,ndiel6))
1007  ABI_ALLOCATE(wfprod,(2,npwdiel))
1008  ABI_ALLOCATE(wfrspa,(2,ndiel4,ndiel5,ndiel6,nspinor,mband))
1009  ABI_ALLOCATE(dummy,(2,1))
1010  wfrspa(:,:,:,:,:,:)=zero
1011  ABI_ALLOCATE(treat_band,(nband_k))
1012  treat_band(:)=.true.
1013  isp1=isp;isp2=isp
1014  if (nspden_eff==2.and.nspinor==2) isp2=isp+1
1015 
1016 !BAND-FFT parallelism
1017  if (mpi_enreg%paral_kgb==1) then
1018    treat_band(:)=.false.
1019 !  We gather the wavefunctions treated by this proc in cg_local
1020    spaceComm=mpi_enreg%comm_band
1021    blocksize=mpi_enreg%nproc_band
1022    nbdblock=nband_k/blocksize
1023    ABI_ALLOCATE(sdispls,(blocksize))
1024    ABI_ALLOCATE(sdisplsloc,(blocksize))
1025    ABI_ALLOCATE(sendcounts,(blocksize))
1026    ABI_ALLOCATE(sendcountsloc,(blocksize))
1027    ABI_ALLOCATE(rdispls,(blocksize))
1028    ABI_ALLOCATE(rdisplsloc,(blocksize))
1029    ABI_ALLOCATE(recvcounts,(blocksize))
1030    ABI_ALLOCATE(recvcountsloc,(blocksize))
1031 !  First gather the kg_k in kg_k_gather_all
1032    npw_k_loc=npw_k_mpi
1033    call xmpi_allgather(npw_k_loc,recvcounts,spaceComm,ier)
1034    rdispls(1)=0
1035    do iproc=2,blocksize
1036      rdispls(iproc)=rdispls(iproc-1)+recvcounts(iproc-1)
1037    end do
1038    ndatarecv=rdispls(blocksize)+recvcounts(blocksize)
1039    ABI_ALLOCATE(kg_k_gather,(3,ndatarecv))
1040    recvcountsloc(:)=recvcounts(:)*3
1041    rdisplsloc(:)=rdispls(:)*3
1042    call xmpi_allgatherv(kg_k_mpi,3*npw_k_loc,kg_k_gather,recvcountsloc(:),rdisplsloc,spaceComm,ier)
1043    ABI_ALLOCATE(npw_per_proc,(mpi_enreg%nproc_fft))
1044    ABI_ALLOCATE(rdispls_all,(mpi_enreg%nproc_fft))
1045    spaceComm=mpi_enreg%comm_fft
1046    call xmpi_allgather(ndatarecv,npw_per_proc,spaceComm,ier)
1047    rdispls_all(1)=0
1048    do iproc=2,mpi_enreg%nproc_fft
1049      rdispls_all(iproc)=rdispls_all(iproc-1)+npw_per_proc(iproc-1)
1050    end do
1051    npw_tot=rdispls_all(mpi_enreg%nproc_fft)+npw_per_proc(mpi_enreg%nproc_fft)
1052    ABI_ALLOCATE(kg_k_gather_all,(3,npw_tot))
1053    call xmpi_allgatherv(kg_k_gather,3*ndatarecv,kg_k_gather_all,3*npw_per_proc(:),3*rdispls_all,spaceComm,ier)
1054 !  At this point kg_k_gather_all contains all the kg
1055    if(allocated(cwavef))  then
1056      ABI_DEALLOCATE(cwavef)
1057    end if
1058    ABI_ALLOCATE(cwavef,(2,npw_k_loc*nspinor*blocksize))
1059    sizemax_per_proc=nband_k/(mpi_enreg%nproc_band*mpi_enreg%nproc_fft)+1
1060    ABI_ALLOCATE(band_loc,(sizemax_per_proc))
1061    ABI_ALLOCATE(cg_local,(2,sizemax_per_proc*npw_tot*nspinor))
1062    iband_loc=0
1063    do ibdblock=1,nbdblock
1064      cwavef(:,1:npw_k_loc*nspinor*blocksize)=&
1065 &     cg_mpi(:,1+(ibdblock-1)*npw_k_loc*nspinor*blocksize+icg_mpi:ibdblock*npw_k_loc*nspinor*blocksize+icg_mpi)
1066      sendcounts(:)=npw_k_loc
1067      do iproc=1,blocksize
1068        sdispls(iproc)=(iproc-1)*npw_k_loc
1069      end do
1070      ABI_ALLOCATE(cwavef_alltoall,(2,ndatarecv*nspinor))
1071      recvcountsloc(:)=recvcounts(:)*2*nspinor
1072      rdisplsloc(:)=rdispls(:)*2*nspinor
1073      sendcountsloc(:)=sendcounts(:)*2*nspinor
1074      sdisplsloc(:)=sdispls(:)*2*nspinor
1075      call timab(547,1,tsec)
1076      spaceComm=mpi_enreg%comm_band
1077      call xmpi_alltoallv(cwavef,sendcountsloc,sdisplsloc,cwavef_alltoall,recvcountsloc,rdisplsloc,spaceComm,ier)
1078      call timab(547,2,tsec)
1079      ABI_ALLOCATE(cwavef_alltoall_gather,(2,npw_tot*nspinor))
1080      blocksize=mpi_enreg%nproc_band
1081      spaceComm=mpi_enreg%comm_fft
1082      call xmpi_allgatherv(cwavef_alltoall,2*nspinor*ndatarecv,cwavef_alltoall_gather,&
1083 &     2*nspinor*npw_per_proc,2*nspinor*rdispls_all,spaceComm,ier)
1084      iproc_fft=modulo(ibdblock-1,mpi_enreg%nproc_fft)
1085      if(mpi_enreg%me_fft==iproc_fft) then !All nproc_band procs of index me_fft will treat these bands
1086        iband_loc=iband_loc+1
1087        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
1088        treat_band(iband)=.true.
1089        band_loc(iband_loc)=iband
1090        cg_local(:,1+(iband_loc-1)*npw_tot*nspinor:iband_loc*npw_tot*nspinor)=cwavef_alltoall_gather(:,1:npw_tot*nspinor)
1091      end if
1092      ABI_DEALLOCATE(cwavef_alltoall_gather)
1093      ABI_DEALLOCATE(cwavef_alltoall)
1094    end do
1095 !  On exit:
1096 !  npw_tot will be npw
1097 !  kg_k_gather_all will be kg_k
1098 !  cg_local will be cg
1099 !  icg will be zero
1100    npw_k=>npw_tot
1101    kg_k=>kg_k_gather_all(:,:)
1102    cg=>cg_local(:,:)
1103    icg=>icg_loc
1104    call sphereboundary(gbound,istwf_k,kg_k,mgfftdiel,npw_k)
1105    ABI_DEALLOCATE(npw_per_proc)
1106    ABI_DEALLOCATE(rdispls_all)
1107    ABI_DEALLOCATE(sendcounts)
1108    ABI_DEALLOCATE(recvcounts)
1109    ABI_DEALLOCATE(sdispls)
1110    ABI_DEALLOCATE(rdispls)
1111    ABI_DEALLOCATE(sendcountsloc)
1112    ABI_DEALLOCATE(sdisplsloc)
1113    ABI_DEALLOCATE(recvcountsloc)
1114    ABI_DEALLOCATE(rdisplsloc)
1115    ABI_DEALLOCATE(kg_k_gather)
1116    ABI_DEALLOCATE(cwavef)
1117 !  Because they will be summed over all procs, and arrive on input, rescale drhode and rhoextrap
1118    if(occopt>=3)drhode(:,:,isp1:isp2)=drhode(:,:,isp1:isp2)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1119    if(extrap==1)rhoextrap(:,:,:,:)=rhoextrap(:,:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1120    do i1=isp1,isp2
1121      susmat(:,:,i1,:,i1)=susmat(:,:,i1,:,i1)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1122    end do
1123 
1124 !  No BAND-FFT parallelism
1125  else ! use argument variables
1126    cg=>cg_mpi
1127    kg_k=>kg_k_mpi
1128    npw_k=>npw_k_mpi
1129    icg=>icg_mpi
1130  end if
1131  iband_loc=0
1132 
1133  call timab(751,2,tsec)
1134  call timab(752,1,tsec)
1135 
1136 !Loop over bands to fft and store Fourier transform of wavefunction
1137  ABI_ALLOCATE(cwavef,(2,npw_k))
1138  do iband=1,nband_k
1139    if(.not. treat_band(iband))  cycle ! I am not treating this band (only for the parallel case)
1140    iband_loc=iband_loc+1
1141 
1142 !  Loop on spinorial components
1143    do ispinor=1,nspinor
1144      iwf=(ispinor-1)*npw_k+(iband_loc-1)*npw_k*nspinor+icg
1145      jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1146 
1147 !    Obtain Fourier transform in fft box
1148      tim_fourwf=8
1149      cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf)
1150      call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
1151 &     istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,&
1152 &     0,mpi_enreg_diel%paral_kgb,tim_fourwf,weight,weight)
1153 
1154      wfrspa(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:)
1155 
1156      if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then
1157 !      In the case of metallic occupation, or if the extrapolation
1158 !      over higher bands is included, must compute the
1159 !      Fourier transform of the density of each band, then
1160 !      generate the part of the susceptibility matrix due
1161 !      varying occupation numbers.
1162 
1163        weight=-two*occ_deavg(iband)*wtk(ikpt)/ucvol
1164        do i3=1,ndiel3
1165          do i2=1,ndiel2
1166            do i1=1,ndiel1
1167              wfraug(1,i1,i2,i3)=wfraug(1,i1,i2,i3)**2+wfraug(2,i1,i2,i3)**2
1168              wfraug(2,i1,i2,i3)=zero
1169            end do
1170          end do
1171 !        If extrapolation, accumulate density in real space
1172          if(extrap==1.and.usepaw==0)then
1173            do i2=1,ndiel2
1174              do i1=1,ndiel1
1175                rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3)
1176              end do
1177            end do
1178          end if
1179        end do
1180 
1181 !      In case of PAW, add compensation charge contribution
1182        if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then
1183          call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,&
1184 &         lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1185 &         ngfftdiel,npwdiel,nspinor,ntypat,1,&
1186 &         pawang,pawtab,ph3d_diel,typat,dummy,wfraug,&
1187 &         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1188          rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:)
1189        end if
1190 
1191 !      Performs the Fourier Transform of the density of the band,
1192 !      and store it in wfprod
1193        tim_fourwf=9
1194        call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1195 &       1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1196 &       ndiel4,ndiel5,ndiel6,3,mpi_enreg_diel%paral_kgb,tim_fourwf,weight,weight)
1197 !      In case of PAW, add compensation charge contribution if not already done
1198        if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then
1199          call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,&
1200 &         lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1201 &         ngfftdiel,npwdiel,nspinor,ntypat,0,&
1202 &         pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1203 &         mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1204        end if
1205 
1206 !      Perform now the summation of terms related to direct change of eigenvalues
1207 !      or extrapolation over higher bands
1208        wght1=zero ; wght2=zero
1209        if(occopt>=3 .and. testocc==1) wght1=doccde(iband+bdtot_index)*wtk(ikpt)/ucvol
1210        if(extrap==1) wght2=two*occ_deavg(iband)*wtk(ikpt)/ucvol
1211        weight=wght1+wght2
1212 
1213        if (abs(weight)>tol12) then
1214          do ipw2=1,npwdiel
1215 !          Only fills lower half of the matrix (here, the susceptibility matrix)
1216 !          Note that wfprod of the first index must behave like a density,
1217 !          so that it is used as generated by fourwf, while wfprod of the
1218 !          second index will be implicitely used to make a scalar product
1219 !          with a potential change, meaning that its complex conjugate must be
1220 !          used. This explains the following signs...
1221            do ipw1=ipw2,npwdiel
1222              susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1223 &             weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1224              susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1225 &             weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1226            end do
1227          end do
1228        end if
1229 
1230        if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then
1231 !        Accumulate product of band densities by their doccde, for the
1232 !        computation of the effect of change of Fermi level.
1233          do ipw=1,npwdiel
1234            drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1
1235            drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1
1236          end do
1237 !        Also accumulate weighted sum of doccde
1238          sumdocc=sumdocc+wght1
1239        end if
1240 
1241 !      End condition of metallic occupancies or extrapolation
1242      end if
1243 
1244 !    End loop on spinorial components
1245    end do
1246 !  End loop on iband
1247  end do
1248 
1249  call timab(752,2,tsec)
1250  call timab(753,1,tsec)
1251 
1252  ABI_DEALLOCATE(cwavef)
1253 
1254 !Stuff for parallelism (bands-FFT)
1255  if(mpi_enreg%paral_kgb==1) then
1256    call xmpi_sum(wfrspa,mpi_enreg%comm_bandfft,ier)
1257    if(occopt>=3) then
1258      call xmpi_sum(drhode(:,:,isp1:isp2),mpi_enreg%comm_bandfft,ier)
1259    end if
1260    if(extrap==1) then
1261      call xmpi_sum(rhoextrap,mpi_enreg%comm_bandfft,ier)
1262    end if
1263    if(occopt>=3) then
1264      call xmpi_sum(sumdocc,mpi_enreg%comm_bandfft,ier)
1265    end if
1266    ABI_ALLOCATE(susmat_mpi,(2,npwdiel,npwdiel))
1267    do i1=isp1,isp2
1268      susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1)
1269      call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier)
1270      susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:)/real(mpi_enreg%nproc_fft*mpi_enreg%nproc_band,dp)
1271    end do
1272    ABI_DEALLOCATE(susmat_mpi)
1273  end if
1274  call timab(753,2,tsec)
1275 
1276 !-- Wavefunctions have been generated in real space ------------------------
1277 !-- Now, compute product of wavefunctions for different bands --------------
1278  call timab(754,1,tsec)
1279 !if (occopt<3) then
1280  tolocc=1.0d-3
1281 !else
1282 !tolocc=1.0d-8
1283 !end if
1284  iproc=-1
1285 
1286  if(nband_k>1)then
1287    do ibd1=1,nband_k-1
1288      do ibd2=ibd1+1,nband_k
1289        iproc=iproc+1
1290        if(modulo(iproc,mpi_enreg%nproc_fft*mpi_enreg%nproc_band) /= me_bandfft) cycle
1291 !      If the occupation numbers are sufficiently different, or
1292 !      if extrapolation is used and the corresponding factor is not zero,
1293 !      then there is a contribution
1294        occdiff=occ(ibd1+bdtot_index)-occ(ibd2+bdtot_index)
1295        if( abs(occdiff)>tolocc      .or. &
1296 &       ( extrap==1 .and.            &
1297 &       ( abs(occ_deavg(ibd1)) + abs(occ_deavg(ibd2)) ) >tolocc ) &
1298 &       ) then
1299 
1300          eigdiff=eigen(ibd1+bdtot_index)-eigen(ibd2+bdtot_index)
1301 !        DEBUG
1302 !        write(std_out,*)' susk : contribution from bands',ibd1,ibd2
1303 !        write(std_out,*)'   occ diff =',occdiff
1304 !        write(std_out,*)'   eig diff =',eigdiff
1305 !        ENDDEBUG
1306 
1307 !        Loop on spinorial components
1308          do ispinor=1,nspinor
1309            jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1310 
1311 !          Store the contribution in wfraug
1312            do i3=1,ndiel3
1313              do i2=1,ndiel2
1314                do i1=1,ndiel1
1315                  wfraug(1,i1,i2,i3)=wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)&
1316 &                 +wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2)
1317                  wfraug(2,i1,i2,i3)=wfrspa(2,i1,i2,i3,ispinor,ibd1)*wfrspa(1,i1,i2,i3,ispinor,ibd2)&
1318 &                 -wfrspa(1,i1,i2,i3,ispinor,ibd1)*wfrspa(2,i1,i2,i3,ispinor,ibd2)
1319                end do
1320              end do
1321            end do
1322 
1323 !          Performs the Fourier Transform of the product, and store it in wfprod
1324            tim_fourwf=19
1325            call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1326 &           1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1327 &           ndiel4,ndiel5,ndiel6,3,mpi_enreg_diel%paral_kgb,tim_fourwf,weight,weight)
1328 
1329 !          In case of PAW, add compensation charge contribution
1330            if (usepaw==1.and.neglect_pawhat==0) then
1331              call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,&
1332 &             lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1333 &             ngfftdiel,npwdiel,nspinor,ntypat,0,&
1334 &             pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1335 &             mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1336            end if
1337 
1338 !          Perform now the summation
1339            wght1=zero ; wght2=zero
1340            if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol
1341            if(extrap==1) wght2=(occ_deavg(ibd1)+occ_deavg(ibd2)) * two*wtk(ikpt)/ucvol
1342            weight=wght1+wght2
1343 
1344 !          DEBUG
1345 !          write(std_out,*)' weight =',weight
1346 !          norm=zero
1347 !          do ipw=1,npwdiel
1348 !          norm=norm+wfprod(1,ipw)**2+wfprod(2,ipw)**2
1349 !          end do
1350 !          write(std_out,*)' norm in reciprocal space  =',norm
1351 !          ENDDEBUG
1352 
1353            if (abs(weight)>tol12) then
1354              do ipw2=1,npwdiel
1355 !              Only fills lower half of the matrix (here, the susceptibility matrix)
1356 !              Note that wfprod of the first index must behave like a density,
1357 !              so that it is used as generated by fourwf, while wfprod of the
1358 !              second index will be implicitely used to make a scalar product
1359 !              with a potential change, meaning that its complex conjugate must be
1360 !              used. This explains the following signs...
1361                do ipw1=ipw2,npwdiel
1362                  susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1363 &                 weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1364                  susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1365 &                 weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1366                end do
1367              end do
1368            end if
1369 
1370 !          End loop on spinorial components
1371          end do
1372 !        End condition of different occupation numbers or extrapolation
1373        end if
1374 !      End internal loop over bands
1375      end do
1376 !    End external loop over bands
1377    end do
1378 !  End condition of having more than one band
1379  end if
1380 
1381  call timab(754,2,tsec)
1382  call timab(755,1,tsec)
1383 
1384  if(mpi_enreg%paral_kgb==1) then
1385    ABI_ALLOCATE(susmat_mpi,(2,npwdiel,npwdiel))
1386    do i1=isp1,isp2
1387      susmat_mpi(:,:,:)=susmat(:,:,i1,:,i1)
1388      call xmpi_sum(susmat_mpi,mpi_enreg%comm_bandfft,ier)
1389      susmat(:,:,i1,:,i1)=susmat_mpi(:,:,:)
1390    end do
1391    ABI_DEALLOCATE(susmat_mpi)
1392    ABI_DEALLOCATE(band_loc)
1393    ABI_DEALLOCATE(treat_band)
1394    ABI_DEALLOCATE(cg_local)
1395    ABI_DEALLOCATE(kg_k_gather_all)
1396  end if
1397 
1398  call destroy_mpi_enreg(mpi_enreg_diel)
1399  ABI_DEALLOCATE(dummy)
1400  ABI_DEALLOCATE(rhoaug)
1401  ABI_DEALLOCATE(wfprod)
1402  ABI_DEALLOCATE(wfraug)
1403  ABI_DEALLOCATE(wfrspa)
1404 
1405  call timab(755,2,tsec)
1406  call timab(750,2,tsec)
1407 
1408 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

PARENTS

      suscep_stat

CHILDREN

      destroy_mpi_enreg,fourwf,init_distribfft_seq,initmpi_seq,pawsushat
      timab

SOURCE

1507 subroutine suskmm(atindx,bdtot_index,cg,cprj_k,doccde,drhode,eigen,extrap,gbound,&
1508 &  gbound_diel,gylmg_diel,icg,ikpt,isp,istwf_k,kg_diel,kg_k,&
1509 &  lmax_diel,mband,mcg,mgfftdiel,mpi_enreg,&
1510 &  natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat,ngfftdiel,nkpt,&
1511 &  npwdiel,npw_k,nspden,nspden_eff,nspinor,nsppol,ntypat,occ,occopt,occ_deavg,paral_kgb,&
1512 &  pawang,pawtab,ph3d_diel,rhoextrap,sumdocc,&
1513 &  susmat,typat,ucvol,usepaw,wtk)
1514 
1515 
1516 !This section has been created automatically by the script Abilint (TD).
1517 !Do not modify the following lines by hand.
1518 #undef ABI_FUNC
1519 #define ABI_FUNC 'suskmm'
1520 !End of the abilint section
1521 
1522  implicit none
1523 
1524 !Arguments ------------------------------------
1525 !scalars
1526  integer,intent(in) :: bdtot_index,extrap,icg,ikpt,isp,istwf_k,lmax_diel,mband,mcg
1527  integer,intent(in) :: mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,neglect_pawhat
1528  integer,intent(in) :: nkpt,npw_k,npwdiel,nspden,nspden_eff,nspinor
1529  integer,intent(in) :: nsppol,ntypat,occopt,paral_kgb,usepaw
1530  real(dp),intent(in) :: ucvol
1531  real(dp),intent(inout) :: sumdocc
1532  type(MPI_type),intent(in) :: mpi_enreg
1533  type(pawang_type),intent(in) :: pawang
1534 !arrays
1535  integer,intent(in) :: atindx(natom),gbound(2*mgfftdiel+8,2)
1536  integer,intent(in) :: gbound_diel(2*mgfftdiel+8,2)
1537  integer,intent(in) :: kg_diel(3,npwdiel),kg_k(3,npw_k),ngfftdiel(18)
1538  integer,intent(in) :: typat(natom)
1539  real(dp),intent(in) :: cg(2,mcg),doccde(mband*nkpt*nsppol)
1540  real(dp),intent(in) :: eigen(mband*nkpt*nsppol)
1541  real(dp),intent(in) :: gylmg_diel(npwdiel,lmax_diel**2,ntypat*usepaw)
1542  real(dp),intent(in) :: occ(mband*nkpt*nsppol),occ_deavg(mband)
1543  real(dp),intent(in) :: ph3d_diel(2,npwdiel,natom*usepaw),wtk(nkpt)
1544  real(dp),intent(inout) :: drhode(2,npwdiel,nspden_eff)
1545  real(dp),intent(inout) :: rhoextrap(ndiel4,ndiel5,ndiel6,nspinor)
1546  real(dp),intent(inout) :: susmat(2,npwdiel,nspden,npwdiel,nspden)
1547  type(pawcprj_type) :: cprj_k(natom,nspinor*nband_k*usepaw)
1548  type(pawtab_type),intent(in) :: pawtab(ntypat*usepaw)
1549 
1550 !Local variables-------------------------------
1551 !scalars
1552  integer :: comm_fft,i1,i2,i3,iband,iband_shift,iband_shift2,ibd1,ibd2,ibdshft1,ibdshft2
1553  integer :: iblk1,iblk2,ipw,ipw1,ipw2,ispinor,iwf,jsp,mblk
1554  integer :: nblk,nbnd_current,nbnd_in_blk,nbnd_in_blk1,ndiel1,ndiel2,ndiel3
1555  integer :: testocc,tim_fourwf
1556  real(dp) :: eigdiff,occdiff,tolocc,weight,wght1,wght2
1557  character(len=500) :: message
1558  type(MPI_type) :: mpi_enreg_diel
1559 !arrays
1560  real(dp) :: tsec(2)
1561  real(dp),allocatable :: cwavef(:,:),dummy(:,:),rhoaug(:,:,:),wfprod(:,:)
1562  real(dp),allocatable :: wfraug(:,:,:,:),wfrspa1(:,:,:,:,:,:)
1563  real(dp),allocatable :: wfrspa2(:,:,:,:,:,:)
1564 
1565 ! *************************************************************************
1566 
1567 !DEBUG
1568 !write(std_out,*)' suskmm : enter '
1569 !if(.true.)stop
1570 !ENDDEBUG
1571 
1572  call timab(760,1,tsec)
1573  call timab(761,1,tsec)
1574 
1575 !Allocations, initializations
1576  ndiel1=ngfftdiel(1) ; ndiel2=ngfftdiel(2) ; ndiel3=ngfftdiel(3)
1577  testocc=1
1578  ABI_ALLOCATE(rhoaug,(ndiel4,ndiel5,ndiel6))
1579  ABI_ALLOCATE(wfraug,(2,ndiel4,ndiel5,ndiel6))
1580  ABI_ALLOCATE(wfprod,(2,npwdiel))
1581  ABI_ALLOCATE(dummy,(2,1))
1582 
1583 !The dielectric stuff is performed in sequential mode.
1584 !Set mpi_enreg_diel accordingly
1585  call initmpi_seq(mpi_enreg_diel)
1586  call init_distribfft_seq(mpi_enreg_diel%distribfft,'c',ngfftdiel(2),ngfftdiel(3),'all')
1587 
1588  comm_fft=mpi_enreg%comm_fft
1589 
1590 !Prepare the blocking : compute the number of blocks,
1591 !the number of bands in each normal block,
1592 !and the number in the first one, usually smaller.
1593 
1594 !Consider that if the number of bands is large, there are at most 8 blocks
1595  nbnd_in_blk=0
1596  if(nband_k>=48)then
1597    mblk=8
1598    nbnd_in_blk=(nband_k-1)/mblk+1
1599 !  If the number of bands is medium, place 6 bands per block
1600  else if(nband_k>=12)then
1601    nbnd_in_blk=6
1602 !  Otherwise, must have at least 2 blocks
1603  else if(nband_k>=2)then
1604    mblk=2
1605    nbnd_in_blk=(nband_k-1)/mblk+1
1606  else
1607    write(message, '(a,a,a,i2,a,a,a)')&
1608 &   '  The number of bands must be larger or equal to 2, in suskmm.',ch10,&
1609 &   '  It is equal to ',nband_k,'.',ch10,&
1610 &   '  Action : choose another preconditioner.'
1611    MSG_ERROR(message)
1612  end if
1613 
1614 !Compute the effective number of blocks, and the number of bands in
1615 !the first block.
1616  nblk=(nband_k-1)/nbnd_in_blk+1
1617  nbnd_in_blk1=nband_k-(nblk-1)*nbnd_in_blk
1618 
1619 !DEBUG
1620 !write(std_out,*)' suskmm : nband_k,nblk,nbnd_in_blk,nbnd_in_blk1 '
1621 !write(std_out,*)nband_k,nblk,nbnd_in_blk,nbnd_in_blk1
1622 !stop
1623 !ENDDEBUG
1624 
1625 !wfrspa1 will contain the wavefunctions of the slow sampling (iblk1)
1626  ABI_ALLOCATE(wfrspa1,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk))
1627 !wfrspa2 will contain the wavefunctions of the rapid sampling (iblk2)
1628  ABI_ALLOCATE(wfrspa2,(2,ndiel4,ndiel5,ndiel6,nspinor,nbnd_in_blk))
1629 
1630  ABI_ALLOCATE(cwavef,(2,npw_k))
1631 
1632  call timab(761,2,tsec)
1633 
1634 !First loop over blocks
1635  do iblk1=1,nblk
1636 
1637    call timab(762,1,tsec)
1638 
1639 !  Initialisation
1640    if(iblk1==1)then
1641 
1642      nbnd_current=nbnd_in_blk1
1643      iband_shift=0
1644 !    Loop over bands to fft and store Fourier transform of wavefunction
1645      do iband=1,nbnd_current
1646 !      Loop on spinorial components
1647        do ispinor=1,nspinor
1648          iwf=(ispinor-1)*npw_k+(iband-1)*npw_k*nspinor+icg
1649 !        Obtain Fourier transform in fft box
1650          tim_fourwf=21
1651          cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf)
1652          call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
1653 &         istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,ndiel4,ndiel5,ndiel6,&
1654 &         0,paral_kgb,tim_fourwf,weight,weight)
1655          wfrspa1(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:)
1656        end do
1657      end do
1658 
1659    else
1660 
1661 !    The Fourier transform of wavefunctions have already been obtained
1662      nbnd_current=nbnd_in_blk
1663      iband_shift=nbnd_in_blk1+(iblk1-2)*nbnd_in_blk
1664 
1665    end if
1666 
1667 !  Loop over bands of this block, to generate band-diagonal
1668    do iband=1,nbnd_current
1669 
1670 !    Loop on spinorial components
1671      do ispinor=1,nspinor
1672        jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1673 
1674        if( (occopt>=3 .and. testocc==1) .or. extrap==1 )then
1675 !        In the case of metallic occupation, or if the extrapolation
1676 !        over higher bands is included, must compute the
1677 !        Fourier transform of the density of each band, then
1678 !        generate the part of the susceptibility matrix due
1679 !        varying occupation numbers.
1680          weight=-two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol
1681          do i3=1,ndiel3
1682            do i2=1,ndiel2
1683              do i1=1,ndiel1
1684                wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,iband)**2&
1685 &               +wfrspa1(2,i1,i2,i3,ispinor,iband)**2
1686                wfraug(2,i1,i2,i3)=zero
1687              end do
1688            end do
1689 !          If extrapolation, accumulate density in real space
1690            if(extrap==1.and.usepaw==0)then
1691              do i2=1,ndiel2
1692                do i1=1,ndiel1
1693                  rhoextrap(i1,i2,i3,ispinor)=rhoextrap(i1,i2,i3,ispinor)+weight*wfraug(1,i1,i2,i3)
1694                end do
1695              end do
1696            end if
1697          end do
1698 
1699 !        In case of PAW, add compensation charge contribution
1700          if (usepaw==1.and.extrap==1.and.neglect_pawhat==0) then
1701            call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,ispinor,ispinor,1,kg_diel,&
1702 &           lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1703 &           ngfftdiel,npwdiel,nspinor,ntypat,1,&
1704 &           pawang,pawtab,ph3d_diel,typat,dummy,wfraug,&
1705 &           mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1706            rhoextrap(:,:,:,ispinor)=rhoextrap(:,:,:,ispinor)+weight*wfraug(1,:,:,:)
1707          end if
1708 
1709 !        Performs the Fourier Transform of the density of the band,
1710 !        and store it in wfprod
1711          tim_fourwf=31
1712          call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1713 &         1,kg_diel,kg_diel,&
1714 &         mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,ndiel4,ndiel5,ndiel6,3,paral_kgb,tim_fourwf,weight,weight)
1715 !        In case of PAW, add compensation charge contribution if not already done
1716          if (usepaw==1.and.extrap==0.and.neglect_pawhat==0) then
1717            call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,iband,iband,1,1,1,kg_diel,&
1718 &           lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1719 &           ngfftdiel,npwdiel,nspinor,ntypat,0,&
1720 &           pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1721 &           mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1722          end if
1723 
1724 !        Perform now the summation of terms related to direct change of eigenvalues
1725 !        or extrapolation over higher bands
1726          wght1=zero ; wght2=zero
1727          if(occopt>=3 .and. testocc==1) wght1=doccde(iband+iband_shift+bdtot_index)*wtk(ikpt)/ucvol
1728          if(extrap==1) wght2=two*occ_deavg(iband+iband_shift)*wtk(ikpt)/ucvol
1729          weight=wght1+wght2
1730 
1731          if (abs(weight)>tol12) then
1732            do ipw2=1,npwdiel
1733 !            Only fills lower half of the matrix (here, the susceptibility matrix)
1734 !            Note that wfprod of the first index must behave like a density,
1735 !            so that it is used as generated by fourwf, while wfprod of the
1736 !            second index will be implicitely used to make a scalar product
1737 !            with a potential change, meaning that its complex conjugate must be
1738 !            used. This explains the following signs...
1739              do ipw1=ipw2,npwdiel
1740                susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1741 &               weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1742                susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1743 &               weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1744              end do
1745            end do
1746          end if
1747 
1748          if( occopt>=3 .and. testocc==1 .and. abs(wght1)>tol12) then
1749 !          Accumulate product of band densities by their doccde, for the
1750 !          computation of the effect of change of Fermi level.
1751            do ipw=1,npwdiel
1752              drhode(1,ipw,jsp)=drhode(1,ipw,jsp)+wfprod(1,ipw)*wght1
1753              drhode(2,ipw,jsp)=drhode(2,ipw,jsp)+wfprod(2,ipw)*wght1
1754            end do
1755 !          Also accumulate weighted sum of doccde
1756            sumdocc=sumdocc+wght1
1757          end if
1758 
1759 !        End condition of metallic occupancies or extrapolation
1760        end if
1761 
1762 !      End loop on spinorial components
1763      end do
1764 !    End loop on iband
1765    end do
1766 
1767    call timab(762,2,tsec)
1768 
1769 !  -- Compute now off-band-diagonal terms ------------------------------------
1770 !  -- Compute product of wavefunctions for different bands, inside the block -
1771 
1772    call timab(763,1,tsec)
1773 
1774 !  if (occopt<3) then
1775    tolocc=1.0d-3
1776 !  else
1777 !  tolocc=1.0d-8
1778 !  end if
1779 
1780    if(nbnd_current>1)then
1781      do ibd1=1,nbnd_current-1
1782        ibdshft1=ibd1+iband_shift
1783        do ibd2=ibd1+1,nbnd_current
1784          ibdshft2=ibd2+iband_shift
1785 
1786 !        If the occupation numbers are sufficiently different, or
1787 !        if extrapolation is used and the corresponding factor is not zero,
1788 !        then there is a contribution
1789          occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index)
1790          if( abs(occdiff)>tolocc      .or. &
1791 &         ( extrap==1 .and.            &
1792 &         ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) &
1793 &         ) then
1794 
1795            eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index)
1796 
1797 !          Loop on spinorial components
1798            do ispinor=1,nspinor
1799              jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1800 
1801 !            Store the contribution in wfraug
1802              do i3=1,ndiel3
1803                do i2=1,ndiel2
1804                  do i1=1,ndiel1
1805                    wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)&
1806 &                   +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2)
1807                    wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa1(1,i1,i2,i3,ispinor,ibd2)&
1808 &                   -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa1(2,i1,i2,i3,ispinor,ibd2)
1809                  end do
1810                end do
1811              end do
1812 
1813 !            Performs the Fourier Transform of the product, and store it in wfprod
1814              tim_fourwf=32
1815              call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1816 &             1,kg_diel,kg_diel, mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1817 &             ndiel4,ndiel5,ndiel6,3,paral_kgb,tim_fourwf,weight,weight)
1818 
1819 !            In case of PAW, add compensation charge contribution
1820              if (usepaw==1.and.neglect_pawhat==0) then
1821                call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibd2,ispinor,ispinor,1,kg_diel,&
1822 &               lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1823 &               ngfftdiel,npwdiel,nspinor,ntypat,0,&
1824 &               pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1825 &               mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1826              end if
1827 
1828 !            Perform now the summation
1829              wght1=zero ; wght2=zero
1830              if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol
1831              if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol
1832              weight=wght1+wght2
1833 
1834              if (abs(weight)>tol12) then
1835                do ipw2=1,npwdiel
1836 !                Only fills lower half of the matrix (here, the susceptibility matrix)
1837 !                Note that wfprod of the first index must behave like a density,
1838 !                so that it is used as generated by fourwf, while wfprod of the
1839 !                second index will be implicitely used to make a scalar product
1840 !                with a potential change, meaning that its complex conjugate must be
1841 !                used. This explains the following signs...
1842                  do ipw1=ipw2,npwdiel
1843                    susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1844 &                   weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1845                    susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1846 &                   weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1847                  end do
1848                end do
1849              end if
1850 
1851 !            End loop on spinorial components
1852            end do
1853 !          End condition of different occupation numbers or extrapolation
1854          end if
1855 !        End internal loop over bands
1856        end do
1857 !      End external loop over bands
1858      end do
1859 !    End condition of having more than one band
1860    end if
1861 
1862 !  Loop on secondary block, with fast varying index, in decreasing order.
1863    if(iblk1/=nblk)then
1864      do iblk2=nblk,iblk1+1,-1
1865        iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk
1866 
1867 !      Loop over bands to fft and store Fourier transform of wavefunction
1868        iband_shift2=nbnd_in_blk1+(iblk2-2)*nbnd_in_blk
1869        do iband=1,nbnd_in_blk
1870 !        Loop on spinorial components
1871          do ispinor=1,nspinor
1872            iwf=(ispinor-1)*npw_k+(iband+iband_shift2-1)*npw_k*nspinor+icg
1873 
1874 !          Obtain Fourier transform in fft box
1875            tim_fourwf=22
1876            cwavef(:,1:npw_k)=cg(:,1+iwf:npw_k+iwf)
1877            call fourwf(1,rhoaug,cwavef,dummy,wfraug,gbound,gbound,&
1878 &           istwf_k,kg_k,kg_k,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,npw_k,1,&
1879 &           ndiel4,ndiel5,ndiel6,0,paral_kgb,tim_fourwf,weight,weight)
1880            wfrspa2(:,:,:,:,ispinor,iband)=wfraug(:,:,:,:)
1881          end do
1882        end do
1883 
1884        do ibd1=1,nbnd_current
1885          ibdshft1=ibd1+iband_shift
1886          do ibd2=1,nbnd_in_blk
1887            ibdshft2=ibd2+iband_shift2
1888 
1889 !          If the occupation numbers are sufficiently different, or
1890 !          if extrapolation is used and the corresponding factor is not zero,
1891 !          then there is a contribution
1892            occdiff=occ(ibdshft1+bdtot_index)-occ(ibdshft2+bdtot_index)
1893            if( abs(occdiff)>tolocc      .or. &
1894 &           ( extrap==1 .and.            &
1895 &           ( abs(occ_deavg(ibdshft1)) + abs(occ_deavg(ibdshft2)) ) >tolocc ) &
1896 &           ) then
1897 
1898              eigdiff=eigen(ibdshft1+bdtot_index) - eigen(ibdshft2+bdtot_index)
1899 
1900 !            Loop on spinorial components
1901              do ispinor=1,nspinor
1902                jsp=isp+ispinor-1;if (nspden_eff==1) jsp=isp
1903 
1904 !              Store the contribution in wfraug
1905                do i3=1,ndiel3
1906                  do i2=1,ndiel2
1907                    do i1=1,ndiel1
1908                      wfraug(1,i1,i2,i3)=wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)&
1909 &                     +wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2)
1910                      wfraug(2,i1,i2,i3)=wfrspa1(2,i1,i2,i3,ispinor,ibd1)*wfrspa2(1,i1,i2,i3,ispinor,ibd2)&
1911 &                     -wfrspa1(1,i1,i2,i3,ispinor,ibd1)*wfrspa2(2,i1,i2,i3,ispinor,ibd2)
1912                    end do
1913                  end do
1914                end do
1915 
1916 !              Performs the Fourier Transform of the product, and store it in wfprod
1917                tim_fourwf=32
1918                call fourwf(1,rhoaug,dummy,wfprod,wfraug,gbound_diel,gbound_diel,&
1919 &               1,kg_diel,kg_diel,mgfftdiel,mpi_enreg_diel,1,ngfftdiel,1,npwdiel,&
1920 &               ndiel4,ndiel5,ndiel6,3,paral_kgb,tim_fourwf,weight,weight)
1921 
1922 !              In case of PAW, add compensation charge contribution
1923                if (usepaw==1.and.neglect_pawhat==0) then
1924                  call pawsushat(atindx,cprj_k,gbound_diel,gylmg_diel,ibd1,ibdshft2,ispinor,ispinor,1,kg_diel,&
1925 &                 lmax_diel,mgfftdiel,natom,nband_k,ndiel4,ndiel5,ndiel6,&
1926 &                 ngfftdiel,npwdiel,nspinor,ntypat,0,&
1927 &                 pawang,pawtab,ph3d_diel,typat,wfprod,dummy,&
1928 &                 mpi_atmtab=mpi_enreg%my_atmtab,comm_atom=mpi_enreg%comm_atom)
1929                end if
1930 
1931 !              Perform now the summation
1932                wght1=zero ; wght2=zero
1933                if(abs(occdiff)>tolocc) wght1= occdiff/eigdiff * two*wtk(ikpt)/ucvol
1934                if(extrap==1) wght2=(occ_deavg(ibdshft1)+occ_deavg(ibdshft2)) * two*wtk(ikpt)/ucvol
1935                weight=wght1+wght2
1936 
1937                if (abs(weight)>tol12) then
1938                  do ipw2=1,npwdiel
1939 !                  Only fills lower half of the matrix (here, the susceptibility matrix)
1940 !                  Note that wfprod of the first index must behave like a density,
1941 !                  so that it is used as generated by fourwf, while wfprod of the
1942 !                  second index will be implicitely used to make a scalar product
1943 !                  with a potential change, meaning that its complex conjugate must be
1944 !                  used. This explains the following signs...
1945                    do ipw1=ipw2,npwdiel
1946                      susmat(1,ipw1,jsp,ipw2,jsp)=susmat(1,ipw1,jsp,ipw2,jsp)+&
1947 &                     weight*(wfprod(1,ipw1)*wfprod(1,ipw2)+wfprod(2,ipw1)*wfprod(2,ipw2))
1948                      susmat(2,ipw1,jsp,ipw2,jsp)=susmat(2,ipw1,jsp,ipw2,jsp)+&
1949 &                     weight*(wfprod(2,ipw1)*wfprod(1,ipw2)-wfprod(1,ipw1)*wfprod(2,ipw2))
1950                    end do
1951                  end do
1952                end if
1953 
1954 !              End loop on spinorial components
1955              end do
1956 !            End condition of different occupation numbers or extrapolation
1957            end if
1958 !          End internal loop over bands
1959          end do
1960 !        End external loop over bands
1961        end do
1962 !      End loop on bloks
1963      end do
1964 
1965 !    Finish the loop on blok with iblk2=iblk1+1, so can use the
1966 !    FFTd wavefunctions for the next iblk1.
1967      do iband=1,nbnd_in_blk
1968        wfrspa1(:,:,:,:,1:nspinor,iband)=wfrspa2(:,:,:,:,1:nspinor,iband)
1969      end do
1970 
1971 !    End condition of iblk1/=nblk
1972    end if
1973 
1974    call timab(763,2,tsec)
1975 
1976 !  End loop on iblk1
1977  end do
1978 
1979 !DEBUG
1980 !write(std_out,*)' suskmm : exit '
1981 !do ipw1=1,npwdiel
1982 !write(std_out,*)ipw1,susmat(1,ipw1,1,ipw1,1),susmat(2,ipw1,1,ipw1,1)
1983 !end do
1984 !write(std_out,*)' suskmm : end of susmat '
1985 !stop
1986 !ENDDEBUG
1987 
1988  call destroy_mpi_enreg(mpi_enreg_diel)
1989  ABI_DEALLOCATE(cwavef)
1990  ABI_DEALLOCATE(dummy)
1991  ABI_DEALLOCATE(rhoaug)
1992  ABI_DEALLOCATE(wfprod)
1993  ABI_DEALLOCATE(wfraug)
1994  ABI_DEALLOCATE(wfrspa1)
1995  ABI_DEALLOCATE(wfrspa2)
1996 
1997  call timab(760,2,tsec)
1998 
1999 end subroutine suskmm