!{\src2tex{textfont=tt}}
!!****f* ABINIT/eig2tot
!! NAME
!! eig2tot
!!
!! FUNCTION
!! This routine calculates the second-order eigenvalues.
!! The output eig2nkq is this quantity for the input k points.
!!
!! COPYRIGHT
!! Copyright (C) 1999-2014 ABINIT group (PB, XG)
!! This file is distributed under the terms of the
!! GNU General Public License, see ~abinit/COPYING
!! or http://www.gnu.org/copyleft/gpl.txt .
!! For the initials of contributors, see ~abinit/doc/developers/contributors .
!!
!! INPUTS
!!  bdeigrf = number of bands for which to calculate the second-order eigenvalues.
!!  clflg(3,mpert)= array on calculated perturbations for eig2rf.
!!  dim_eig2nkq = 1 if eig2nkq is to be computed.
!!  eigbrd(2,mband*nsppol,nkpt,3,natom,3,natom) = broadening factors for the 
!!            electronic eigenvalues (optional).
!!  eigen0(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all K-points: 
!!            <k,n'|H(0)|k,n'> (hartree).
!!  eigenq(nkpt_rbz*mband*nsppol) = 0-order eigenvalues at all shifted K-points:
!!            <k+Q,n'|H(0)|k+Q,n'> (hartree).
!!  eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert) = matrix of first-order: 
!!            <k+Q,n'|H(1)|k,n> (hartree) (calculated in cgwf3).
!!  eig2nkq(2,mband*nsppol,nkpt,3,natom,3,natom*dim_eig2nkq) = second derivatives of
!!            the electronic eigenvalues.
!!  elph2_imagden = imaginary part of the denominator of the sum-over-state expression
!!            for the electronic eigenenergy shift due to second-order electron-phonon
!!            interation.
!!  ieig2rf = integer for calculation type.
!!  indsym(4,nsym,natom) = indirect indexing array for atom labels
!!            (not used yet, but will be used with symmetries).
!!  mband = maximum number of bands.
!!  mpert = maximum number of perturbations.
!!  natom = number of atoms in the unit cell.
!!  npert = number of phonon perturbations, without taking into account directions:
!!            natom. 
!!  nsym = number of symmetries (not used yet).
!!  mpi_enreg = informations about MPI parallelization.
!!  nkpt_rbz = number of k-points for each perturbation.
!!  nsppol = 1 for unpolarized, 2 for spin-polarized.
!!  smdelta = integer controling the calculation of electron lifetimes.
!!  symq(4,2,nsym) = 1 if symmetry preserves present qpoint. From symq3 (not used yet).
!!  symrec(3,3,nsym) = 3x3 matrices of the group symmetries (reciprocal space)
!!            (not used yet).
!!  symrel(3,3,nsym) = array containing the symmetries in real space (not used yet).
!!  timrev = 1 if time-reversal preserves the q wavevector; 0 otherwise 
!!            (not in use yet).
!!  dtset = OPTIONAL, dataset structure containing the input variable of the
!!            calculation. This is required to use the k-interpolation routine.
!!  eigenq_fine(mband_fine,mkpt_fine,nsppol_fine) = OPTIONAL, 0-order eigenvalues
!!            at all shifted K-points: <k+Q,n'|H(0)|k+Q,n'> (hartree) of the
!!            fine grid. This information is read from the WF dense k-grid file.  
!!  hdr_fine = OPTIONAL, header of the WF file of the fine k-point grid. This
!!            variable is required for the k-interpolation routine.  
!!  hdr0     = OPTIONAL, header of the GS WF file of the corse k-point grid. This
!!            variable is required for the k-interpolation routine.  
!!
!! OUTPUT
!!  eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= diagonal part of the 
!!            second-order eigenvalues: E^{(2),diag}_{k,q,j}.
!!  eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)= OPTIONAL, array containing the
!!            electron lifetimes.
!!
!! PARENTS
!!      respfn
!!
!! CHILDREN
!!      distrb2,ioddb8_out,kptfine_av,outbsd,outfan,psddb8,smeared_delta,timab
!!      xmpi_sum
!!
!! SOURCE

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

#include "abi_common.h"

subroutine eig2tot(dtfil,xred,psps,pawtab,displ,natom,gprimd,bdeigrf,clflg,dim_eig2nkq,eigen0,eigenq,eigen1,eig2nkq,&
&  elph2_imagden,esmear,ieig2rf,mband,mpert,npert,mpi_enreg,&
&  nkpt_rbz,nsppol,smdelta,phfrq,dtset,occ_rbz,eigbrd,eigenq_fine,hdr_fine,hdr0) 


 use defs_basis
 use defs_datatypes
 use defs_abitypes
 use m_profiling
 use m_xmpi
 use m_errors

 use m_pawtab,   only : pawtab_type
 use m_ddb,      only : psddb8

!This section has been created automatically by the script Abilint (TD).
!Do not modify the following lines by hand.
#undef ABI_FUNC
#define ABI_FUNC 'eig2tot'
 use interfaces_18_timing
 use interfaces_32_util
 use interfaces_51_manage_mpi
 use interfaces_72_response, except_this_one => eig2tot
!End of the abilint section

 implicit none

!Arguments ------------------------------------
!scalars
 integer,intent(in) :: bdeigrf,dim_eig2nkq,ieig2rf,mband,mpert,natom,nkpt_rbz
 integer,intent(in) :: npert,nsppol,smdelta
 real(dp),intent(in) :: elph2_imagden,esmear
 type(MPI_type),intent(inout) :: mpi_enreg
 type(datafiles_type), intent(in) :: dtfil
 type(pseudopotential_type), intent(inout) :: psps
!arrays
 real(dp) :: gprimd(3,3)
 integer,intent(in) :: clflg(3,mpert)
 real(dp),intent(in) :: displ(2*3*natom*3*natom)
 real(dp),intent(in) :: eigen0(nkpt_rbz*mband*nsppol)
 real(dp),intent(in) :: eigen1(nkpt_rbz*2*nsppol*mband**2,3,mpert)
 real(dp),intent(in) :: eigenq(nkpt_rbz*mband*nsppol)
 real(dp),intent(in) :: phfrq(mpert)
 real(dp),intent(in) :: occ_rbz(mband*nkpt_rbz*nsppol)
 real(dp),intent(inout) :: eig2nkq(2,mband*nsppol,nkpt_rbz,3,npert,3,npert*dim_eig2nkq)
 type(dataset_type), intent(in) :: dtset
 real(dp),intent(in) :: xred(3,natom)
 real(dp),intent(inout),optional :: eigbrd(2,mband*nsppol,nkpt_rbz,3,npert,3,npert)
 real(dp),intent(in),pointer,optional :: eigenq_fine(:,:,:)
 type(pawtab_type), intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
 type(hdr_type),intent(in),optional :: hdr_fine,hdr0

!Local variables-------------------------------
!tolerance for non degenerated levels
!scalars
 integer :: band2tot_index,band_index,bandtot_index,choice,fullinit,iband,idir1,idir2
 integer :: iatom1,iatom2,ikpt,imode,ipert1,ipert2,isppol,jband,nkpt_sub,ikpt2,nblok,unitout,vrsddb
!integer :: ipw
 character(len=fnlen) :: dscrpt
 integer :: master,me,spaceworld,ierr
! real(dp),parameter :: etol=1.0d-6
 real(dp),parameter :: etol=1.0d-7
!real(dp),parameter :: etol=zero 
 real(dp) :: ar,ai,deltae,den,eig1_i1,eig1_i2,eigen_corr
 real(dp) :: eig1_r1,eig1_r2,eig2_diai,den_av
 real(dp) :: eig2_diar,eigbrd_i,eigbrd_r,tolwfr,wgt_int
 character(len=500) :: message
 logical :: test_do_band
!arrays
 integer, allocatable :: nband_rbz(:)
 integer,pointer      :: kpt_fine_sub(:)
 real(dp)             :: tsec(2)
 real(dp),allocatable :: center(:)
 real(dp) :: eigen(mband*nsppol),eigen_prime(mband*nsppol)
 real(dp),allocatable :: fan(:,:,:,:,:)
 real(dp),allocatable :: smdfun(:,:)
 real(dp),pointer     :: wgt_sub(:)
 real(dp) :: d_at1_dir1_re,d_at1_dir1_im
 real(dp) :: d_at2_dir2_re,d_at2_dir2_im
 real(dp) :: displ2cart(2,3,3),displ2red(2,3,3),tmp_displ2(2,3,3)

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

!Init parallelism
 master =0
 spaceworld=mpi_enreg%comm_cell
 me=mpi_enreg%me_kpt
!DEBUG
!write(std_out,*)' eig2tot : enter '
!write(std_out,*)' mband=',mband
!write(std_out,*)' nsppol=',nsppol
!write(std_out,*)' nkpt_rbz=',nkpt_rbz
!write(std_out,*)' npert=',npert
!ENDDEBUG

!Init interpolation method
 if(present(eigenq_fine))then
   ABI_ALLOCATE(center,(3))
 end if

 call timab(148,1,tsec)

 if(nsppol==2)then
   message = 'nsppol=2 is not yet treated in eig2tot. Sorry for this ...'
   MSG_ERROR(message)
 end if        

 band2tot_index =0
 bandtot_index=0
 band_index=0

 if(xmpi_paral==1) then
   ABI_ALLOCATE(mpi_enreg%proc_distrb,(nkpt_rbz,mband,nsppol))
   ABI_ALLOCATE(nband_rbz,(nkpt_rbz))
   if (allocated(mpi_enreg%my_kpttab)) then
     ABI_DEALLOCATE(mpi_enreg%my_kpttab)
   end if
   ABI_ALLOCATE(mpi_enreg%my_kpttab,(nkpt_rbz))
!  Assume the number of bands is the same for all k points.
   nband_rbz(:)=mband
   call distrb2(mband,nband_rbz,nkpt_rbz,mpi_enreg%nproc_cell,nsppol,mpi_enreg)
 end if

 if(ieig2rf == 4 ) then
   ABI_ALLOCATE(fan,(2,mband*nsppol,nkpt_rbz,mband*nsppol,3*natom))
   fan(:,:,:,:,:) = zero
   eigbrd(:,:,:,:,:,:,:)=zero
   eigen_corr = 0
 end if
 
 do isppol=1,nsppol
   do ikpt =1,nkpt_rbz

     if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,mband,isppol,me)) then
       band2tot_index = band2tot_index + 2*mband**2
       bandtot_index = bandtot_index + mband
       cycle
     end if

     if(present(eigenq_fine))then
       write(std_out,*) 'Start of the energy denominator interpolation method.'
       nkpt_sub = 0
!      center is the k+q point around which we will average the kpt_fine
       center = hdr0%kptns(:,ikpt)+ dtset%qptn(:) 

       call kptfine_av(center,dtset%qptrlatt,hdr_fine%kptns,hdr_fine%nkpt,&
&       kpt_fine_sub,nkpt_sub,wgt_sub)
       write(std_out,'(a,3f8.4,a,i3)') 'Number of k-points of the fine grid &
&       around the k+Q point ',center,' is:',nkpt_sub
       write(std_out,'(a,f10.5)') 'The sum of the weights of the k-points is: ',SUM(wgt_sub)
     end if

     if(smdelta >0) then   !broadening
       if(.not.allocated(smdfun))  then
         ABI_ALLOCATE(smdfun,(mband,mband))
       end if
       smdfun(:,:) = zero
       do iband=1,mband
         eigen(iband) = eigen0(iband+bandtot_index)
         eigen_prime(iband) =eigenq(iband+bandtot_index)
       end do
       if(esmear>tol6) then
         call smeared_delta(eigen,eigen_prime,esmear,mband,smdelta,smdfun)
       end if
     end if

     ipert1=1 ! Suppose all perturbations lead to the same number of planewaves

     do iband=1,bdeigrf

!      If the k point and band belong to me, compute the contribution
       test_do_band=.true.
       if(mpi_enreg%proc_distrb(ikpt,iband,isppol)/=me)test_do_band=.false.
       
       if(test_do_band)then
!        ------------------------------------------------------------------------------------------------------!
!        ------- ieig2rf ==3 : Non dynamic traditional AHC theory with Sternheimer (computed in eig2stern.F90)-!
!        ------------------------------------------------------------------------------------------------------!
         if(ieig2rf==3 .or. ieig2rf==4 ) then
           do ipert1=1,npert
             do idir1=1,3
               if(clflg(idir1,ipert1)==0) cycle
               do ipert2=1,npert
                 do idir2=1,3
                   if(clflg(idir2,ipert2)==0)cycle
                   eig2_diar = zero ; eig2_diai = zero ; eigbrd_r = zero ; eigbrd_i = zero
                   do jband=1,mband
                     eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
                     eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,ipert2)
                     eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,ipert1)
                     eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,ipert2) !the negative sign is from the CC
!                    If no interpolation, fallback on to the previous
!                    implementation
                     if(.not. present(eigenq_fine))then
                       deltae=eigenq(jband+bandtot_index)-eigen0(iband+bandtot_index)
                     end if
                     ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
                     ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2

!                    Sum over all active space to retrieve the diagonal gauge
!                    if(abs(deltae)>etol) then ! This is commented because
!                    there is no problem with divergencies with elph2_imag != 0
                     if( present(eigenq_fine))then
                       den_av = zero
                       wgt_int = zero
                       do ikpt2=1,nkpt_sub
                         deltae=eigenq_fine(jband,kpt_fine_sub(ikpt2),1)&
&                         -eigen0(iband+bandtot_index)
                         den_av = den_av-(wgt_sub(ikpt2)*deltae)/(deltae**2+elph2_imagden**2)
                         wgt_int = wgt_int+wgt_sub(ikpt2)
                       end do
                       den = den_av/wgt_int
                     else
                       if(abs(elph2_imagden) < etol) then  
                         if(abs(deltae)>etol) then
                           den=-one/(deltae**2+elph2_imagden**2)
                         else
                           den= zero
                         end if
                       else
                         den=-one/(deltae**2+elph2_imagden**2)
                       end if
                     end if
                     
                     if( present(eigenq_fine))then
                       eig2_diar=eig2_diar+ar*den
                       eig2_diai=eig2_diai+ai*den
                     else
                       eig2_diar=eig2_diar+ar*deltae*den
                       eig2_diai=eig2_diai+ai*deltae*den
                     end if

                     if(present(eigbrd))then
                       if(smdelta >0) then   !broadening
                         eigbrd_r = eigbrd_r + ar*smdfun(iband,jband)
                         eigbrd_i = eigbrd_i + ai*smdfun(iband,jband)
                       end if
                     end if
                   end do !jband

!                  Store the contribution
                   eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = &
&                   eig2nkq(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diar
                   eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = &
&                   eig2nkq(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) + eig2_diai 

                   if(present(eigbrd))then
                     if(smdelta >0) then   !broadening
                       eigbrd(1,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_r
                       eigbrd(2,iband+band_index,ikpt,idir1,ipert1,idir2,ipert2) = eigbrd_i
                     end if
                   end if

                 end do !idir2
               end do !ipert2
             end do  !idir1
           end do   !ipert1
         end if !ieig2rf 3

!        -------------------------------------------------------------------------------------------!
!        ------- ieig2rf ==4  Dynamic AHC using second quantization and Sternheimer from eig2stern -!
!        -------------------------------------------------------------------------------------------!   
         if(ieig2rf ==4 ) then
           do imode=1,3*natom
             if(phfrq(imode)>etol)then
               eigbrd_r = zero ; eigbrd_i = zero
               do jband=1,mband 
                 do iatom1=1,natom
                   do iatom2=1,natom
                     do idir1=1,3
                       do idir2=1,3
!                        Compute the mean cartesian displacements
                         d_at1_dir1_re=displ(1 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1))))
                         d_at1_dir1_im=displ(2 + 2*(idir1-1 +3*(iatom1-1 +natom*(imode-1))))
                         d_at2_dir2_re=displ(1 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1))))
                         d_at2_dir2_im=displ(2 + 2*(idir2-1 +3*(iatom2-1 +natom*(imode-1))))
                         displ2cart(1,idir1,idir2)=(d_at1_dir1_re*d_at2_dir2_re+ &
&                         d_at1_dir1_im*d_at2_dir2_im )
                         displ2cart(2,idir1,idir2)=(d_at1_dir1_im*d_at2_dir2_re- &
&                         d_at1_dir1_re*d_at2_dir2_im )
!                        DBSP
!                        if dtset%qptn(1) < etol .and. dtset%qptn(2) < etol .and. dtset%qptn(3) < etol
!                        END
                       end do
                     end do
!                    Switch to reduced coordinates in two steps
                     tmp_displ2(:,:,:)=zero
                     do idir1=1,3
                       do idir2=1,3
                         tmp_displ2(:,:,idir1)=tmp_displ2(:,:,idir1)+displ2cart(:,:,idir2)*gprimd(idir2,idir1)
                       end do
                     end do
                     displ2red(:,:,:)=zero
                     do idir1=1,3
                       do idir2=1,3
                         displ2red(:,idir1,:)=displ2red(:,idir1,:)+tmp_displ2(:,idir2,:)*gprimd(idir2,idir1)
                       end do
                     end do
                     do idir1=1,3
                       do idir2=1,3
                         eig1_r1 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir1,iatom1)
                         eig1_r2 = eigen1(2*jband-1+(iband-1)*2*mband+band2tot_index,idir2,iatom2)
                         eig1_i1 = eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir1,iatom1)
                         eig1_i2 = - eigen1(2*jband+(iband-1)*2*mband+band2tot_index,idir2,iatom2)
                         ar=eig1_r1*eig1_r2-eig1_i1*eig1_i2
                         ai=eig1_r1*eig1_i2+eig1_i1*eig1_r2
                         
                         fan(1,iband+band_index,ikpt,jband,imode) = fan(1,iband+band_index,ikpt,jband,imode)-&
&                         (ar*displ2red(1,idir1,idir2)-ai*displ2red(2,idir1,idir2))
                         fan(2,iband+band_index,ikpt,jband,imode) = fan(2,iband+band_index,ikpt,jband,imode)-&
&                         (ar*displ2red(2,idir1,idir2)+ai*displ2red(1,idir1,idir2))
                         
                       end do !idir2
                     end do !idir1
                   end do ! iatom2
                 end do ! iatom1
               end do !jband
             end if ! phfrq(imode)
           end do !imode
         end if !ieig2rf 4
       end if ! Selection of processor
     end do !iband

     band2tot_index = band2tot_index + 2*mband**2
     bandtot_index = bandtot_index + mband

     if(present(eigenq_fine))then
       ABI_DEALLOCATE(kpt_fine_sub) ! Deallocate the variable
       ABI_DEALLOCATE(wgt_sub)
     end if

   end do    !ikpt
   band_index = band_index + mband
 end do !isppol

!Accumulate eig2nkq and/or eigbrd
 if(xmpi_paral==1) then
   if(ieig2rf == 3) then
     call xmpi_sum(eig2nkq,spaceworld,ierr)
   end if
   if(ieig2rf == 4) then
     call xmpi_sum(eig2nkq,spaceworld,ierr)
     call xmpi_sum(fan,spaceworld,ierr)
   end if
   if(present(eigbrd) .and. (ieig2rf == 3 .or. ieig2rf == 4))then
     if(smdelta >0) then
       call xmpi_sum(eigbrd,spaceworld,ierr)
     end if
   end if
   ABI_DEALLOCATE(nband_rbz)
   ABI_DEALLOCATE(mpi_enreg%proc_distrb)
   ABI_DEALLOCATE(mpi_enreg%my_kpttab)
 end if 

 if(ieig2rf > 2) then
   write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGR2D.'
   write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
   do idir1=1,3
     do idir2=1,3
       ar=eig2nkq(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
       ai=eig2nkq(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
       write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai 
     end do ! idir2
   end do ! idir1
 end if 

 if(present(eigbrd))then
   if(smdelta >0) then   !broadening
     write(ab_out,'(a)')' '
     write(ab_out,'(a)')' Components of second-order derivatives of the electronic energy, EIGI2D.'
     write(ab_out,'(a)')' For automatic tests, printing the matrix for the first k-point, first band, first atom.'
     do idir1=1,3
       do idir2=1,3
         ar=eigbrd(1,1,1,idir1,1,idir2,1) ; if(abs(ar)<tol10)ar=zero
         ai=eigbrd(2,1,1,idir1,1,idir2,1) ; if(abs(ai)<tol10)ai=zero
         write (ab_out,'(4i4,2es20.10)') idir1,1,idir2,1,ar,ai
       end do
     end do !nband
   end if
 end if

 if(allocated(smdfun))  then
   ABI_DEALLOCATE(smdfun)
 end if
 if(present(eigenq_fine))then
   ABI_DEALLOCATE(center)
 end if

!print _EIGR2D file for this perturbation

 unitout = dtfil%unddb
 vrsddb=100401
 dscrpt=' Note : temporary (transfer) database '
!tolwfr must be initialized here, but it is a dummy value
 tolwfr=1.0_dp
 call ioddb8_out (dscrpt,dtfil%fnameabo_eigr2d,dtset%natom,dtset%mband,&
& dtset%nkpt,dtset%nsym,dtset%ntypat,dtfil%unddb,vrsddb,&
& dtset%acell_orig(1:3,1),dtset%amu_orig(:,1),dtset%dilatmx,dtset%ecut,dtset%ecutsm,&
& dtset%intxc,dtset%iscf,dtset%ixc,dtset%kpt,dtset%kptnrm,&
& dtset%natom,dtset%nband,dtset%ngfft,dtset%nkpt,dtset%nspden,dtset%nspinor,&
& dtset%nsppol,dtset%nsym,dtset%ntypat,occ_rbz,dtset%occopt,dtset%pawecutdg,&
& dtset%rprim_orig(1:3,1:3,1),dtset%sciss,dtset%spinat,dtset%symafm,dtset%symrel,&
& dtset%tnons,tolwfr,dtset%tphysel,dtset%tsmear,&
& dtset%typat,dtset%usepaw,dtset%wtk,xred,psps%ziontypat,dtset%znucl)

 nblok=1 ; fullinit=1 ; choice=2
 call psddb8 (choice,psps%dimekb,psps%ekb,fullinit,psps%indlmn,&
& psps%lmnmax,nblok,dtset%ntypat,dtfil%unddb,pawtab,&
& psps%pspso,psps%usepaw,psps%useylm,vrsddb)

 call outbsd(bdeigrf,dtset,eig2nkq,dtset%natom,nkpt_rbz,unitout)

!print _FAN file for this perturbation
 if(ieig2rf == 4 ) then
   unitout = dtfil%unddb
   vrsddb=100401
   dscrpt=' Note : temporary (transfer) database '
!  tolwfr must be initialized here, but it is a dummy value
   tolwfr=1.0_dp
   call ioddb8_out (dscrpt,dtfil%fnameabo_fan,dtset%natom,dtset%mband,&
&   dtset%nkpt,dtset%nsym,dtset%ntypat,dtfil%unddb,vrsddb,&
&   dtset%acell_orig(1:3,1),dtset%amu_orig(:,1),dtset%dilatmx,dtset%ecut,dtset%ecutsm,&
&   dtset%intxc,dtset%iscf,dtset%ixc,dtset%kpt,dtset%kptnrm,&
&   dtset%natom,dtset%nband,dtset%ngfft,dtset%nkpt,dtset%nspden,dtset%nspinor,&
&   dtset%nsppol,dtset%nsym,dtset%ntypat,occ_rbz,dtset%occopt,dtset%pawecutdg,&
&   dtset%rprim_orig(1:3,1:3,1),dtset%sciss,dtset%spinat,dtset%symafm,dtset%symrel,&
&   dtset%tnons,tolwfr,dtset%tphysel,dtset%tsmear,&
&   dtset%typat,dtset%usepaw,dtset%wtk,xred,psps%ziontypat,dtset%znucl)

   nblok=1 ; fullinit=1 ; choice=2
   call psddb8 (choice,psps%dimekb,psps%ekb,fullinit,psps%indlmn,&
&   psps%lmnmax,nblok,dtset%ntypat,dtfil%unddb,pawtab,&
&   psps%pspso,psps%usepaw,psps%useylm,vrsddb)
   
   call outfan(bdeigrf,dtset,fan,nkpt_rbz,unitout)

   ABI_DEALLOCATE(fan)
 end if
!print _EIGI2D file for this perturbation
 if(smdelta>0) then
   unitout = dtfil%unddb
   vrsddb=100401
   dscrpt=' Note : temporary (transfer) database '
!  tolwfr must be initialized here, but it is a dummy value
   tolwfr=1.0_dp
   call ioddb8_out (dscrpt,dtfil%fnameabo_eigi2d,dtset%natom,dtset%mband,&
&   dtset%nkpt,dtset%nsym,dtset%ntypat,dtfil%unddb,vrsddb,&
&   dtset%acell_orig(1:3,1),dtset%amu_orig(:,1),dtset%dilatmx,dtset%ecut,dtset%ecutsm,&
&   dtset%intxc,dtset%iscf,dtset%ixc,dtset%kpt,dtset%kptnrm,&
&   dtset%natom,dtset%nband,dtset%ngfft,dtset%nkpt,dtset%nspden,dtset%nspinor,&
&   dtset%nsppol,dtset%nsym,dtset%ntypat,occ_rbz,dtset%occopt,dtset%pawecutdg,&
&   dtset%rprim_orig(1:3,1:3,1),dtset%sciss,dtset%spinat,dtset%symafm,dtset%symrel,&
&   dtset%tnons,tolwfr,dtset%tphysel,dtset%tsmear,&
&   dtset%typat,dtset%usepaw,dtset%wtk,xred,psps%ziontypat,dtset%znucl)

   nblok=1 ; fullinit=1 ; choice=2
   call psddb8 (choice,psps%dimekb,psps%ekb,fullinit,psps%indlmn,&
&   psps%lmnmax,nblok,dtset%ntypat,dtfil%unddb,pawtab,&
&   psps%pspso,psps%usepaw,psps%useylm,vrsddb)

   call outbsd(bdeigrf,dtset,eigbrd,dtset%natom,nkpt_rbz,unitout)
 end if !smdelta


 call timab(148,2,tsec)

!DEBUG
!write(std_out,*)' eig2tot: exit'
!ENDDEBUG

end subroutine eig2tot
!!***

