Logo Search packages:      
Sourcecode: abinit version File versions  Download package

csigme.F90

!{\src2tex{textfont=tt}}
!!****f* ABINIT/csigme
!! NAME
!! csigme
!!
!! FUNCTION
!! Calculating diagonal matrix elements or all matrix elements of self-energy operator
!!
!! COPYRIGHT
!! Copyright (C) 1999-2007 ABINIT group (FB, GMR, VO, LR, RWG, MG)
!! 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
!!  efermi=Fermi energy
!!  en(sp%nk,sp%nb,sp%nsppol)=KS energies for k-points, bands and spin
!!  ep=epsilonm1_parameters (see the definition of this structured datatype)
!!  er=epsilonm1_results (see the definition of this structured datatype)
!!  grottb(npwvec,2,nop)= contains the index of (IR) G where I is the identity or the inversion
!!   operation and G is one of the npwvec vectors in reciprocal space
!!  gvec(3,npwvec)=integer coordinates of each plane wave in reciprocal space
!!  igfft(npwvec,5,5,5)=index of G-G0 planewaves (see cigfft routine)
!!  lt_k datatype containing information on the little group
!!  jk=k point where to calculate the matrix element (might be redundant with sp)
!!  kbz(3,sp%nkbz)=k-point coordinates, full Brillouin zone
!!  kibz(3,sp%nk)=k-point coordinates, irreducible Brillouin zone
!!  ktab(sp%nkbz)= table giving for each k-point in the BZ (kBZ), the corresponding
!!   irreducible point (kIBZ), where kBZ= (IR) kIBZ and I is the inversion or the identity operation
!!  ktabi(sp%nkbz)= for each k-point in the BZ defines whether inversion has to be
!!   considered in the relation kBZ=(IR) kIBZ (1 => only R; -1 => -R)
!!  ktabr(nr,sp%nkbz)= index of R**-1 r in the FFT array where kBZ=(IR) kIBZ
!!  minbnd,maxbnd=indices of the lowest and highest bands for which the calculation
!!    of the matrix element is to be done
!!  mpi_enreg=informations about MPI parallelization (to be completed)
!!  ngfft1,ngfft1a,ngfft2,ngfft3=FFT grid dimensions
!!  nop=number of symmetry operations
!!  nr=number of points of FFT grid
!!  occ(sp%nk,sp%nb,sp%nsppol)=occupation numbers, for each k point in IBZ, and each band
!!  op(3,3,nop)=symmetry operations in reciprocal space
!!  q=q points at which matrix elements must be computed
!!  qbz(3,sp%nqbz)=coordinates of all q points in BZ
!!  qpg(sp%npwx,sp%nq)=|q+G| norm
!!  qtab(sp%nqbz)= table giving for each q-point in the BZ (qBZ), the corresponding
!!   irreducible point (qIBZ), where qBZ= (IR) qIBZ and I is the inversion or the identity operation
!!  qtabi(sp%nqbz)= for each q-point in the BZ defines whether inversion has to be
!!   considered in the relation qBZ=(IR) qIBZ (1 => only R; -1 => -R)
!!  qtabo(sp%nqbz)= the symmetry operation R that takes qIBZ to each qBZ
!!  sp=sigma_parameters (see the definition of this structured datatype)
!!  ucvol=unit cell volume
!!  wfr(nr,sp%nb,sp%nk,sp%nsppol)=wavefunctions in real space
!!
!! OUTPUT
!!  sr=sigma_results (see the definition of this structured datatype)
!!
!! NOTES
!!  The treatment of the divergence of Gygi+Baldereschi (PRB 1986) is included.
!!  The calculation of energy derivative is based on finite elements.
!!
!! TODO
!! One should eliminate the automatic arrays, and use allocatable ones
!!
!! PARENTS
!!      sigma
!!
!! CHILDREN
!!      calc_coh,calc_sig_noppm,calc_sig_ppm,cgemv,cggfft,diago_hamilt,distrb2
!!      findqg0,mpi_barrier,prep_coh,rho_tw_g,timab,xcomm_init,xmaster_init
!!      xme_init,xsum_mpi
!!
!! SOURCE

#if defined HAVE_CONFIG_H
#include "config.h"
#endif

subroutine csigme(sp,ep,jk,minbnd,maxbnd,&
& q,gvec,npwvec,op,nop,kibz,ktab,ktabr,ktabi,kbz,&
& ngfft1,ngfft1a,ngfft2,ngfft3,igfft,nr,en,occ,qbz,qtab,qtabi,&
& qtabo,qpg,ucvol,grottb,wfr,efermi,&
& er,sr,mpi_enreg,ppmodel,eigpot,b1,b2,b3,nbnds_per_proc,min_band_proc,max_band_proc,&
&  parallelism_is_on_kpoints,parallelism_is_on_bands,lt_k)

 use defs_basis
 use defs_datatypes

!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifdef HAVE_FORTRAN_INTERFACES
 use interfaces_01manage_mpi
 use interfaces_15gw, except_this_one => csigme
 use interfaces_lib01hidempi
#else
 use defs_xfuncmpi
#endif
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: jk,max_band_proc,maxbnd,min_band_proc,minbnd
 integer,intent(in) :: nbnds_per_proc,ngfft1,ngfft1a,ngfft2,ngfft3,nop,npwvec
 integer,intent(in) :: nr,ppmodel
 real(dp),intent(in) :: efermi,ucvol
 logical,intent(in) :: parallelism_is_on_bands,parallelism_is_on_kpoints
 type(MPI_type),intent(inout) :: mpi_enreg
 type(epsilonm1_parameters),intent(in) :: ep
 type(epsilonm1_results),intent(in) :: er
 type(little_group),intent(in) :: lt_k
 type(sigma_parameters),intent(in) :: sp
 type(sigma_results),intent(inout) :: sr
!arrays
 integer,intent(in) :: grottb(npwvec,2,nop),gvec(3,npwvec),ktab(sp%nkbz)
 integer,intent(in) :: ktabi(sp%nkbz),ktabr(nr,sp%nkbz),qtab(sp%nqbz)
 integer,intent(in) :: qtabi(sp%nqbz),qtabo(sp%nqbz)
 integer,target :: igfft(npwvec,5,5,5)
 real(dp),intent(in) :: b1(3),b2(3),b3(3)
 real(dp),intent(in) :: en(sp%nk,min_band_proc:max_band_proc,sp%nsppol)
 real(dp),intent(in) :: kbz(3,sp%nkbz),kibz(3,sp%nk)
 real(dp),intent(in) :: occ(sp%nk,min_band_proc:max_band_proc,sp%nsppol)
 real(dp),intent(in) :: op(3,3,nop),q(3,sp%nq),qbz(3,sp%nqbz)
 real(dp),intent(in) :: qpg(sp%npwx,sp%nq)
 complex,intent(in) :: eigpot(sp%npwc,sp%npwc,sp%nq)
 complex,intent(in) :: wfr(nr,min_band_proc:max_band_proc,sp%nk,sp%nsppol)

!Local variables ------------------------------
!(To be ordered)
!Warning : there are automatic arrays
! bz_geometry_factor: (see gwa.pdf, appendix A.4)
! sphere=7.79, fcc=7.44, sc=6.188, bcc=6.946, wz=5.255 (Lu Fu-Fa)
! grid on the omega real axis
!For the calculation of the kinetic and
!Coulomb contributions to the bandgap energy.
!FBruneval: added variables
! real(dp),parameter :: tol_empty=0.01 !below this value the state is assumed empty
!scalars
 integer :: conv,ff,i,ib,ierr,ig,ig1,iggp,igp,ii,ii1,iik,iiq,ik,ikbz,ikibz,io
 integer :: ioe0j,iop,iopk,iopq,iq,iq0bz,iqbz,iqibz,is,isless,j,jb,jb_proc_rank
 integer :: jik,jj,jkbz,jkibz,k_index1,k_index2,kb,master,mband,me,me_loc
 integer :: nbdblock,nkpt,nomega,npwc1,npwc2,nrb,nsize,spaceComm,theta_fact
 integer :: tim_fourdp,twofm1,zb
 real(dp),parameter :: bz_geometry_factor=7.44,domegareal=0.5/Ha_eV
 real(dp) :: a,b,den,den_coh,den_sex,domegasi,e0i,e0j,fact_sp,i_sz,imdsigmee0
 real(dp) :: inv_den,kincontrib,omegasi,otw,q0vol,redsigmee0,reomegame0i,smrt
 real(dp) :: sqrt_i_sz,theta_mu_minus_e0i,tol_empty,twofm1_zcut,z2
 complex :: ccoh,csex,ct,dct,dsigc,full_den,idelta,ieta,inv_full_den,num,numf
 complex :: rhotwgdp_igp,sigc,sigcohme,sigxme,twofm1_idelta,z0,zz
 complex,external :: cdotc
 logical :: am_master,cohsex,kinclb
 character(len=500) :: message
!arrays
 integer :: g0(3),g01(3),irottb(nr,nop)
 integer,allocatable :: nband_t(:)
 integer,pointer :: igfftg0(:)
 real(dp) :: e0pde(sp%nomegasrd,sp%nb,sp%nsppol),ki(3),kj(3),kjmki(3),kk(3)
 real(dp) :: qbzpg(sp%npwx),scme(sp%nomegasrd),tsec(2)
 real(dp),allocatable :: omegame0i(:),otq(:,:)
 complex :: rhotwg(sp%npwx),rhotwg_ki(sp%npwx,minbnd:maxbnd)
 complex :: rhotwgdpcc(sp%npwx),rhotwgp(sp%npwx),sigsex(sp%npwc)
 complex,allocatable :: botsq(:,:),eig(:,:),eig_1(:),epsm1q(:,:,:)
 complex,allocatable :: epsm1q_trcc(:,:,:),ket(:,:),ket1(:,:),ket2(:,:)
 complex,allocatable :: rhotwg1(:),rhotwg_1(:,:,:),rhotwgdpccc(:),sigccoh(:,:)
 complex,allocatable :: sigckin(:),sigcme2(:,:),sigcme_3(:),sigcme_new(:)
 complex,allocatable :: sigcsex(:,:),sigctmp(:),sigma_phi(:,:),sum(:)
 complex,allocatable :: wfr_jb(:,:,:),wfr_temp(:,:)
!no_abirules
#if defined MPI || defined MPI_FFT
           complex,allocatable::xsum1(:,:,:),xsum2(:,:)
           real(dp),allocatable::xsum3(:,:)
#endif
!MG060926 added (or modified) variables
 integer :: isppol
!added spin dimension
 complex, allocatable :: wfg2(:,:,:,:)
 complex,allocatable :: sigcme_tmp(:,:,:,:),sigxme_tmp(:,:,:)
!END MG
!MG for symmetrization
integer :: wtqp,wtqm,nqbzk,indx,ncplx(sp%nsppol),wcplx,icx,sumcxtab
integer :: cplx(maxbnd-minbnd+1,sp%nsppol)
integer,allocatable :: cxtab(:,:,:)
real(dp) :: ek0
complex,allocatable :: dummy_xme(:,:,:),dummy_cme(:,:,:,:)
integer, allocatable :: ndwn(:,:),nup(:,:)
!ENDMG
!

! *************************************************************************

#ifdef VMS
!DEC$ ATTRIBUTES ALIAS:'CDOTC' :: cdotc
!DEC$ ATTRIBUTES ALIAS:'CGEMV' :: cgemv
#endif

!DEBUG
! write(6,*)' csigme : enter ',minbnd,maxbnd,sp%nb
!ENDDEBUG

!These lines added by Shaltaf for parallelization 10/08/05
!parallization is on k points only, one must include the bands
!too in the future
!Init mpi_comm
!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
#ifndef HAVE_FORTRAN_INTERFACES
 real(dp) :: linfit
#endif
!End of the abilint section

 if(sr%nomega/=sp%nomegasr) then
  write(message,'(4a)')ch10,&
& ' csigme : BUG- ',ch10,&
& ' sr%nomega/=sp%nomegasr '
  call wrtout(6,message,'COLL')
  call leave_new('COLL')
 end if


 call xcomm_init(mpi_enreg,spaceComm) !Init mpi_comm
 call xme_init(mpi_enreg,me)          !Init me
 call xmaster_init(mpi_enreg,master)  !Init master
 !End of parallelization initialisation

 if(mod(sp%gwcalctyp,10)==0) write(message,'(a)')' standard GW with PPM'
 if(mod(sp%gwcalctyp,10)==2) write(message,'(a)')' standard GW without PPM'
 if(mod(sp%gwcalctyp,10)==5) write(message,'(a)')' Hartree-Fock calculation'
 if(mod(sp%gwcalctyp,10)==6) write(message,'(a)')' Screened Exchange calculation'
 if(mod(sp%gwcalctyp,10)==7) write(message,'(a)')' COHSEX calculation'
 if(mod(sp%gwcalctyp,10)==8) write(message,'(a)')' model GW with PPM'
 if(mod(sp%gwcalctyp,10)==9) write(message,'(a)')' model GW without PPM'
 call wrtout(6,message,'COLL')

 if(sp%gwcalctyp<10) then
  write(message,'(a)')' Perturbative Calculation'
 else if(sp%gwcalctyp<20) then
  write(message,'(a)')' Self-Consistent on Energies only'
 else
  write(message,'(a)')' Self-Consistent on Energies and Wavefunctions'
 end if
 call wrtout(6,message,'COLL')

 call timab(421,1,tsec)
 call timab(422,1,tsec)

 !These lines added by GMR on 04/05/05 for the calculation
 !of the decomposition of the correlation part of sigma into
 !Coulomb-hole (coh) and screened-exchange (sex)
 cohsex=.false.
 if(mod(sp%splitsigc,2)==1) cohsex=.true.
 !Kinetic (kin) and Coulomb (clb)
 kinclb=.false.
 if(int(sp%splitsigc/2)==1) kinclb=.true.

 tim_fourdp=2
 idelta=cmplx(0,sp%zcut)

 !analytic integration of q**-2 over the volume element:
 ! $\int_V d^3q 1/q^2 = bz_geometric_factor V^(-2/3)$
 q0vol=(8* pi*pi*pi)/(ucvol*sp%nkbz) ! V_BZ / N_k
 i_sz = bz_geometry_factor * q0vol**(-two_thirds)
 sqrt_i_sz = sqrt(i_sz)

 nrb=ngfft1*ngfft2*ngfft3 !FFT factor

 !print info
  write(message,'(3a,3f8.3,2a,i3,a,i3,a)')&
& ' calculating <nk|sigma|nk>',ch10,      &
& ' k = ',(kbz(i,jk),i=1,3),ch10,         &
& ' bands n = from ',minbnd,' to ',maxbnd,ch10
  call wrtout(6,message,'COLL')

 jkbz=jk
 jkibz=ktab(jkbz)
 jik=(3-ktabi(jkbz))/2
 kj(:)=kbz(:,jkbz)

 !MG added nsspol as last dimension
 !Diagonal elements of $\Sigma_x$ and $\Sigma_c(\omega)$
 sr%sigxme(:,jkibz,:)=0
 sr%sigcme(:,jkibz,:,:)=0

 !MG060914  since we have to evaluate \Sigma around the KS\QP eigenvalue that can
 !depend on the spin, a the spin index is needed
 !anyway I do not understand why we are storing the frequencies for ALL the bands, it is a waste
!MG NOTE  it is better to force an even sp%nomegasrd so we have the KS\QP
!eigenvalue in the middle of the interval

! Rshaltaf:
!1) Not every proc may have the wavefunctions that correspond to the block of bands (minbamd:maxband)
!2) every Proc has to verify that it contains the correspondent bands, if not
!3) it must communicate with a proc that does have them
 allocate(wfr_jb(nr,minbnd:maxbnd,sp%nsppol))
 do jb=minbnd,maxbnd
  if(mpi_enreg%nproc>1.and.parallelism_is_on_bands)then
   jb_proc_rank=minval(abs(mpi_enreg%proc_distrb(jb,:,:)))
   if(me==jb_proc_rank)then
    wfr_jb(:,jb,:)=wfr(:,jb,jkibz,:)
   end if
   call xcast_mpi(wfr_jb,jb_proc_rank,spaceComm,ierr)
  else
   wfr_jb(:,jb,:)=wfr(:,jb,jkibz,:) ! kpoint case
  end if
 end do

 ioe0j=sp%nomegasrd/2+1
 do is=1,sp%nsppol
  do jb=1,sp%nb
   do io=1,sp%nomegasrd
    e0pde(io,jb,is)=sr%egw(jb,jkibz,is)+sp%deltae*(io-ioe0j)
    sr%omegasrd(jb,jkibz,io,is) = e0pde(io,jb,is)
   end do
  end do
 end do

 nomega=sp%nomegasr+sp%nomegasrd
 !sigcme2 is the array used to accumulate diagonal matrix elements over k-points and ALL the bands!
 allocate(sigcme2(nomega,sp%nb),sigcme_3(nomega))
 allocate(omegame0i(nomega))
 !these are the arrays used to accumulate the matrix elements over k-points and bands
 allocate(sigcme_tmp(nomega,sp%nb,sp%nb,sp%nsppol),sigxme_tmp(sp%nb,sp%nb,sp%nsppol))

 allocate(ket(sp%npwc,nomega),sigctmp(nomega))
 allocate(ket1(sp%npwc,nomega),ket2(sp%npwc,nomega))

 !zero matrix elements
 sigcme2(:,:)=zero
 !MG added spin variable
 sigxme_tmp(:,:,:)=zero
 sigcme_tmp(:,:,:,:)=zero

 !These lines added by GMR on 04/05/05 for the calculation
 !of the decomposition of the correlation part of sigma into
 !Coulomb-hole (coh) and screened-exchange (sex)
 !MG added new dimension for spin
 if (cohsex) then
  allocate(sigccoh(sp%nb,sp%nsppol))
  allocate(sigcsex(sp%nb,sp%nsppol))
  sigccoh(:,:)=zero !zero matrix elements
  sigcsex(:,:)=zero !zero matrix elements
 end if

 !These lines added by YMN on 07/04/04 for the calculation
 !of the kinetic and Coulomb contributions to the bandgap energy.
!MG FIXME this array is not used anymore
 if (kinclb) then
  allocate(sigckin(sp%nb))
  sigckin(:) = 0._dp
 end if

!MG Added for symmetrization
 if (lt_k%sym_flag/=0) then

  nqbzk=0
  do ii=1,lt_k%nkbz
   nqbzk=nqbzk+lt_k%ibzq(ii)
  end do

  !evaluate number of complexes and number of bands in each complex
  allocate (cxtab(maxbnd-minbnd+1,sp%nsppol,maxbnd-minbnd+1))
  cxtab=0

  do isppol=1,sp%nsppol
   do ib=1,maxbnd-minbnd+1
    do jb=1,maxbnd-minbnd+1
     !the tolerance is a little bit arbitrary (~0.001 eV)
     !could be reduced, in particular in the proximity of accidental degeneracies
     if (abs(en(jkibz,ib-1+minbnd,isppol)-en(jkibz,jb-1+minbnd,isppol))<0.00005) then
      cxtab(ib,isppol,jb)=1
     end if
    end do
   end do
  end do ! isppol
  !write(*,*)cxtab
 end if
!ENDMG

 call timab(422,2,tsec)

 if (lt_k%sym_flag/=0) then
  write(message,'(a,i6,a)')' calculation status (',nqbzk,' to be completed):'
 else
  write(message,'(a,i6,a)')' calculation status (',sp%nkbz,' to be completed):'
 end if
 call wrtout(6,message,'COLL')

#if defined MPI
           call MPI_BARRIER(spaceComm,ierr)
#endif

 allocate(eig(sp%npwc,sp%npwc))

 !Prepare the static COHSEX calculation:
 !calculate the Fourier transform of \phi_j^*(r)\phi_k(r)
 if(mod(sp%gwcalctyp,10)==7) then
  !MG060926 now wfg2 is allocated only in case of COHSEX, added spin dimension
  allocate (wfg2(nr,minbnd:maxbnd,minbnd:maxbnd,sp%nsppol))
  wfg2=zero
  !MG071009 FIXME I think there is a BUG in prep_coh because we are not taking into
  !               account the time reversal symmetry.
  !However now prep_coh also works in the case nsppol
  call prep_coh(sp%nb,nr,sp%nsppol,wfr(:,:,jkibz,:),ktabr(:,jkibz),minbnd,maxbnd,&
&                ngfft1,ngfft2,ngfft3,tim_fourdp,wfg2,mpi_enreg)
 end if

!allocate(rhotwg_1(sp%npwx,minbnd:maxbnd,sp%nb))

 if(ppmodel==1.or.ppmodel==2)then
  npwc1=sp%npwc;npwc2=sp%npwc
 else if(ppmodel==4)then
  npwc1=sp%npwc;npwc2=1
 else ! for ppmodel3 and noppmodel
  npwc1=1;npwc2=1
 end if

 allocate(botsq(sp%npwc,npwc1),otq(sp%npwc,npwc2))
 if( mod(sp%gwcalctyp,10)==2 .or. mod(sp%gwcalctyp,10)==6 .or. &
&    mod(sp%gwcalctyp,10)==7 .or. mod(sp%gwcalctyp,10)==9     ) then
  allocate(epsm1q(sp%npwc,sp%npwc,ep%nomega),epsm1q_trcc(sp%npwc,sp%npwc,ep%nomega) )
 end if

 !MG this is for the correct normalization of theta_mu_minus_e0,
 !   since in case of nsppol==2 occ $\in [0,1]$
 if (sp%nsppol==1) then
  fact_sp=0.5
  tol_empty=0.01 !below this value the state is assumed empty
 else if (sp%nsppol==2) then
  fact_sp=1.0
  tol_empty=0.005 !to be consistent and obtain similar results if a metallic
                  !spin unpolarized system is  treated using nsppol==2
 else
  write(message,'(3a)')&
&  ' csigme : BUG ',ch10,&
&  ' nsppol must be 1 or 2'
  call wrtout(6,message,'COLL')
  call leave_new('COLL')
 end if

!DEBUG
!write(6,*)' csigme : lt_k%sym_flag=',lt_k%sym_flag
!write(6,*)' csigme : sp%nb ',sp%nb
!ENDDEBUG

 do ikbz=1,sp%nkbz !loop over k_i in BZ

  !This line added by Shaltaf for parallelization
  if(parallelism_is_on_kpoints)then
   if(minval(abs(mpi_enreg%proc_distrb(ikbz,:,:)-mpi_enreg%me))/=0) cycle
  end if
  write(6,*)ikbz
  call timab(423,1,tsec)

  ikibz=ktab(ikbz)
  iik=(3-ktabi(ikbz))/2
  ki(:)=kbz(:,ikbz)

  kjmki(:)=kj(:)-ki(:)
  call findqg0(iq,g0,kjmki,sp%nqbz,qbz) !identify q and G0 where q+G0=k_j-k_i

!MG SYMMETRIZATION here there was a problem with gfort
! the following coding is safer
  if (lt_k%sym_flag/=0) then
   if (lt_k%ibzq(iq)/=1) cycle !only q in the IBZ_k
   wtqp=0
   wtqm=0
   do jj=1,lt_k%nop
    wtqp=wtqp+lt_k%wtksym(1,jj,iq)
    wtqm=wtqm+lt_k%wtksym(2,jj,iq)
   end do
  end if
!END SYMMETRIZATION

  write(message,'(a,2i6)')' csigme : ik=',ikbz
  call wrtout(06,message,'COLL')

  ikibz=ktab(ikbz)
  iik=(3-ktabi(ikbz))/2
  ki(:)=kbz(:,ikbz)

  kjmki(:)=kj(:)-ki(:)
  call findqg0(iq,g0,kjmki,sp%nqbz,qbz) ! identify q and G0 where q+G0=k_j-k_i
  iqibz=qtab(iq)
  iopq=qtabo(iq)
  iiq=(3-qtabi(iq))/2
  igfftg0=>igfft(:,g0(1)+3,g0(2)+3,g0(3)+3)

  if(mod(sp%gwcalctyp,10)==0.or.mod(sp%gwcalctyp,10)==8) then

!  calculate plasmonpole parameters for this qBZ (converting from qIBZ to qBZ)
!  MG NOTE here we are not taking into account a possible umklapp vector!!!!
!  The code however stops before in surot.F90
   if(ppmodel==1.or.ppmodel==2)then
    do j=1,sp%npwc
     do i=1,sp%npwc
      botsq(grottb(i,iiq,iopq),grottb(j,iiq,iopq))=er%bigomegatwsq(i,j,iqibz)
      otq(grottb(i,iiq,iopq),grottb(j,iiq,iopq))=er%omegatw(i,j,iqibz)
     end do
    end do
   else if(ppmodel==3)then
    do ii=1,sp%npwc         !! ii: DM bands
     otq(ii,1)=er%omegatw(ii,1,iqibz)
     botsq(ii,1)=er%bigomegatwsq(ii,1,iqibz)
     do j=1,sp%npwc
      eig(grottb(j,iiq,iopq),ii)=eigpot(j,ii,iqibz)
     end do
    end do
   else if(ppmodel==4)then
    do ii=1,sp%npwc     !! ii: DM bands
     otq(ii,1)=er%omegatw(ii,1,iqibz)
     do j=1,sp%npwc
      botsq(grottb(j,iiq,iopq),ii)=er%bigomegatwsq(j,ii,iqibz) ! ii is the band index of DM
     end do
    end do
   end if ! ppmodel

   if(iiq==2)then
    botsq(:,:)=conjg(botsq(:,:))
    eig(:,:)=conjg(eig(:,:))
   end if

  else if(mod(sp%gwcalctyp,10)==2.or.mod(sp%gwcalctyp,10)==9) then

   !numerical integration or model GW using the contour deformation technique
   do io=1,ep%nomega
    do i=1,sp%npwc
     do j=1,sp%npwc
      epsm1q(grottb(i,iiq,iopq),grottb(j,iiq,iopq),io)=er%epsm1(i,j,io,iqibz)
     end do
    end do
   end do
   if(iiq==2) epsm1q(:,:,:)=conjg(epsm1q(:,:,:))
   !minus exchange contribution
   do io=1,ep%nomega
    do i=1,sp%npwc
     epsm1q(i,i,io)=epsm1q(i,i,io)-1.
    end do
   end do

  else if(mod(sp%gwcalctyp,10)==6.or.mod(sp%gwcalctyp,10)==7) then

   !SEX or COHSEX, only \omega=0 is needed here
   do i=1,sp%npwc
    do j=1,sp%npwc
     epsm1q(grottb(i,iiq,iopq),grottb(j,iiq,iopq),1)=er%epsm1(i,j,1,iqibz)
    end do
   end do
   if(iiq==2) epsm1q(:,:,1)=conjg(epsm1q(:,:,1))
   !minus exchange contribution
   do i=1,sp%npwc
    epsm1q(i,i,1)=epsm1q(i,i,1)-1.
   end do

  else if(mod(sp%gwcalctyp,10)==9) then

   !for model GW create transpose(conjg(epsm1q))
   do io=1,ep%nomega
    epsm1q_trcc(:,:,io)=transpose(conjg(epsm1q(:,:,io)))
   end do

  end if ! gwcalctyp

! Set up table of |q+G| in BZ
! MG NOTE Also in this case we are not taking into account umklapp processes
! The code however stops before in surot.F90
  do ig=1,sp%npwx
   qbzpg(grottb(ig,iiq,iopq))=qpg(ig,iqibz)
  end do

  call timab(423,2,tsec)

! the following lines commented by Rshaltaf and moved down for
! saving memory rhotwg_k(npwx,min:max,nb) is a huge matrix
! dont you think that Fabien??
! FB061013: This is now fixed with the small array rhotwg_ki(npwx,min:max)

  do ib=1,sp%nb

   if(parallelism_is_on_bands)then
    if(minval(abs(mpi_enreg%proc_distrb(ib,:,:)-mpi_enreg%me))/=0) cycle
   end if
!  MG060926 added inner loop over spin

   do isppol=1,sp%nsppol
!   MG060926 NOTE  if nsppol==2 ==>  0<=occ<= 1,
!   tol_empty is set fo 0.1 if nsspol==1 or 0.05 if nsppol==2
!   in order to be consistent and obtain the same results if a metallic
!   spin unpolarized system is treated using nsppol==2

!   Skip empty state ib for HF, SEX, and COHSEX
    if ( occ(ikibz,ib,isppol)<tol_empty.and.&
&      (mod(sp%gwcalctyp,10)==5.or.mod(sp%gwcalctyp,10)==6 .or. mod(sp%gwcalctyp,10)==7) ) cycle

    do jb=minbnd,maxbnd
!    Calculate all rho-twiddle(G) at once
     call rho_tw_g(sp%npwx,nr,nrb,ngfft1a,ngfft1,ngfft2,ngfft3,igfftg0,&
&     wfr(:,ib,ikibz,isppol),iik,ktabr(:,ikbz),wfr_jb(:,jb,isppol),jik,ktabr(:,jkbz),&
&     rhotwg_ki(:,jb),tim_fourdp)
!    call rho_tw_g(sp%npwx,nr,nrb,ngfft1a,ngfft1,ngfft2,ngfft3,igfftg0,wfr(:,ib,ikibz,isppol),&
!&    iik,ktabr(:,ikbz),wfr(:,kb,jkibz,isppol),jik,ktabr(:,jkbz),rhotwgp(:),tim_fourdp)
!    Divide rho-twiddle(G) by |q+G|
     rhotwg_ki(:,jb)=rhotwg_ki(:,jb)/qbzpg(:)

!    Treat analytically the case q->0
     if(ikbz==jkbz) then
      if(ib==jb) then
       rhotwg_ki(1,jb)=sqrt_i_sz
      else
       rhotwg_ki(1,jb)=0 !  ignore always
      end if
     end if
    end do !jb

    !MG using coefficient fact_sp defined outside the loops
    theta_mu_minus_e0i=fact_sp*occ(ikibz,ib,isppol)

    !MG starting point to evaluate the derivative of sigma and the spectral function, it depends on spin
    e0i=en(ikibz,ib,isppol)

    !MG this is for the spectral function, now we use e0i=en(ikibz,ib,isppol)
    !FIXME the interval is not centered on eoi !!!!!!!!!!!! WHY?
    if(sp%nomegasr>0) omegame0i(1:sr%nomega)=real(sr%omega(1:sr%nomega))-e0i

    do kb=minbnd,maxbnd

     call timab(424,1,tsec)

     do io=sp%nomegasr+1,nomega
      !MG  calculating frequencies $\omega$ - $\epsilon_in$ to evaluate the derivative of $\Sigma$
      !    here we have to take into account that the KS energy could depend on spin
      !    NOTE it is crazy because we have stored e_KS+ Delta \omega in sr%omegasrd and now
      !    we have to subtract e_KS
      omegame0i(io) = real(sr%omegasrd(kb,jkibz,io-sp%nomegasr,isppol)) - e0i
     end do ! io

     rhotwgp(:)=rhotwg_ki(:,kb)

     !First calculate the ket  \Sigma | \phi_k >
     ket(:,:)=czero
     ket1(:,:)=czero
     ket2(:,:)=czero

     if(mod(sp%gwcalctyp,10)==0) then !GW calculation WITH Plasmon-Pole Model

      call calc_sig_ppm(sp%npwc,nomega,rhotwgp,botsq,otq,omegame0i,sp%zcut,&
&      theta_mu_minus_e0i,ket,ppmodel,eig,sigcme_3,sp%npwx,npwc1,npwc2)
      if(ppmodel==3.or.ppmodel==4)then
!      MG APPLYING SYMMETRIES, I am not completely sure but the automatic tests are
!      correct FIXME
       if (lt_k%sym_flag==0) then
        sigcme2(:,kb)=sigcme2(:,kb)+sigcme_3(:) !only relevant if ppm3 or ppm4
       else
!       only relevant if ppm3 or ppm4
        sigcme2(:,kb)=sigcme2(:,kb)+&
&        (wtqp+wtqm)*real(sigcme_3(:))+(wtqp-wtqm)*(0.,1.)*aimag(sigcme_3(:))
       end if
      end if

     else if(mod(sp%gwcalctyp,10)==2) then !GW calculation numerical integration with contour deformation

      call calc_sig_noppm(sp%npwc,nomega,ep%nomega,ep%nomegaer,ep%nomegaei,rhotwgp(1:sp%npwc),&
&      er%omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket)

     else if(mod(sp%gwcalctyp,10)==6.or.mod(sp%gwcalctyp,10)==7) then
      !for SEX or COHSEX, here only the SEX part is considered
      call cgemv('N',sp%npwc,sp%npwc,(1.,0.),epsm1q(:,:,1),sp%npwc,rhotwgp(1:sp%npwc),1,(0.,0.),sigsex,1)
      sigsex(:)=-theta_mu_minus_e0i*sigsex(:)
      !static SEX is energy independent
      do io=1,nomega
       ket(:,io)=sigsex(:)
      end do

     else if(mod(sp%gwcalctyp,10)==8) then !MODEL GW calculation WITH Plasmon-Pole Model

      !calculate \Sigma(E_k) | k > to obtain < j | \Sigma(E_k) | k >
      allocate(sigcme_new(nomega))
      call calc_sig_ppm(sp%npwc,nomega,rhotwgp(1:sp%npwc),botsq,otq,omegame0i,sp%zcut,&
&      theta_mu_minus_e0i,ket1,ppmodel,eig,sigcme_new,sp%npwx,npwc1,npwc2)
      !calculate < k | \Sigma(E_k) to obtain < k | \Sigma(E_k) | j >^*
      if(sp%gwcalctyp==28) then
       call calc_sig_ppm(sp%npwc,nomega,rhotwgp(1:sp%npwc),conjg(transpose(botsq(:,:))),&
&       transpose(otq(:,:)),omegame0i,sp%zcut,theta_mu_minus_e0i,ket2,&
&       ppmodel,eig,sigcme_3,sp%npwx,npwc1,npwc2)
       ket(:,:)=(ket1(:,:)+ket2(:,:))*0.5
      else
       ket(:,:)=ket1(:,:)
      end if
      deallocate(sigcme_new)

     else if(mod(sp%gwcalctyp,10)==9) then !MODEL GW calculation numerical integration

      !calculate \Sigma(E_k) | k > to obtain < j | \Sigma(E_k) | k >
      call calc_sig_noppm(sp%npwc,nomega,ep%nomega,ep%nomegaer,ep%nomegaei,rhotwgp(1:sp%npwc),&
&                       er%omega,epsm1q,omegame0i,theta_mu_minus_e0i,ket1)
      !calculate < k | \Sigma(E_k) to obtain < k | \Sigma(E_k) | j >^*

      if(sp%gwcalctyp==29) then
       call calc_sig_noppm(sp%npwc,nomega,ep%nomega,ep%nomegaer,ep%nomegaei,rhotwgp(1:sp%npwc),&
&                        er%omega,epsm1q_trcc,omegame0i,theta_mu_minus_e0i,ket2)
       ket(:,:)=(ket1(:,:)+ket2(:,:))*0.5
      else
       ket(:,:)=ket1(:,:)
      end if

     end if !of if (sp%gwcalctyp...)

     !FBruneval 05/12/15
     !we could use hermitianity for HF, SEX, COHSEX to reduce next loop to:   do jb=minbnd,kb
     do jb=minbnd,maxbnd

      !check self-consistent or not? if not ==> only diagonal elements
      if(sp%gwcalctyp<20.and.jb/=kb) cycle
      rhotwg(:)=rhotwg_ki(:,jb)

      !Second calculate the bra ket < \phi_j | \Sigma | \phi_k >
      do io=1,nomega
       sigctmp(io) = cdotc(sp%npwc,rhotwg(1:sp%npwc),1,ket(:,io),1)
      end do

!DEBUG
!if(ikbz==1)then
!       write(6,*)' write ket =',ket(:,1)
!       write(6,*)' write ket1 =',ket1(:,1)
!       write(6,*)' write ket2 =',ket2(:,1)
!       write(6,*)' otq=',otq(1,1)
!       write(6,*)' botsq=',botsq(1,1)
!endif
!ENDDEBUG

      !Set up bare exchange matrix element
      !must be multiplied later by 4*PI/(ucvol*nkbz)
      sigxme=-cdotc(sp%npwx,rhotwg,1,rhotwgp,1)*theta_mu_minus_e0i
      !MG060926 added spin index

!MG SYMMETRIZING
      if (lt_k%sym_flag==0) then
       sigxme_tmp(jb,kb,isppol) = sigxme_tmp(jb,kb,isppol) + sigxme
      else
       sigxme_tmp(jb,kb,isppol) = sigxme_tmp(jb,kb,isppol) + &
&       (wtqp+wtqm)*real(sigxme)+(wtqp-wtqm)*(0.,1.)*aimag(sigxme)
      end if
!END MG

      call timab(424,2,tsec)
      call timab(425,1,tsec)

      !the static COH is calculated only once, since it is independent with respect to
      !the band index ib used to sum over bands
      if(mod(sp%gwcalctyp,10)==7.and.ib==1) then
       !MG060926 added spin
       call calc_coh(nr,sp%npwx,sp%npwc,ngfft1a,ngfft1,ngfft2,ngfft3,&
&                  gvec,epsm1q(:,:,1),qbzpg,i_sz,iqibz,jb,kb,wfg2(:,jb,kb,isppol),sigcohme)
       !static COH is energy independent
       !MG FIXME Check this part, I cannot understand why we are accumulating in sigctmp
       sigctmp(:)=sigctmp(:)+sigcohme
      end if

      !MG060926 added spin to accumulate
      !must be multiplied later by 4*PI/(ucvol*nkbz)
      !SYMMETRIZING
      if (lt_k%sym_flag==0) then
       sigcme_tmp(:,jb,kb,isppol) = sigcme_tmp(:,jb,kb,isppol) + sigctmp(:)
      else
       sigcme_tmp(:,jb,kb,isppol) = sigcme_tmp(:,jb,kb,isppol) + &
&      (wtqp+wtqm)*real(sigctmp(:))+(wtqp-wtqm)*(0.,1.)*aimag(sigctmp(:))
      end if

      !These lines added by GMR on 04/05/05 for the calculation
      !of the decomposition of the correlation part of sigma into
      !Coulomb-hole (coh) and screened-exchange (sex)
      !MG FIXME the case nsppol==2 is not well tested
      !   Do not understand why > 0.5
      if(theta_mu_minus_e0i>0.5) then
       ff=1
      else
       ff=0
      end if

      twofm1=2*ff-1
      twofm1_zcut=twofm1*sp%zcut
      twofm1_idelta=twofm1*idelta
      !Now, introduce rhotwgdpcc, for speed reasons
      rhotwgdpcc(:)=conjg(rhotwg(:))
      reomegame0i=omegame0i(ioe0j)

      if ((cohsex).and.(io == sr%nomega+ioe0j)) then
       ccoh=0
       csex=0
       do igp = 1, sp%npwc
        rhotwgdp_igp=rhotwg(igp)
        do ig = 1, sp%npwc
         if(ppmodel==3.and.ig/=igp)cycle
         otw=dble(otq(ig,igp)) ! in principle otw -> otw - ieta
         num = rhotwgdpcc(ig)*botsq(ig,igp)*rhotwgdp_igp
         den_coh = reomegame0i-otw
         if(den_coh**2>sp%zcut**2)then
          ccoh = ccoh + num / (den_coh * otw)
         else ! if den is small
          ccoh = ccoh + num*cmplx(den_coh,twofm1_zcut)/ &
&                      ((den_coh**2+twofm1_zcut**2)*otw)
         end if
         den_sex = reomegame0i**2-otw**2
         if(den_sex**2>sp%zcut**2)then
          csex = csex - 2 * ff * num / den_sex
         else ! if den is small
          csex = csex - 2*ff*num*cmplx(den_sex,twofm1_zcut)/ &
&                      ((den_sex**2+twofm1_zcut**2)*otw)
         end if
        end do ! ig
       end do ! igp

       !Must be multiplied later by 4*pi/(ucvol*nkbz).
       sigccoh(jb,isppol) = sigccoh(jb,isppol) + ccoh / two
       sigcsex(jb,isppol) = sigcsex(jb,isppol) + csex / two

      end if ! cohsex

!     Added by YMN on 07/04/04.
!     Calculation of the kinetic contribution to the bandgap energy.

!     if ((kinetic).and.(io == sr%nomega+ioe0j)) then
!      print *,' omega = ',sr%omegasrd(jb,jkibz,io-sp%nomegasr)
!      do igp = 1, sp%npwc
!       do ig = 1, sp%npwc
!        dct = zero
!        do ig1 = 1, sp%npwc
!         call fkin(ff,kincontrib,reomegame0i,otq(ig,ig1),otq(ig1,igp),sp%zcut)
!         dct = dct+botsq(ig,ig1)*botsq(ig1,igp)*kincontrib
!        end do
!        ct = ct-rhotwgdpcc(ig)*dct*rhotwg(igp)
!       end do
!      end do
!      sigckin(jb) = sigckin(jb)-ct*half  !Must be multiplied later by 4*pi/(ucvol*nkbz).
!     end if ! kinetic

      call timab(425,2,tsec)

     end do !jb used to calculate matrix elements of $\Sigma$

     !shaltaf (030406): this has to be done in a clean way later
     if(mod(sp%gwcalctyp,10)==0) then
      if(ppmodel==3.or.ppmodel==4)then
       !MG060926 added spin variable
       sigcme_tmp(:,kb,kb,isppol)=sigcme2(:,kb)
      end if
     end if

    end do !kb to calculate matrix elements of $\Sigma$
   end do !end loop over spin
  end do !ib
 end do !ikbz

!Now we have all the diagonal (off-diagonal) matrix elements for the input k-point

 deallocate(sigcme2,sigcme_3)

! These lines added by Shaltaf for parallelization 10/08/05
! modefied by shaltaf 01/03/06
! Rewritten using upper level primitives by XG 20060628

!MG FIXME FIXME FIXME FIXME
!MG I have modified the following lines inserting a loop over spins
!since there is no subroutine in xsum_mpi.F90 that deals with a complex arrays
!whose shape is dummy(:,:,:,:)
 do is=1,sp%nsppol
  call xsum_mpi(sigcme_tmp(:,:,:,is),spaceComm,ierr)
  call xsum_mpi(sigxme_tmp(:,:,is),spaceComm,ierr)
!END FIXME
 end do


!#          if defined MPI
!  if(mpi_enreg%nproc>1)then       ! collect data from proccessors
!  call MPI_BARRIER(spaceComm,ierr)
!  allocate(xsum1(nomega,sp%nb,sp%nb))
!  call MPI_REDUCE(sigcme_tmp,xsum1,nomega*sp%nb*sp%nb,MPI_COMPLEX,&
!          &  MPI_SUM,master,spaceComm,ierr)
!  sigcme_tmp(:,:,:)=xsum1(:,:,:)
!  deallocate(xsum1)
!  allocate(xsum2(sp%nb,sp%nb))
!  call MPI_REDUCE(sigxme_tmp,xsum2,sp%nb*sp%nb,MPI_COMPLEX,&
!          &  MPI_SUM,master,spaceComm,ierr)
!  sigxme_tmp(:,:)=xsum2(:,:)
!  deallocate(xsum2)
!  end if
!#          end if

! parallelization ends here

 call timab(426,1,tsec)

 !MG060926 added spin
 !multiply by constants
 sigxme_tmp(:,:,:)=(4*pi/(ucvol*sp%nkbz))*sigxme_tmp(:,:,:)
 sigcme_tmp(:,:,:,:)=(4*pi/(ucvol*sp%nkbz))*sigcme_tmp(:,:,:,:)

!MG in case of symmetrization we have to sum over complexes
!FIXME, actually only diagonal terms are considered
 if (lt_k%sym_flag/=0) then
  allocate (dummy_xme(-minbnd+maxbnd+1,-minbnd+maxbnd+1,sp%nsppol))
  allocate (dummy_cme(nomega,-minbnd+maxbnd+1,-minbnd+maxbnd+1,sp%nsppol))
  dummy_cme=zero
  dummy_xme=zero

  !diagonal elements
  do is=1,sp%nsppol
   do ib=1,maxbnd-minbnd+1
    sumcxtab=0
    do jb=1,maxbnd-minbnd+1
     if (cxtab(ib,is,jb)==1) then
      dummy_xme(ib,ib,is)=dummy_xme(ib,ib,is)+sigxme_tmp(minbnd-1+jb,minbnd-1+jb,is)
      !NOTE that frequencies should be equal, this is another reason to use a strict
      !criterion for the tollerance on eigenvalues
      dummy_cme(:,ib,ib,is)=dummy_cme(:,ib,ib,is)+sigcme_tmp(:,minbnd-1+jb,minbnd-1+jb,is)
     end if
     sumcxtab=sumcxtab+cxtab(ib,is,jb)
    end do
    dummy_xme(ib,ib,is)=dummy_xme(ib,ib,is)/sumcxtab
    dummy_cme(:,ib,ib,is) =dummy_cme(:,ib,ib,is)/sumcxtab
    !write(*,*) dummy_xme(ib,ib,is),dummy_cme(1,ib,ib,is)
   end do
  end do

  !off diagonal elements NOT SURE !FIXME
  !sigxme_tmp and sigc_me contain matrix element <n|Sigma |m>
  !in this case we sum only over the right complex for
  !the upper right part of the matrix is modified summing over the complex of |m>,
  !the lower left part is obtained summing over the complex of <n|

   allocate (ndwn(maxbnd-minbnd+1,maxbnd-minbnd+1))
   allocate (nup(maxbnd-minbnd+1,maxbnd-minbnd+1))

   do is=1,sp%nsppol
   ndwn=0
   nup=0

    do ib=1,maxbnd-minbnd+1
     do jb=1,maxbnd-minbnd+1

      if (jb>ib) then
       do kb=1,maxbnd-minbnd+1
        do zb=1,maxbnd-minbnd+1
         !if (zb>kb .and. cxtab(jb,is,zb)==1 .and. cxtab(ib,is,kb)==1) then
         if (zb>kb .and. cxtab(jb,is,zb)==1) then
          dummy_xme(ib,jb,is)=dummy_xme(ib,jb,is)+sigxme_tmp(minbnd-1+kb,minbnd-1+zb,is)
          dummy_cme(:,ib,jb,is)=dummy_cme(:,ib,jb,is)+sigcme_tmp(:,minbnd-1+kb,minbnd-1+zb,is)
          nup(ib,jb)=nup(ib,jb)+1
         end if
        end do
       end do
       if (nup(ib,jb)==0) nup(ib,jb)=1
       dummy_xme(ib,jb,is)=dummy_xme(ib,jb,is)/nup(ib,jb)
       dummy_cme(:,ib,jb,is)=dummy_cme(:,ib,jb,is)/nup(ib,jb)
      end if

      if (jb<ib) then
       do kb=1,maxbnd-minbnd+1
        do zb=1,maxbnd-minbnd+1
         !if (zb<kb .and. cxtab(ib,is,kb)==1 .and. cxtab(jb,is,kb)==1) then
         if (zb<kb .and. cxtab(ib,is,kb)==1) then
          dummy_xme(ib,jb,is)=dummy_xme(ib,jb,is)+sigxme_tmp(minbnd-1+kb,minbnd-1+zb,is)
          dummy_cme(:,ib,jb,is)=dummy_cme(:,ib,jb,is)+sigcme_tmp(:,minbnd-1+kb,minbnd-1+zb,is)
          ndwn(ib,jb)=ndwn(ib,jb)+1
         end if
        end do
       end do
       if (ndwn(ib,jb)==0) ndwn(ib,jb)=1
       dummy_xme(ib,jb,is)=dummy_xme(ib,jb,is)/ndwn(ib,jb)
       dummy_cme(:,ib,jb,is)=dummy_cme(:,ib,jb,is)/ndwn(ib,jb)
      end if

     end do
    end do

   end do

   !copy values
   do is=1,sp%nsppol
    do ib=1,maxbnd-minbnd+1
     do jb=1,maxbnd-minbnd+1
!     if (ib/=jb) cycle
     sigxme_tmp(ib-1+minbnd,jb-1+minbnd,is)=dummy_xme(ib,jb,is)
     sigcme_tmp(:,ib-1+minbnd,jb-1+minbnd,is)=dummy_cme(:,ib,jb,is)
    end do
   end do
  end do

  deallocate (cxtab,dummy_xme,dummy_cme)
  deallocate (ndwn,nup)
 end if
!END MG

 !if non self-consistent, erase all off-diagonal elements
 if(sp%gwcalctyp<20) then
  do jb=minbnd,maxbnd
   do kb=minbnd,maxbnd
    if(jb==kb) cycle
    !MG060923 added new dimension for spin
    sr%hhartree(jb,kb,jkibz,:)=0
   end do
  end do
 end if

 !MG the subroutine has been modified to account for spin
 call diago_hamilt(minbnd,maxbnd,sp%nsppol,sr%hhartree(minbnd:maxbnd,minbnd:maxbnd,jkibz,:),&
&  sigxme_tmp(minbnd:maxbnd,minbnd:maxbnd,:),sigcme_tmp(sr%nomega+ioe0j,minbnd:maxbnd,minbnd:maxbnd,:),&
&  sr%eigvec_qp(minbnd:maxbnd,minbnd:maxbnd,jkibz,:),sr%en_qp_diago(minbnd:maxbnd,jkibz,:))

 !store final results
 !MG added loop over spin
 do is=1,sp%nsppol
  do jb=minbnd,maxbnd

  !bare exchange diagonal element
  sr%sigxme(jb,jkibz,is) = real(sigxme_tmp(jb,jb,is))
  !correlation matrix elements calculated at energy E0
  sr%sigcmee0(jb,jkibz,is)=sigcme_tmp(sr%nomega+ioe0j,jb,jb,is)
  !the values of SigC are linearly interpolated
  !and the fit parameter a is assumed as derivative dSigC/dE
  scme(:)=real(sigcme_tmp(sr%nomega+1:sr%nomega+sp%nomegasrd,jb,jb,is))
  e0pde(:,jb,is)=sr%omegasrd(jb,jkibz,:,is)

  if(sp%nomegasrd==1) then
   smrt=0.
   a=0.
  else
   smrt=linfit(sp%nomegasrd,e0pde(:,jb,is),scme(:),a,b)
  end if

  if( smrt>0.1/Ha_eV) then
   write(message,'(3a)')&
&   ' csigme : WARNING -',ch10,&
&   ' the values of Re Sig_c are not linear '
   call wrtout(6,message,'COLL')
   write(*,*) smrt,a ; write(*,*)'jb = ',jb,scme(:)*Ha_eV
  end if

  sr%dsigmee0(jb,jkibz,is)=cmplx(a,0)
  !calculate Z = ( 1 - dSigma/domega(E0) )^-1
  sr%ze0(jb,jkibz,is) = 1 / ( 1 - sr%dsigmee0(jb,jkibz,is) )
  !calculate DeltaE_GW = E - E0 = (Sigma(E0) - V_xc) / (1 - dSigma/domega)
  sr%degw(jb,jkibz,is) = ( sr%sigxme(jb,jkibz,is)+sr%sigcmee0(jb,jkibz,is)-sr%e0(jb,jkibz,is)+sr%hhartree(jb,jb,jkibz,is) ) / &
&                        ( 1 - sr%dsigmee0(jb,jkibz,is) )

  sr%egw(jb,jkibz,is) = sr%e0(jb,jkibz,is) + sr%degw(jb,jkibz,is)
  !calculate Sigma(E) = Sigma(E0) + (E-E0) dSigma/dE
  sr%sigmee(jb,jkibz,is) = sr%sigxme(jb,jkibz,is)+sr%sigcmee0(jb,jkibz,is) + &
&                        sr%degw(jb,jkibz,is) * sr%dsigmee0(jb,jkibz,is)
  !spectra of Sigma
  do io=1,sr%nomega
   sr%sigcme(jb,jkibz,io,is) = sigcme_tmp(io,jb,jb,is)
   sr%sigxcme(jb,jkibz,io,is) = sr%sigxme(jb,jkibz,is) + sr%sigcme(jb,jkibz,io,is)
  end do

  do io=1,sp%nomegasrd
   sr%sigcmesrd(jb,jkibz,io,is) = sigcme_tmp(sr%nomega+io,jb,jb,is)
   sr%sigxcmesrd(jb,jkibz,io,is) = sr%sigxme(jb,jkibz,is) + sr%sigcmesrd(jb,jkibz,io,is)
  end do

  write (6,'(i5,9f8.3)') jb,sr%e0(jb,jkibz,is)*Ha_eV,&
&  sr%vxcme(jb,jkibz,is)*Ha_eV,sr%sigxme(jb,jkibz,is)*Ha_eV,&
&  real(sr%sigcmee0(jb,jkibz,is))*Ha_eV,
  real(sr%ze0(jb,jkibz,is)),real(sr%dsigmee0(jb,jkibz,is)),&
  real(sr%sigmee(jb,jkibz,is))*Ha_eV,&
  real(sr%degw(jb,jkibz,is))*Ha_eV,real(sr%egw(jb,jkibz,is))*Ha_eV
  write (6,'(i5,9f8.3)') jb,0.0,0.0,0.0,
  aimag(sr%sigcmee0(jb,jkibz,is))*Ha_eV,&
  aimag(sr%ze0(jb,jkibz,is)),aimag(sr%dsigmee0(jb,jkibz,is)),&
  aimag(sr%sigmee(jb,jkibz,is))*Ha_eV,&
  aimag(sr%degw(jb,jkibz,is))*Ha_eV,aimag(sr%egw(jb,jkibz,is))*Ha_eV

  end do !jb
 end do !is

 !These lines added by GMR on 04/05/05 for the calculation
 !of the decomposition of the correlation part of sigma into
 !Coulomb-hole (coh) and screened-exchange (sex)
 if (cohsex) then
!FIXME Check this part
  sigccoh(:,:)=(4*pi/(ucvol*sp%nkbz))*sigccoh(:,:)
  sigcsex(:,:)=(4*pi/(ucvol*sp%nkbz))*sigcsex(:,:)
  write (6,'(/a)') " COH and SEX decomposition of Sig_c(E0) (eV):"
  write (6,'(a5,3a10)') " band","     Sig_c","   Sig_coh","   Sig_sex"

  do is=1,sp%nsppol
   do jb = minbnd,maxbnd

    write (6,'(i5,3f10.5," (Re)")') jb,real(sr%sigcmee0(jb,jkibz,is))*Ha_eV, 
                                     real(sigccoh(jb,is))*Ha_eV, &
                                     real(sigcsex(jb,is))*Ha_eV

    write (6,'(5x,3f10.5," (Im)")') aimag(sr%sigcmee0(jb,jkibz,is))*Ha_eV, &
&                                  aimag(sigccoh(jb,is))*Ha_eV, &
&                                  aimag(sigcsex(jb,is))*Ha_eV
   end do
  end do
 end if

 !These lines added by YMN on 07/04/04 for the calculation
 !of the kinetic and Coulomb contributions to the bandgap energy.
!FIXME This part should be commented because sigckin is not calculated anymore
 if (kinclb) then
  sigckin(:) = (four*pi/(ucvol*sp%nkbz))*sigckin(:)
  write (6,'(/,a)') " Kinetic and Coulomb contributions to Sig_c(E0) (eV):"
  write (6,'(a5,3a10)') " band","     Sig_c","   Sig_c^k","   Sig_c^c"
  do jb = minbnd,maxbnd
!MG FIXME using only spin up in all the terms in the following part
   write (6,'(i5,3f10.5," (Re)")') jb,real(sr%sigcmee0(jb,jkibz,1))*Ha_eV, 
                                     real(sigckin(jb))*Ha_eV, &
                                    (real(sr%sigcmee0(jb,jkibz,1))-real(sigckin(jb)))*Ha_eV

   write (6,'(5x,3f10.5," (Im)")') aimag(sr%sigcmee0(jb,jkibz,1))*Ha_eV, &
&                                  aimag(sigckin(jb))*Ha_eV, &
&                                 (aimag(sr%sigcmee0(jb,jkibz,1))-aimag(sigckin(jb)))*Ha_eV
  end do

  deallocate(sigckin)
 end if

 call timab(426,2,tsec)
 call timab(421,2,tsec)

 deallocate(omegame0i,sigcme_tmp,ket,ket1,ket2,otq,botsq,wfr_jb)
 ! on suggestion by Deyu Lu (PMA)
 if(allocated(sigctmp)) deallocate(sigctmp )
 if(allocated(sigccoh )) deallocate(sigccoh )
 if(allocated(sigcsex )) deallocate(sigcsex )
 if(allocated(eig )) deallocate(eig )
 if(allocated(epsm1q)) deallocate(epsm1q)
 if(allocated(epsm1q_trcc)) deallocate(epsm1q_trcc)
 if(allocated(rhotwg_1 )) deallocate(rhotwg_1 )
!MT 2006-0902  ??
 if(associated(mpi_enreg%proc_distrb)) deallocate(mpi_enreg%proc_distrb)
! if(allocated(mpi_enreg%proc_distrb)) deallocate(mpi_enreg%proc_distrb)

 if(mod(sp%gwcalctyp,10)==7) deallocate (wfg2)

!DEBUG
!write(6,*)' csigme : exit '
!ENDDEBUG

end subroutine csigme
!!***

Generated by  Doxygen 1.6.0   Back to index