The ABINIT Users Mailing List ( CLOSED )

["Ludwig, Christian" <ludwigc@uni-mainz.de>]

Dear Matteo,

 cd 68_gw && make clean && make ends in giving the same error message as 
before.

> Could you provide the CPP preprocessed source for hdr_vs_dtset and 
> m_sigma_results

******************* hdr_vs_dtset_cpp.f90: ****************************

subroutine hdr_vs_dtset(Hdr,Dtset)
 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_errors
 use m_crystal, only : print_symmetries
 use interfaces_14_hidewrite
 use interfaces_68_gw, except_this_one => hdr_vs_dtset
 implicit none
 type(Hdr_type) :: Hdr
 type(Dataset_type),intent(inout) :: Dtset
 logical,parameter :: NO_FATAL=.FALSE.,IS_FATAL=.TRUE.
 integer :: ii,ik,jj,isym,ierr
 logical :: test
 logical :: tsymrel,ttnons,tsymafm
 character(len=500) :: msg
 ! === Check basic dimensions ===
 ierr=0
 call compare_and_change_int(Hdr%natom, Dtset%natom, 'natom', IS_FATAL,ierr)
 call compare_and_change_int(Hdr%nkpt, Dtset%nkpt, 'nkpt', IS_FATAL,ierr)
 call compare_and_change_int(Hdr%npsp, Dtset%npsp, 'npsp', IS_FATAL,ierr)
 call compare_and_change_int(Hdr%nspden, Dtset%nspden, 'nspden', 
IS_FATAL,ierr)
 call 
compare_and_change_int(Hdr%nspinor,Dtset%nspinor,'nspinor',IS_FATAL,ierr)
 call compare_and_change_int(Hdr%nsppol, Dtset%nsppol, 'nsppol', 
IS_FATAL,ierr)
 call compare_and_change_int(Hdr%nsym, Dtset%nsym, 'nsym', IS_FATAL,ierr)
 call compare_and_change_int(Hdr%ntypat, Dtset%ntypat, 'ntypat', 
IS_FATAL,ierr)
 call compare_and_change_int(Hdr%usepaw, Dtset%usepaw, 'usepaw', 
IS_FATAL,ierr)
 call compare_and_change_int(Hdr%usewvl, Dtset%usewvl, 'usewvl', 
IS_FATAL,ierr)
 ! === The number of fatal errors must be zero ===
 if (ierr/=0) then
  write(msg,'(3a)')&
& ' Cannot continue, basic dimensions reported in the header do not agree 
with input file. ',ch10,&
& ' Check consistency between the content of the external file and the input 
file. '
  call msg_hndl(msg,"ERROR" ,"COLL","<stdin>",48)
 end if
 test=ALL(ABS(Hdr%xred-Dtset%xred_orig(:,1:Dtset%natom))<tol6)
 call assert(test,'Mismatch in xred',"<stdin>",52)
 test=ALL(Hdr%typat==Dtset%typat(1:Dtset%natom))
 call assert(test,'Mismatch in typat',"<stdin>",55)
 !
 ! * Check if the lattice from the input file agrees with that read from the 
KSS file
 if ( (ANY(ABS(Hdr%rprimd-Dtset%rprimd_orig)>tol6)) ) then
  write(msg,'(6a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& ' real lattice vectors read from Header ',ch10,&
& ' differ from the values specified in the input file'
  call wrtout(std_out,msg,'COLL')
  write(msg,'(3a,3(3es16.6),3a,3(3es16.6),3a)')ch10,&
& ' rprimd from Hdr file   = ',ch10,(Hdr%rprimd(:,jj),jj=1,3),ch10,&
& ' rprimd from input file = 
',ch10,(Dtset%rprimd_orig(:,jj),jj=1,3),ch10,ch10,&
& '  Modify the lattice vectors in the input file '
  call wrtout(std_out,msg,'COLL') ; call leave_new('COLL')
 end if
 ! === Check symmetry operations ===
 tsymrel=(ALL(Hdr%symrel==Dtset%symrel(:,:,1:Dtset%nsym)))
 if (.not.tsymrel) then
  write(msg,'(6a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& ' real space symmetries read from Header ',ch10,&
& ' differ from the values inferred from the input file'
  call wrtout(std_out,msg,'COLL')
  tsymrel=.FALSE.
 end if
 ttnons=ALL(ABS(Hdr%tnons-Dtset%tnons(:,1:Dtset%nsym))<tol6)
 if (.not.ttnons) then
  write(msg,'(6a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& ' fractional translations read from Header ',ch10,&
& ' differ from the values inferred from the input file'
  call wrtout(std_out,msg,'COLL')
  ttnons=.FALSE.
 end if
 tsymafm=ALL(Hdr%symafm==Dtset%symafm(1:Dtset%nsym))
 if (.not.tsymafm) then
  write(msg,'(6a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& ' AFM symmetries read from Header ',ch10,&
& ' differ from the values inferred from the input file'
  call wrtout(std_out,msg,'COLL')
  tsymafm=.FALSE.
 end if
 if (.not.(tsymrel.and.ttnons.and.tsymafm)) then
  write(msg,'(a)')' Header '
  call wrtout(std_out,msg,'COLL')
  call print_symmetries(Hdr%nsym,Hdr%symrel,Hdr%tnons,Hdr%symafm)
  write(msg,'(a)')' Dtset  '
  call wrtout(std_out,msg,'COLL')
  call print_symmetries(Dtset%nsym,Dtset%symrel,Dtset%tnons,Dtset%symafm)
  write(msg,'(a)')' Check symmetry operations '
  call msg_hndl(msg,"ERROR" ,"COLL","<stdin>",110)
 end if
 ! * Check if the k-points from the input file agrees with that read from the 
KSS file
 if ( (ANY(ABS(Hdr%kptns(:,:)-Dtset%kpt(:,1:Dtset%nkpt))>tol6)) ) then
  write(msg,'(9a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& '  k-points read from Header ',ch10,&
& '  differ from the values specified in the input file',ch10,&
& '  k-points from Hdr file                        | k-points from input file 
',ch10
  call wrtout(std_out,msg,'COLL')
  do ik=1,Dtset%nkpt
   write(msg,'(3(3es16.6,3x))')Hdr%kptns(:,ik),Dtset%kpt(:,ik)
   call wrtout(std_out,msg,'COLL')
  end do
  write(msg,'(a)')' Modify the k-mesh in the input file '
  call msg_hndl(msg,"ERROR" ,"COLL","<stdin>",126)
 end if
 if (ANY(ABS(Hdr%wtk(:)-Dtset%wtk(1:Dtset%nkpt))>tol6)) then
  write(msg,'(9a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& '  k-point weights read from Header ',ch10,&
& '  differ from the values specified in the input file',ch10,&
& '  Hdr file  |  File ',ch10
  call wrtout(std_out,msg,'COLL')
  do ik=1,Dtset%nkpt
   write(msg,'(2(f11.5,1x))')Hdr%wtk(ik),Dtset%wtk(ik)
   call wrtout(std_out,msg,'COLL')
  end do
  write(msg,'(a)')' Check the k-mesh and the symmetries of the system. '
  call msg_hndl(msg,"ERROR" ,"COLL","<stdin>",141)
 end if
 ! Check istwfk storage
 if ( (ANY(Hdr%istwfk(:)/=Dtset%istwfk(1:Dtset%nkpt))) ) then
  write(msg,'(9a)')ch10,&
& ' hdr_vs_dtset : ERROR - ',ch10,&
& '  istwfk read from Header ',ch10,&
& '  differ from the values specified in the input file',ch10,&
& '  Hdr | input ',ch10
  call wrtout(std_out,msg,'COLL')
  do ik=1,Dtset%nkpt
   write(msg,'(i5,3x,i5)')Hdr%istwfk(ik),Dtset%istwfk(ik)
   call wrtout(std_out,msg,'COLL')
  end do
  write(msg,'(a)')' Modify istwfk in the input file '
  call msg_hndl(msg,"ERROR" ,"COLL","<stdin>",157)
 end if
CONTAINS !===========================================================
 subroutine compare_and_change_int(int_exp,int_found,info_string,isfatal,ierr)
 use interfaces_14_hidewrite
 implicit none
 integer,intent(in) :: int_exp
 integer,intent(inout) :: int_found,ierr
 character(len=*),intent(in) :: info_string
 logical,intent(in) :: isfatal
 logical :: leq
 character(len=500) :: msg
 leq=(int_exp==int_found)
 if (.not.leq) then
  write(msg,'(4a,i6,a,i6)')ch10,&
& ' hdr_vs_dtset : WARNING - Mismatch in '//TRIM(info_string),ch10,&
& '  Expected = ',int_exp,' Found = ',int_found
  call wrtout(std_out,msg,'COLL')
  if (isfatal) then
   ! Increase ierr to signal we should stop in the caller.
   ierr=ierr+1
  else
   ! Change the value found to reflect the expected one.
   int_found=int_exp
   write(msg,'(a)')' The value found has been set to the expected one.'
   call wrtout(std_out,msg,'COLL')
  end if
 end if
 end subroutine compare_and_change_int
end subroutine hdr_vs_dtset


******************* m_sigma_results_cpp.f90: **************************

MODULE m_sigma_results
 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_errors
 implicit none
 private
 public :: &
& write_sigma_results_header, &
& write_sigma_results, &
& print_Sigma_perturbative, &
& print_Sigma_QPSC, &
& nullify_sigma_results, &
& init_sigma_results, &
& destroy_sigma_results, &
& allocate_sigma_results, &
& etsf_dump_QP, &
& abi_etsf_get_QP
CONTAINS 
!===================================================================================================================
======
subroutine write_sigma_results_header(Sp,Er,Cryst,Kmesh,Qmesh)
 use defs_basis
 use interfaces_14_hidewrite
 implicit none
 type(Bz_mesh_type),intent(in) :: Kmesh,Qmesh
 type(Crystal_structure),intent(in) :: Cryst
 type(Epsilonm1_results),intent(in) :: Er
 type(Sigma_parameters),intent(in) :: Sp
 integer :: mod10
 character(len=500) :: msg

 write(msg,'(a)')' SIGMA fundamental parameters:'
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 mod10=MOD(Sp%gwcalctyp,10)
 SELECT CASE (mod10)
 CASE (0)
  write(msg,'(a,i2)')' PLASMON POLE MODEL ',Sp%ppmodel
 CASE (1)
  write(msg,'(a)')' ANALYTIC CONTINUATION'
 CASE (2)
  write(msg,'(a)')' CONTOUR DEFORMATION'
CASE (5)
  write(msg,'(a)')' Hartree-Fock'
 CASE (6)
  write(msg,'(a)')' Screened Exchange'
 CASE (7)
  write(msg,'(a)')' COHSEX'
 CASE (8)
  write(msg,'(a,i2)')' MODEL GW with PLASMON POLE MODEL ',Sp%ppmodel
 CASE (9)
  write(msg,'(a)')' MODEL GW without PLASMON POLE MODEL'
 CASE DEFAULT
  write(msg,'(a,i3)')' Wrong value for Sp%gwcalctyp = ',Sp%gwcalctyp
  call msg_hndl(msg,"BUG" ,"COLL","<stdin>",76)
 END SELECT
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of plane-waves for SigmaX         ',Sp%npwx
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of plane-waves for SigmaC and W   ',Sp%npwc
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of plane-waves for wavefunctions  ',Sp%npwwfn
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of bands                          ',Sp%nbnds
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of independent spin polarizations ',Sp%nsppol
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of spinorial components           ',Sp%nspinor
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of k-points in IBZ                ',Kmesh%nibz
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of q-points in IBZ                ',Qmesh%nibz
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of symmetry operations            ',Cryst%nsym
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of k-points in BZ                 ',Kmesh%nbz
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of q-points in BZ                 ',Qmesh%nbz
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of frequencies for dSigma/dE      ',Sp%nomegasrd
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,f12.2)')' frequency step for dSigma/dE [eV]        
',Sp%deltae*Ha_eV
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
write(msg,'(a,i12)')' number of omega for Sigma on real axis   ',Sp%nomegasr
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,f12.2)')' max omega for Sigma on real axis  [eV]   
',Sp%maxomega_r*Ha_eV
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,f12.2)')' zcut for avoiding poles [eV]             
',Sp%zcut*Ha_eV
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 if (Sp%soenergy>0.1d-4) then
  write(msg,'(a,f12.2)')' scissor energy [eV]                      
',Sp%soenergy*Ha_eV
  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 end if
 if (MOD(Sp%gwcalctyp,10)==1) then
  write(msg,'(a,i12)')' number of imaginary frequencies for Sigma',Sp%nomegasi
  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
  write(msg,'(a,f12.2)')' max omega for Sigma on imag axis  [eV]   
',Sp%omegasimax*Ha_eV
  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 end if
 write(msg,'(2a)')ch10,' EPSILON^-1 parameters (SCR file):'
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 !write(std_out,*) titem1(2)(1:79)
 write(msg,'(a,i12)')' dimension of the eps^-1 matrix on file   ',Er%Hscr%npwe
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' dimension of the eps^-1 matrix used      ',Er%npwe
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of plane-waves for wavefunctions  
',Er%Hscr%npwwfn_used
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of bands                          
',Er%Hscr%nbnds_used
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of q-points in IBZ                ',Qmesh%nibz
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of frequencies                    ',Er%nomega
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of real frequencies               ',Er%nomega_r
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(a,i12)')' number of imag frequencies               ',Er%nomega_i
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 write(msg,'(3a)')ch10,' matrix elements of self-energy operator (all in 
[eV])',ch10
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
 if (Sp%gwcalctyp<10) then
  write(msg,'(a)')' Perturbative Calculation'
 else if (Sp%gwcalctyp<20) then
  write(msg,'(a)')' Self-Consistent on Energies only'
 else
  write(msg,'(a)')' Self-Consistent on Energies and Wavefunctions'
 end if
 call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')

end subroutine write_sigma_results_header
subroutine write_sigma_results(ikcalc,ikibz,Sp,Sr,Kmesh,KS_BSt)
 use defs_basis
 use m_gwdefs, only : unt_gw, unt_sig, unt_sgr, unt_sgm
 use interfaces_14_hidewrite
 implicit none
 integer,intent(in) :: ikcalc,ikibz
 type(Bandstructure_type),intent(in) :: KS_BSt
 type(BZ_mesh_type),intent(in) :: Kmesh
 type(Sigma_parameters),intent(in) :: Sp
 type(Sigma_results),intent(in) :: Sr
 integer :: ib,io,is
 character(len=500) :: msg
 character(len=12) :: tag_spin(2)
 real(dp),pointer :: ks_energy(:,:,:)
 !unt_gw = 21 ! File with GW corrections.
 !unt_sig = 22 ! Self-energy as a function of frequency.
 !unt_sgr = 23 ! Derivative wrt omega of the Self-energy.
 !unt_sgm = 20 ! Sigma on the Matsubara axis.
 !TODO pass the object instead of the pointer
 ks_energy => KS_BSt%eig(:,:,:)
 tag_spin=(/'            ','            '/)
 if (Sr%nsppol==2) tag_spin=(/',  SPIN UP  ',',  SPIN DOWN'/)
 do is=1,Sr%nsppol
  write(msg,'(2a,3f8.3,a)')ch10,' k = ',Sp%xkcalc(:,ikcalc),tag_spin(is)
  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
  write(msg,'(a)')&
& ' Band     E0 <VxcLDA>   SigX SigC(E0)      Z dSigC/dE  Sig(E)    E-E0      
 E'
  if (Sr%usepawu/=0) then
   write(msg,'(a)')&
& ' Band     E0 <VxcLDA>   <H_U>  SigX SigC(E0)      Z dSigC/dE  Sig(E)    
E-E0       E'
  end if
  if (Sp%gwcalctyp>=10) then
   write(msg,'(2a)')&
& ' Band     E_lda   <Vxclda>   E(N-1)  <Hhartree>   SigX  SigC[E(N-1)]',&
& '    Z     dSigC/dE  Sig[E(N)]  DeltaE  E(N)_pert E(N)_diago'
  end if
  call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
  write(unt_gw,'(3f10.6)')Sp%xkcalc(:,ikcalc)
  write(unt_gw,'(i4)')Sp%maxbnd(ikcalc)-Sp%minbnd(ikcalc)+1
  write(unt_sig,'("# k = ",3f10.6)')Sp%xkcalc(:,ikcalc)
  write(unt_sig,'("# b = ",2i10)')Sp%minbnd(ikcalc),Sp%maxbnd(ikcalc)
  write(unt_sgr,'("# k = ",3f10.6)')Sp%xkcalc(:,ikcalc)
  write(unt_sgr,'("# b = ",2i10)')Sp%minbnd(ikcalc),Sp%maxbnd(ikcalc)
  do ib=Sp%minbnd(ikcalc),Sp%maxbnd(ikcalc)
   if (Sp%gwcalctyp>=10) then
    call print_Sigma_QPSC(Sr,ikibz,ib,is,Kmesh,ks_energy,unit=ab_out)
    call 
print_Sigma_QPSC(Sr,ikibz,ib,is,Kmesh,ks_energy,unit=std_out,prtvol=1)
   else
    !
    ! === When using non-ppm, write out also the imaginary part in ab_out
    SELECT CASE(Sp%gwcalctyp)
    CASE(1,2)
     call print_Sigma_perturbative(Sr,ikibz,ib,is,unit=ab_out,prtvol=1)
    CASE DEFAULT
     call print_Sigma_perturbative(Sr,ikibz,ib,is,unit=ab_out)
    END SELECT
    call print_Sigma_perturbative(Sr,ikibz,ib,is,unit=std_out,prtvol=1)
   end if
   write(unt_gw,'(i6,3f9.4)') &
& ib, &
& REAL (Sr%egw(ib,ikibz,is)) *Ha_eV,&
& REAL (Sr%degw(ib,ikibz,is))*Ha_eV,&
& AIMAG(Sr%egw(ib,ikibz,is)) *Ha_eV
  end do !ib
  if (Sr%e0gap(ikibz,is)**2+Sr%egwgap(ikibz,is)**2+Sr%degwgap(ikibz,is)**2 > 
tol10) then
   ! Output the direct gap for each spin
   ! If all the gaps are zero, this means that it could not be computed in 
the calling routine
   write(msg,'(2a,f8.3)')ch10,' E^0_gap       ',Sr%e0gap(ikibz,is)*Ha_eV
   call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
   write(msg,'(a,f8.3)') ' E^GW_gap      ',Sr%egwgap(ikibz,is)*Ha_eV
   call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
   write(msg,'(a,f8.3,a)') ' DeltaE^GW_gap ',Sr%degwgap(ikibz,is)*Ha_eV,ch10
   call wrtout(std_out,msg,'COLL') ; call wrtout(ab_out,msg,'COLL')
  end if
  ! === Output spectral function ===
  do io=1,Sr%nomega_r
   write(unt_sig,'(100(e11.5,2x))')&
& REAL(Sr%omega_r(io))*Ha_eV,&
& (REAL(Sr%sigxcme(ib,ikibz,io,is))*Ha_eV,&
& AIMAG(Sr%sigxcme(ib,ikibz,io,is))*Ha_eV,&
& one/pi*ABS(AIMAG(Sr%sigcme(ib,ikibz,io,is)))&
& /( 
(REAL(Sr%omega_r(io)-Sr%hhartree(ib,ib,ikibz,is)-Sr%sigxcme(ib,ikibz,io,is)))**2&
& +(AIMAG(Sr%sigcme(ib,ikibz,io,is)))**2) /Ha_eV,&
& ib=Sp%minbnd(ikcalc),Sp%maxbnd(ikcalc))
  end do
  do ib=Sp%minbnd(ikcalc),Sp%maxbnd(ikcalc)
   write(unt_sgr,'("# ik, ib",2i5)')ikibz,ib
   do io=1,Sr%nomega4sd
    write(unt_sgr,'(100(e11.5,2x))') &
& REAL (Sr%omega4sd (ib,ikibz,io,is)) *Ha_eV,&
& REAL (Sr%sigxcme4sd(ib,ikibz,io,is)) *Ha_eV,&
& AIMAG(Sr%sigxcme4sd(ib,ikibz,io,is)) *Ha_eV
   end do
  end do
  if (MOD(sp%gwcalctyp,10)==1) then
   ! For AC write matrix elements of sigma along the imaginary axis
   do ib=Sp%minbnd(ikcalc),Sp%maxbnd(ikcalc)
    write(unt_sgm,'("# ik, ib",2i5)')ikibz,ib
    do io=1,Sr%nomega_i
     write(unt_sgm,'(3(e11.5,2x))') &
& AIMAG(Sr%omega_i(io)) *Ha_eV,&
& REAL (Sr%sigxcmesi(ib,ikibz,io,is))*Ha_eV,&
& AIMAG(Sr%sigxcmesi(ib,ikibz,io,is))*Ha_eV
    end do
   end do
  end if
 end do !is
end subroutine write_sigma_results
subroutine 
print_Sigma_perturbative(Sr,ik_ibz,iband,isp,unit,prtvol,mode_paral,witheader)
 use defs_basis
 use interfaces_14_hidewrite
 implicit none
 integer,intent(in) :: iband,ik_ibz,isp
 integer,optional,intent(in) :: prtvol,unit
 character(len=4),optional,intent(in) :: mode_paral
 logical,optional,intent(in) :: witheader
 type(Sigma_results),intent(in) :: Sr
 integer :: unt,verbose
 logical :: ltest
 character(len=4) :: mode
 character(len=500) :: msg
 unt=std_out ; if (PRESENT(unit)) unt=unit
 verbose=0 ; if (PRESENT(prtvol)) verbose=prtvol
 mode='COLL' ; if (PRESENT(mode_paral)) mode=mode_paral
 if (PRESENT(witheader)) then
  if (witheader) then
   write(msg,'(a)')' Band     E0 <VxcLDA>   SigX SigC(E0)      Z dSigC/dE  
Sig(E)    E-E0       E '
   call wrtout(unt,msg,mode)
  end if
 end if
 if (Sr%usepawu==0) then
  if (Sr%nsig_ab/=1) then
   write(msg,'(i5,9f8.3)') &
& iband, &
& Sr%e0 (iband,ik_ibz,1)*Ha_eV, &
& SUM(Sr%vxcme (iband,ik_ibz,:))*Ha_eV, &
& SUM(Sr%sigxme (iband,ik_ibz,:))*Ha_eV, &
& REAL(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV,&
& REAL(Sr%ze0 (iband,ik_ibz,1)), &
& REAL(SUM(Sr%dsigmee0(iband,ik_ibz,:))), &
& REAL(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV,&
& REAL(Sr%degw (iband,ik_ibz,1))*Ha_eV, &
& REAL(Sr%egw (iband,ik_ibz,1))*Ha_eV
    call wrtout(unt,msg,mode)
   if (verbose/=0) then
    write(msg,'(i5,9f8.3)') &
& iband, &
& zero, &
& zero, &
& zero, &
& AIMAG(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV,&
& AIMAG(Sr%ze0 (iband,ik_ibz,1)), &
& AIMAG(SUM(Sr%dsigmee0(iband,ik_ibz,:))), &
& AIMAG(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV,&
& AIMAG(Sr%degw (iband,ik_ibz,1))*Ha_eV, &
& AIMAG(Sr%egw (iband,ik_ibz,1))*Ha_eV
     call wrtout(unt,msg,mode)
   end if
  else
   write(msg,'(i5,9f8.3)') &
& iband, &
& Sr%e0 (iband,ik_ibz,isp)*Ha_eV, &
& Sr%vxcme (iband,ik_ibz,isp)*Ha_eV, &
& Sr%sigxme (iband,ik_ibz,isp)*Ha_eV, &
& REAL(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,&
& REAL(Sr%ze0 (iband,ik_ibz,isp)), &
& REAL(Sr%dsigmee0(iband,ik_ibz,isp)), &
& REAL(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,&
& REAL(Sr%degw (iband,ik_ibz,isp))*Ha_eV,&
& REAL(Sr%egw (iband,ik_ibz,isp))*Ha_eV
    call wrtout(unt,msg,mode)
   if (verbose/=0) then
    write(msg,'(i5,9f8.3)') &
& iband, &
& zero, &
& zero, &
& zero, &
& AIMAG(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,&
& AIMAG(Sr%ze0 (iband,ik_ibz,isp)), &
& AIMAG(Sr%dsigmee0(iband,ik_ibz,isp)), &
& AIMAG(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,&
& AIMAG(Sr%degw (iband,ik_ibz,isp))*Ha_eV,&
& AIMAG(Sr%egw (iband,ik_ibz,isp))*Ha_eV
     call wrtout(unt,msg,mode)
   end if
  end if
 else
  ! === PAW+U+GW calculation ===
  ltest=(Sr%nsig_ab==1)
  call assert(ltest,'LDA+U with spinor not implemented',&
& "<stdin>",392)
  write(msg,'(i5,10f8.3)') &
& iband, &
& Sr%e0 (iband,ik_ibz,isp)*Ha_eV, &
& Sr%vxcme (iband,ik_ibz,isp)*Ha_eV, &
& Sr%vUme (iband,ik_ibz,isp)*Ha_eV, &
& Sr%sigxme (iband,ik_ibz,isp)*Ha_eV, &
& REAL(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,&
& REAL(Sr%ze0 (iband,ik_ibz,isp)), &
& REAL(Sr%dsigmee0(iband,ik_ibz,isp)), &
& REAL(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,&
& REAL(Sr%degw (iband,ik_ibz,isp))*Ha_eV,&
& REAL(Sr%egw (iband,ik_ibz,isp))*Ha_eV
   call wrtout(unt,msg,mode)
  if (verbose/=0) then
   write(msg,'(i5,10f8.3)') &
& iband, &
& zero, &
& zero, &
& zero, &
& zero, &
& AIMAG(Sr%sigcmee0(iband,ik_ibz,isp))*Ha_eV,&
& AIMAG(Sr%ze0 (iband,ik_ibz,isp)), &
& AIMAG(Sr%dsigmee0(iband,ik_ibz,isp)), &
& AIMAG(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV,&
& AIMAG(Sr%degw (iband,ik_ibz,isp))*Ha_eV,&
& AIMAG(Sr%egw (iband,ik_ibz,isp))*Ha_eV
    call wrtout(unt,msg,mode)
  end if
 end if
end subroutine print_Sigma_perturbative
subroutine 
print_Sigma_QPSC(Sr,ik_ibz,iband,isp,Kmesh,ks_energy,unit,prtvol,mode_paral)
 use defs_basis
 use interfaces_14_hidewrite
 implicit none
 integer,intent(in) :: iband,ik_ibz,isp
 integer,intent(in),optional :: prtvol,unit
 character(len=4),intent(in),optional :: mode_paral
 type(BZ_mesh_type),intent(in) :: Kmesh
 type(Sigma_results),intent(in) :: Sr
 real(dp),intent(in) :: ks_energy(Sr%nbnds,Kmesh%nibz,Sr%nsppol)
 integer :: unt,verbose
 character(len=4) :: mode
 character(len=500) :: msg
 unt =std_out ; if (PRESENT(unit)) unt =unit
 verbose=0 ; if (PRESENT(prtvol)) verbose=prtvol
 mode ='COLL' ; if (PRESENT(mode_paral)) mode =mode_paral
 if (Sr%usepawu==0) then
  if (Sr%nsig_ab/=1) then
   write(msg,'(i5,12(2x,f8.3))') &
& iband, &
& ks_energy (iband,ik_ibz,1)*Ha_eV, &
         !sr%e0 (iband,ikibz,1)*Ha_eV, &
& SUM(Sr%vxcme (iband,ik_ibz,:))*Ha_eV, &
& Sr%e0 (iband,ik_ibz,1)*Ha_eV, &
& REAL(SUM(Sr%hhartree(iband,iband,ik_ibz,:)))*Ha_eV,&
& SUM(Sr%sigxme (iband,ik_ibz,:))*Ha_eV, &
& REAL(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV, &
& REAL(Sr%ze0 (iband,ik_ibz,1)), &
& REAL(SUM(Sr%dsigmee0(iband,ik_ibz,:))), &
& REAL(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV, &
& REAL(Sr%degw (iband,ik_ibz,1))*Ha_eV, &
& REAL(Sr%egw (iband,ik_ibz,1))*Ha_eV, &
& Sr%en_qp_diago (iband,ik_ibz,1)*Ha_eV
   call wrtout(unt,msg,mode)
   write(msg,'(i5,12(2x,f8.3))') &
& iband, &
& zero, &
          !sr%e0 (iband,ikibz,isp)*Ha_eV, &
& zero, &
& zero, &
& AIMAG(SUM(Sr%hhartree(iband,iband,ik_ibz,:)))*Ha_eV,&
& zero, &
& AIMAG(SUM(Sr%sigcmee0(iband,ik_ibz,:)))*Ha_eV, &
& AIMAG(Sr%ze0 (iband,ik_ibz,1)), &
& AIMAG(SUM(Sr%dsigmee0(iband,ik_ibz,:))), &
& AIMAG(SUM(Sr%sigmee (iband,ik_ibz,:)))*Ha_eV, &
& AIMAG(Sr%degw (iband,ik_ibz,1))*Ha_eV, &
& AIMAG(Sr%egw (iband,ik_ibz,1))*Ha_eV, &
& zero
   if (verbose/=0) call wrtout(unt,msg,mode)
  else
   write(msg,'(i5,12(2x,f8.3))') &
& iband, &
& ks_energy (iband,ik_ibz,isp)*Ha_eV, &
         !sr%e0 (iband,ikibz,isp)*Ha_eV, &
& Sr%vxcme (iband,ik_ibz,isp)*Ha_eV, &
& Sr%e0 (iband,ik_ibz,isp)*Ha_eV, &
& REAL(Sr%hhartree (iband,iband,ik_ibz,isp))*Ha_eV,&
& Sr%sigxme (iband,ik_ibz,isp)*Ha_eV, &
& REAL(Sr%sigcmee0 (iband,ik_ibz,isp))*Ha_eV, &
& REAL(Sr%ze0 (iband,ik_ibz,isp)), &
& REAL(Sr%dsigmee0 (iband,ik_ibz,isp)), &
& REAL(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV, &
& REAL(Sr%degw (iband,ik_ibz,isp))*Ha_eV, &
& REAL(Sr%egw (iband,ik_ibz,isp))*Ha_eV, &
& Sr%en_qp_diago(iband,ik_ibz,isp)*Ha_eV
   call wrtout(unt,msg,mode)
   write(msg,'(i5,12(2x,f8.3))') &
& iband, &
& zero, &
          !sr%e0 (iband,ikibz,isp)*Ha_eV, &
& zero, &
& zero, &
& AIMAG(Sr%hhartree (iband,iband,ik_ibz,isp))*Ha_eV,&
& zero, &
& AIMAG(Sr%sigcmee0 (iband,ik_ibz,isp))*Ha_eV, &
& AIMAG(Sr%ze0 (iband,ik_ibz,isp)), &
& AIMAG(Sr%dsigmee0 (iband,ik_ibz,isp)), &
& AIMAG(Sr%sigmee (iband,ik_ibz,isp))*Ha_eV, &
& AIMAG(Sr%degw (iband,ik_ibz,isp))*Ha_eV, &
& AIMAG(Sr%egw (iband,ik_ibz,isp))*Ha_eV, &
& zero
   if (verbose/=0) call wrtout(unt,msg,mode)
  end if
 else
  ! === PAW+U+GW calculation ===
  msg="PAW+U+GW not yet implemented"
  call msg_hndl(msg,"ERROR" ,"COLL","<stdin>",525)
 end if
end subroutine print_Sigma_QPSC
subroutine nullify_sigma_results(Sr)
 use defs_basis
 implicit none
 type(Sigma_results),intent(inout) :: Sr
 nullify(Sr%maxbnd )
 nullify(Sr%minbnd )
 !nullify(Sr%ame )
 nullify(Sr%degwgap )
 nullify(Sr%egwgap )
 nullify(Sr%en_qp_diago)
 nullify(Sr%e0 )
 nullify(Sr%e0gap )
 nullify(Sr%omega_r )
 nullify(Sr%xkcalc )
 nullify(Sr%sigxme )
 nullify(Sr%vxcme )
 nullify(Sr%vUme )
 nullify(Sr%degw )
 nullify(Sr%dsigmee0 )
 nullify(Sr%egw )
 nullify(Sr%eigvec_qp )
 nullify(Sr%hhartree )
 nullify(Sr%sigcme )
 nullify(Sr%sigmee )
 nullify(Sr%sigcmee0 )
 nullify(Sr%sigcmesi )
 nullify(Sr%sigcme4sd )
 nullify(Sr%sigxcme )
 nullify(Sr%sigxcmesi )
 nullify(Sr%sigxcme4sd )
 nullify(Sr%ze0 )
 nullify(Sr%omega_i )
 nullify(Sr%omega4sd )
end subroutine nullify_sigma_results
subroutine init_sigma_results(Sp,nkibz,usepawu,Sr)
 use defs_basis
 implicit none
 integer,intent(in) :: nkibz,usepawu
 type(Sigma_parameters),intent(in) :: Sp
 type(Sigma_results),intent(inout) :: Sr
 !character(len=500) :: msg
 integer :: b1gw,b2gw,mod10
 call nullify_sigma_results(Sr)
 ! === Copy important dimensions ===
 mod10=MOD(Sp%gwcalctyp,10)
 Sr%nkcalc =Sp%nkcalc
 Sr%gwcalctyp =Sp%gwcalctyp
 Sr%deltae =Sp%deltae
 Sr%maxomega4sd=Sp%maxomega4sd
 Sr%maxomega_r =Sp%maxomega_r
 Sr%scissor_ene=Sp%soenergy
 !FIXME this should be done in allocate_sigma_results
 allocate(Sr%minbnd(Sr%nkcalc),Sr%maxbnd(Sr%nkcalc))
 Sr%minbnd=Sp%minbnd ; Sr%maxbnd=Sp%maxbnd
 allocate(Sr%xkcalc(3,Sr%nkcalc))
 Sr%xkcalc=Sp%xkcalc
 Sr%b1gw =Sp%minbdgw ! * min and Max GW band index over k and spin.
 Sr%b2gw =Sp%maxbdgw ! Used to dimension arrays.
 Sr%nbnds =Sp%nbnds
 Sr%nkibz =nkibz
 Sr%nsppol =Sp%nsppol
 Sr%nsig_ab =Sp%nsig_ab
 Sr%nomega_r =Sp%nomegasr !FIXME change name
 Sr%nomega_i =Sp%nomegasi
 Sr%nomega4sd=Sp%nomegasrd
 Sr%usepawu =usepawu
 !======================================================
 ! === Allocate arrays in the sigma_results datatype ===
 !======================================================
 b1gw=Sr%b1gw
 b2gw=Sr%b2gw
 !TODO use this routine
 ! hhartree(b1,b2,k,s)= <b1,k,s|T+v_{loc}+v_{nl}+v_{H}|b2,k,s>
 allocate(Sr%hhartree(b1gw:b2gw,b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 Sr%hhartree=czero
 ! === QP amplitudes and energies ===
 allocate(Sr%en_qp_diago(Sr%nbnds,Sr%nkibz,Sr%nsppol)) ; 
Sr%en_qp_diago(:,:,:)=zero
 allocate(Sr%eigvec_qp(Sr%nbnds,Sr%nbnds,Sr%nkibz,Sr%nsppol)) ; 
Sr%eigvec_qp(:,:,:,:)=czero
 ! Dont know if it is better to do this here or in the sigma
 ! * Initialize with KS wavefunctions and energies
 !do ib=1,Sr%nbnds
 ! Sr%en_qp_diago(ib,:,:)=en(:,ib,:)
 ! Sr%eigvec_qp(ib,ib,:,:)=cone
 !end do
 allocate(Sr%vxcme (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%vUme (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%sigxme (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%sigcme (b1gw:b2gw,Sr%nkibz,Sr%nomega_r,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%sigxcme (b1gw:b2gw,Sr%nkibz,Sr%nomega_r,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%sigcmee0(b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%ze0 (b1gw:b2gw,Sr%nkibz,Sr%nsppol))
 allocate(Sr%dsigmee0(b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%sigmee (b1gw:b2gw,Sr%nkibz,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%degw (b1gw:b2gw,Sr%nkibz,Sr%nsppol))
 allocate(Sr%e0 (Sr%nbnds,Sr%nkibz,Sr%nsppol))
 allocate(Sr%egw(Sr%nbnds,Sr%nkibz,Sr%nsppol))
 allocate(Sr%e0gap (Sr%nkibz,Sr%nsppol))
 allocate(Sr%degwgap(Sr%nkibz,Sr%nsppol))
 allocate(Sr%egwgap (Sr%nkibz,Sr%nsppol))
 !allocate(Sr%ame(Sr%nbnds,Sr%nkibz,Sr%nomega_r))
 !
 ! === These quantities are used to evaluate $\Sigma(E)$ around the KS\QP 
eigenvalue ===
 allocate(Sr%omega4sd (b1gw:b2gw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol))
 allocate(Sr%sigcme4sd (b1gw:b2gw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol*Sr%nsig_ab))
 allocate(Sr%sigxcme4sd(b1gw:b2gw,Sr%nkibz,Sr%nomega4sd,Sr%nsppol*Sr%nsig_ab))
 !TODO Find better treatment
 ! Mesh along the real axis.
 if (Sr%nomega_r>0) then
  allocate(Sr%omega_r(Sr%nomega_r))
  Sr%omega_r(:)=Sp%omega_r(:)
 end if
 Sr%e0 =zero
 Sr%egw =czero
 Sr%e0gap =zero
 Sr%sigcme =czero
 Sr%sigxme =czero
 Sr%sigxcme =czero
 Sr%sigcmee0 =czero
 Sr%ze0 =czero
 Sr%dsigmee0 =czero
 Sr%sigmee =czero
 Sr%omega4sd =czero
 Sr%sigcme4sd =czero
 Sr%sigxcme4sd=czero
 Sr%degw =czero
 ! === Analytical Continuation ===
 if (mod10==1) then
  ! FIXME omegasi should not be in Sp% here we should construct the mesh
  allocate(Sr%omega_i(Sr%nomega_i)) ; Sr%omega_i=Sp%omegasi
  allocate(Sr%sigcmesi (b1gw:b2gw,Sr%nkibz,Sr%nomega_i,Sr%nsppol*Sr%nsig_ab))
  allocate(Sr%sigxcmesi(b1gw:b2gw,Sr%nkibz,Sr%nomega_i,Sr%nsppol*Sr%nsig_ab))
  Sr%sigcmesi =czero
  Sr%sigxcmesi=czero
 end if
end subroutine init_sigma_results
subroutine destroy_sigma_results(Sr)
 use defs_basis
 implicit none
 type(Sigma_results),intent(inout) :: Sr

 if (associated(Sr%maxbnd )) deallocate(Sr%maxbnd)
 if (associated(Sr%minbnd )) deallocate(Sr%minbnd)
 !if (associated(Sr%ame )) deallocate(Sr%ame)
 if (associated(Sr%degwgap )) deallocate(Sr%degwgap)
 if (associated(Sr%egwgap )) deallocate(Sr%egwgap)
 if (associated(Sr%en_qp_diago)) deallocate(Sr%en_qp_diago)
 if (associated(Sr%e0 )) deallocate(Sr%e0)
 if (associated(Sr%e0gap )) deallocate(Sr%e0gap)
 if (associated(Sr%omega_r )) deallocate(Sr%omega_r)
 if (associated(Sr%xkcalc )) deallocate(Sr%xkcalc)
 if (associated(Sr%sigxme )) deallocate(Sr%sigxme)
 if (associated(Sr%vxcme )) deallocate(Sr%vxcme)
 if (associated(Sr%vUme )) deallocate(Sr%vUme)
 if (associated(Sr%degw )) deallocate(Sr%degw)
 if (associated(Sr%dsigmee0 )) deallocate(Sr%dsigmee0)
 if (associated(Sr%egw )) deallocate(Sr%egw)
 if (associated(Sr%eigvec_qp )) deallocate(Sr%eigvec_qp)
 if (associated(Sr%hhartree )) deallocate(Sr%hhartree)
 if (associated(Sr%sigcme )) deallocate(Sr%sigcme)
 if (associated(Sr%sigmee )) deallocate(Sr%sigmee)
 if (associated(Sr%sigcmee0 )) deallocate(Sr%sigcmee0)
 if (associated(Sr%sigcmesi )) deallocate(Sr%sigcmesi)
 if (associated(Sr%sigcme4sd )) deallocate(Sr%sigcme4sd)
 if (associated(Sr%sigxcme )) deallocate(Sr%sigxcme)
 if (associated(Sr%sigxcmesi )) deallocate(Sr%sigxcmesi)
 if (associated(Sr%sigxcme4sd )) deallocate(Sr%sigxcme4sd)
 if (associated(Sr%ze0 )) deallocate(Sr%ze0)
 if (associated(Sr%omega_i )) deallocate(Sr%omega_i)
 if (associated(Sr%omega4sd )) deallocate(Sr%omega4sd)

end subroutine destroy_sigma_results
subroutine 
allocate_sigma_results(Sr,b1gw,b2gw,nbnds,nkibz,nkcalc,nsppol,nsig_ab,nomega_r,nomega_i,nomega4sd,&
& omega_r,omega_i) ! Optional
 use defs_basis
 implicit none
 integer,intent(in) :: b1gw,b2gw,nkibz,nsppol,nsig_ab,nbnds
 integer,intent(in) :: nomega_r,nomega_i,nomega4sd,nkcalc
 type(Sigma_results),intent(inout) :: Sr
 complex(dpc),optional,intent(in) :: omega_r(:),omega_i(:)
 integer :: ii
 call nullify_sigma_results(Sr)
 Sr%nkcalc=nkcalc
 allocate(Sr%minbnd(Sr%nkcalc),Sr%maxbnd(Sr%nkcalc))
 allocate(Sr%xkcalc(3,Sr%nkcalc))
 ! hhartree(b1,b2,k,s)= <b1,k,s|T+v_{loc}+v_{nl}+v_{H}|b2,k,s>
 allocate(Sr%hhartree(b1gw:b2gw,b1gw:b2gw,nkibz,nsppol*nsig_ab))
 Sr%hhartree=czero
 ! === QP amplitudes and energies ===
 allocate(Sr%en_qp_diago(nbnds,nkibz,nsppol))
 allocate(Sr%eigvec_qp(nbnds,nbnds,nkibz,nsppol))
 Sr%en_qp_diago=zero
 Sr%eigvec_qp =czero
 allocate(Sr%vxcme (b1gw:b2gw,nkibz,nsppol*nsig_ab))
 allocate(Sr%vUme (b1gw:b2gw,nkibz,nsppol*nsig_ab))
 allocate(Sr%sigxme (b1gw:b2gw,nkibz,nsppol*nsig_ab))
 allocate(Sr%sigcme (b1gw:b2gw,nkibz,nomega_r,nsppol*nsig_ab))
 allocate(Sr%sigxcme (b1gw:b2gw,nkibz,nomega_r,nsppol*nsig_ab))
 allocate(Sr%sigcmee0(b1gw:b2gw,nkibz,nsppol*nsig_ab))
 allocate(Sr%ze0 (b1gw:b2gw,nkibz,nsppol))
 allocate(Sr%dsigmee0(b1gw:b2gw,nkibz,nsppol*nsig_ab))
 allocate(Sr%sigmee (b1gw:b2gw,nkibz,nsppol*nsig_ab))
 allocate(Sr%degw (b1gw:b2gw,nkibz,nsppol))
 allocate(Sr%e0 (nbnds,nkibz,nsppol))
 allocate(Sr%egw(nbnds,nkibz,nsppol))
 allocate(Sr%e0gap (nkibz,nsppol))
 allocate(Sr%degwgap(nkibz,nsppol))
 allocate(Sr%egwgap (nkibz,nsppol))
 !allocate(Sr%ame(nbnds,nkibz,nomega_r))
 ! === These quantities are used to evaluate $\Sigma(E)$ around the KS\QP 
eigenvalue ===
 allocate(Sr%omega4sd (b1gw:b2gw,nkibz,nomega4sd,nsppol))
 allocate(Sr%sigcme4sd (b1gw:b2gw,nkibz,nomega4sd,nsppol*nsig_ab))
 allocate(Sr%sigxcme4sd(b1gw:b2gw,nkibz,nomega4sd,nsppol*nsig_ab))
 ! Mesh along the real axis.
 if (nomega_r>0) then
  allocate(Sr%omega_r(nomega_r))
  if (PRESENT(omega_r)) then
   ii=assert_eq(SIZE(omega_r),SIZE(Sr%omega_r),&
& 'DIM omega_r=/Sr%omega_r',"<stdin>",813)
   Sr%omega_r(:)=omega_r(:)
  end if
 end if
 ! === Analytical Continuation ===
 !if (mod10==1) then
  if (nomega_i>0) then
  ! FIXME omegasi should not be in Sp% here we should construct the mesh
  allocate(Sr%omega_i(nomega_i))
  !; Sr%omega_i=Sp%omegasi FIXME this has to be done outside
  allocate(Sr%sigcmesi (b1gw:b2gw,nkibz,nomega_i,nsppol*nsig_ab))
  allocate(Sr%sigxcmesi(b1gw:b2gw,nkibz,nomega_i,nsppol*nsig_ab))
  Sr%omega_i =czero
  Sr%sigcmesi =czero
  Sr%sigxcmesi=czero
  if (PRESENT(omega_i)) then
   ii=assert_eq(SIZE(Sr%omega_i),SIZE(omega_i),&
& 'DIM Sr%omega_i /= omega_i',"<stdin>",831)
   Sr%omega_i=omega_i
  end if
 end if
 Sr%e0 =zero
 Sr%egw =czero
 Sr%e0gap =zero
 Sr%sigcme =czero
 Sr%sigxme =czero
 Sr%sigxcme =czero
 Sr%sigcmee0 =czero
 Sr%ze0 =czero
 Sr%dsigmee0 =czero
 Sr%sigmee =czero
 Sr%omega4sd =czero
 Sr%sigcme4sd =czero
 Sr%sigxcme4sd=czero
 Sr%degw =czero
end subroutine allocate_sigma_results
subroutine etsf_dump_QP(Sr,QP_BSt,KS_BSt,Hdr,Cryst,Kmesh,filapp)
 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_numeric_tools, only : c2r
 use m_crystal, only : abi_crystal_put
 use m_electrons, only : abi_bands_put
 use interfaces_14_hidewrite
 use interfaces_16_hideleave
 use interfaces_61_ionetcdf
 implicit none
 type(Bandstructure_type),intent(in) :: QP_BSt,KS_BSt
 type(BZ_mesh_type),intent(in) :: Kmesh
 type(Crystal_structure),intent(in) :: Cryst
 type(Sigma_results),intent(in) :: Sr
 type(Hdr_type),intent(inout) :: Hdr
 character(len=fnlen),intent(in) :: filapp
 character(len=500) :: msg
 write(msg,'(4a)')ch10,&
& ' etsf_dump_QP : ERROR - ',ch10,&
& '  ETSF-IO support is not activated. '
 call wrtout(std_out,msg,'COLL')
 !??call leave_new('COLL')
end subroutine etsf_dump_QP
subroutine abi_etsf_get_QP(Sr,KS_BSt,Hdr,Cryst,filapp)
 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_crystal
 use m_numeric_tools, only : r2c
 use m_electrons, only : abi_bands_read
 use interfaces_14_hidewrite
 use interfaces_16_hideleave
 use interfaces_61_ionetcdf
 implicit none
 type(Bandstructure_type),intent(out) :: KS_BSt
 type(Crystal_structure),intent(out) :: Cryst
 type(Sigma_results),intent(out) :: Sr
 type(Hdr_type),intent(out) :: Hdr
 character(len=fnlen),intent(in) :: filapp
 character(len=500) :: msg
 write(msg,'(4a)')ch10,&
& ' abi_etsf_get_QP : ERROR - ',ch10,&
& '  ETSF-IO support is not activated. '
 call wrtout(std_out,msg,'COLL')
 !??call leave_new('COLL')
end subroutine abi_etsf_get_QP
END MODULE m_sigma_results


> the error message reported by XLF
I am not sure exactly which error message you mean. I gave you the output 
from make. How else do I generate error messages?

> as well as the version of the compiler.
The XLF version is 11.1

If the problem is indeed a compiler error, would using another version help? 
I have 10.1 at hand. This is a downgrade, but should I try it?



________________________________________
Von: Matteo Giantomassi [Matteo.Giantomassi@uclouvain.be]
Gesendet: Mittwoch, 1. Juli 2009 16:12
An: forum@abinit.org
Betreff: Re: AW: [abinit-forum] still problems with compilation on ibm power

> Hi Matteo,
>
> I did what you suggested and the compilation went a little further. But
> then I got the same kind of error for the file m_sigma_results. I did the
> same for this file, enclosing the executable part and commenting the
> calls. Is this ok or is m_sigma_results needed for something important?

 Dear Christian,

 It's ok: m_sigma_results is only needed in order to run GW calculations.
 The code should compile and the features of abinit you are interested in
 will be available despite the fact that you had to comment all the calls
to the
 procedure contained in the module.

> The compilation seemed to be fine, up to a point where I got
> ar: hdr_vs_dtset.o: No such file or directory         (not surprising,
> since it was not compiled)

 That's strange because the object file hdr_vs_dtset.o should be present
 in 68_gw once make has completed the compilation of the sources in 68_gw.

 Actually hdr_vs_dtset.o will contain a fake routine but it should be there
 to avoid problems at compile-time. What happens if you do the following:

 cd 68_gw && make clean && make


> make[3]: *** [lib68_gw.a] Error 1
>(I did not find this lib68_gw anywhere)
>

 The library lib68_gw  should be produced by ar (the command below) but
 ar exits with a non-zero status error due to
 the missing hdr_vs_dtset.o.

 BTW: It's difficult to understand the reason of your problem without
 having the preprocessed source files. I suspect your compilation
 problem is due to an internal compiler error and not to some kind
 of syntax error in the GW code. Could you provide the CPP preprocessed
 source for hdr_vs_dtset and m_sigma_results, the
 error message reported by XLF as well as the version of the compiler.

 Thanks in advance.

 Best Regards,
 Matteo Giantomassi


> Now what? The complete output of the compilation is
>
> 1501-510  Compilation successful for file split_sigc_cpp.f90.
> ABINIT WRAPPER END
> rm -f lib68_gw.a
> ar rc lib68_gw.a assemblychi0_sym.o assemblychi0sfq0.o
> assemblychi0q0_sym.o bands_classification.o bstruct_init.o calc_coh.o
> calc_density.o calc_ffm.o calc_rpa_functional.o calc_sig_noppm.o
> calc_sig_ppm.o calc_sig_ppm_comp.o calc_vHxc_braket.o calc_wf_qp.o cchi0.o
> cchi0q0.o cppm1par.o cppm2par.o cppm3par.o cppm4par.o cqratio.o csigme.o
> cutoff_m_elem.o debug_tools.o dosym.o dotproductqrc.o eps1_tc.o
> fourdp_6d.o get_bands_sym_GW.o gwcompleteness.o gw_tools.o hdr_init.o
> hdr_vs_dtset.o kramerskronig.o linear_optics_paw.o make_epsm1_driver.o
> mlwfovlp_qp.o m_bz_mesh.o m_commutator_vkbr.o m_coulombian.o m_crystal.o
> m_dyson_solver.o m_electrons.o m_finite_cylinder.o m_fft_mesh.o
> m_gsphere.o m_gwannier.o m_gwdefs.o m_io_screening.o m_little_group.o
> m_ppmodel.o m_qparticles.o m_screening.o m_sigma_results.o m_spectra.o
> m_wfs.o outeps.o pawr.o paw_mknewh0.o paw_mkrhox.o paw_mkrhox_spl.o
> paw_rho_tw_g.o paw_symcprj.o print_psps.o q0fit.o rho_tw_g.o setup_qmesh.o
> setup_screening.o setup_si
>  gma.o set_gwdistrb.o symmetrize_afm_chi0.o symf12.o spectral.o
> split_sigc.o interfaces_68_gw.o
> ar: hdr_vs_dtset.o: No such file or directory
> make[3]: *** [lib68_gw.a] Error 1
> make[3]: Leaving directory
> `/usr/local/abinit/abinit-5.8.3/build/src/68_gw'
> make[2]: *** [all-recursive] Error 1
> make[2]: Leaving directory `/usr/local/abinit/abinit-5.8.3/build/src'
> make[1]: *** [all-recursive] Error 1
> make[1]: Leaving directory `/usr/local/abinit/abinit-5.8.3/build'
> make: *** [all] Error 2