! FDMNES subroutines

! Filling of FDM matrix and solving of the linear system of equation.

subroutine Mat_FDM(Adimp,Ampl_dft,Atom_axe,Base_hexa,Basereelt,Cal_xanes,cgrad, &
                      clapl,Dist_coop,distai,E_comp,Ecinetic_out,Eclie_out,Eimag,Eneg,Enervide, &
                      gradvr,ia_coop,iaabsi,iaprotoi,iato,ibord,icheck,ie,igrph,irep_util,isbord, &
                      iso,ispinin,isrt,isvois,ivois,karact,Lato,Lmaxa,lmaxso,Lso,mato,mpinodes,mpirank0,mso, &
                      n_atom_coop,n_atom_proto,nab_coop,natome,nbm,nbm_c,nbord,nbordf,nbtm, &
                      ngrph,nim,nicm,nlm_rixs,nlmagm,nlmmax,nlmomax,nlmsa,nlmsam,nlmso_reg,nphiato,nphiato_c,npoint,npr,npsom, &
                      nsort,nsortf,nso1,nspin,nspino,nspinp,numia,nvois,Phiato,Phiato_c,poidsa,poidso,Posi,R_rydb, &
                      Recop,Relativiste,Repres_comp,RIXS_gen,Rmtg,Rot_Spinor,Rsort,rvol,Rydberg,Spin_axis_free,Spinorbite, &
                      State_all_r,Sym_cubic,Tau_coop,Taull,Time_fill,Time_tria,V0bd_out,Vr,xyz,Ylm_comp, &
                      Ylmato,Ylms,Ylmso)

  use declarations
  implicit none
  
  integer:: i, ia, iaabsi, iab, ib, ich, icheck, ie, igrph, ii, ipra, iprb, irep, is1, is2, isg, isp, ispinin, j, &
    jj, jsp, L, L1, L2, Lm, Lm01, Lm01c, Lm02, Lm02c, Lm1, Lm2, Lmax_rixs, Lmaxso, Lmf, Lms, Lmw, m, m1, m2, mpierr, &
    MPI_host_num_for_mumps, mpinodes, mpirank0, n_atom_coop, n_atom_proto, &
    nab_coop, natome, nbm, nbm_c, nbtm, ngrph, nim, nicm, nligne, nligne_i, nligneso, nlm_rixs, &
    nlmagm, nlmmax, nlmomax, nlmsam, nlmso, nlmso_eq, nlmso_i, nlmso_reg, nlmw, nphiato, nphiato_c,npoint, npr, npsom, &
    nsort, nsort_c, nsort_r, nsortf, nso1, nspin, nspino, nspinp, nspinr, nvois

  integer, dimension(0:npoint):: new
  integer, dimension(n_atom_coop):: ia_coop
  integer, dimension(natome):: ianew, ianewL, iaprotoi, Lmaxa, nbord, nbordf, nlmsa
  integer, dimension(nsort):: isrt
  integer, dimension(npsom):: numia
  integer, dimension(nso1,ngrph):: Lso, mso, iso
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato
  integer, dimension(nbtm,natome):: ibord, isbord
  integer, dimension(nrepm,2):: irep_util
  integer, dimension(npsom,nvois):: ivois, isvois
  integer, dimension(:), allocatable:: lb1, lb2, newinv

  character(len=Length_name):: Matrix_info

  complex(kind=db):: ampl1, ampl2, cfac
  complex(kind=db), dimension(nopsm,nrepm):: karact
  complex(kind=db), dimension(nlm_rixs,nspinp,nlmomax,nspinp,0:mpinodes-1):: Ampl_dft
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,natome):: Taull
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,nab_coop,nspino):: Tau_coop
  complex(kind=db), dimension(nphiato_c,nlmagm,nspinp,nspino,natome):: Phiato_c
  complex(kind=db), dimension(nbm_c,nlmmax,nspinp,nspino,natome):: Ylms

  complex(kind=db), dimension(:,:), allocatable:: norm, norm_t
  complex(kind=db), dimension(:,:,:), allocatable:: Bessel, Neuman
  complex(kind=db), dimension(:,:,:,:), allocatable:: taull_tem

  logical:: Base_hexa, Basereel, Basereelt, Cal_xanes, Cal_comp, E_comp, Eneg, Recop, &
    Relativiste, Repres_comp, RIXS_gen, Rydberg, Spin_axis_free, Spinorbite, State_all_r, Stop_job, Sym_cubic, Ylm_comp
  logical, dimension(natome):: Atom_axe

  real(kind=db):: adimp, Dist, Eclie_out, Eimag, Enervide, R_rydb, Rsort, Time_fill, Time_tria 

  real(kind=db), dimension(2):: Dist_coop
  real(kind=db), dimension(3,3):: Rot_Spinor 
  real(kind=db), dimension(nspin):: Ecinetic_out, V0bd_out
  real(kind=db), dimension(nopsm,nspino):: Kar, Kari
  real(kind=db), dimension(nvois):: cgrad
  real(kind=db), dimension(0:nvois):: clapl
  real(kind=db), dimension(npoint,nspin):: Vr
  real(kind=db), dimension(nbm,natome):: poidsa
  real(kind=db), dimension(nim):: rvol
  real(kind=db), dimension(nsortf):: poidso
  real(kind=db), dimension(4,npsom):: xyz
  real(kind=db), dimension(nsort,nlmomax):: Ylmso
  real(kind=db), dimension(0:n_atom_proto):: Rmtg
  real(kind=db), dimension(natome):: distai
  real(kind=db), dimension(3,natome):: posi
  real(kind=db), dimension(nicm,3,nspin):: gradvr
  real(kind=db), dimension(nbtm,nlmmax,natome):: Ylmato
  real(kind=db), dimension(nphiato,nlmagm,nspinp,nspino,natome):: Phiato

  real(kind=db), dimension(:,:), allocatable :: smi, smr
  real(kind=db), dimension(:,:,:), allocatable :: Besselr, Neumanr

  call MPI_Comm_Size(MPI_COMM_MUMPS,MPI_host_num_for_mumps,mpierr)

  Stop_job = .true.

  do isp = 1,nspino
    irep = abs( irep_util(igrph,isp) )
! There is no up-down coupling in the following case
    if( irep == 0 ) irep = abs( irep_util(igrph,3-isp) )
    Kar(:,isp) = real( karact(:,irep), db )
    if( Repres_comp ) then
      if( irep_util(igrph,isp) > 0 ) then
        Kari(:,isp) = aimag( karact(:,irep) )
      else
        Kari(:,isp) = - aimag( karact(:,irep) )
      endif
    else
      Kari(:,isp) = 0._db
    endif
  end do
  Basereel = Basereelt
  if( E_comp .or. Ylm_comp ) then
    Cal_comp = .true.
    Basereel = .false.
  else
    Cal_comp = .false.
  endif

  if( .not. Basereel ) Cal_comp = .true.

  if( E_comp ) then
    nsort_c = nsort
    nsort_r = 0
  else
    nsort_c = 0
    nsort_r = nsort
  endif
  
  if( Spinorbite ) then
    nspinr = nspin
  else
    nspinr = 1
  endif
  
  allocate( Bessel(nsort_c,0:Lmaxso,nspinr) )
  allocate( Neuman(nsort_c,0:Lmaxso,nspinr) )
  allocate( Besselr(nsort_r,0:Lmaxso,nspinr) )
  allocate( Neumanr(nsort_r,0:Lmaxso,nspinr) )

  jsp = 0
  do isp = 1,nspin
    if( .not. Spinorbite .and. isp /= ispinin ) cycle
    jsp = jsp + 1
! Calculation of the wave function in the outer sphere
    call phiso(Adimp,Bessel,Besselr,E_comp,Ecinetic_out(isp),Eclie_out,Eimag, &
       Eneg,Enervide,icheck,jsp,isrt,Lmaxso,mpirank0,Neuman,Neumanr,npsom,nsort,nsort_c,nsort_r,nspinr, &
       R_Rydb,Rsort,Rydberg,V0bd_out(isp),xyz)
  end do

  nlmso_eq = nlmso_reg
  nlmso = nlmso_reg
  nligne = nspino * npr + sum( nlmsa(1:natome) ) + nlmso

  allocate( lb1(nligne) ); allocate( lb2(nligne) )
  allocate( newinv(nligne) )

  ich = icheck
  call newind(distai,ianew,ianewL,ibord,ich,isrt,ivois,lb1,lb2,mpirank0,natome, &
              nbordf,nbtm,new,newinv,nligne,nlmsa,nlmso_eq,npoint,nsort,nsortf,nspino,npsom,numia, &
              nvois,xyz)

  if( Cal_comp ) then
    nligne_i = nligne
    nlmso_i = nlmso
  else
    nligne_i = 0
    nlmso_i = 0
  endif
  allocate( smi(nlmso_i,nligne_i) )
  allocate( smr(nlmso,nligne) )
  smi(:,:) = 0._db
  smr(:,:) = 0._db

  if( icheck > 1 .or. icheck > 2 ) write(3,110) igrph

  nligneso = nligne - nlmso_eq

  if( nlmso > 0 ) &
    call Mat_solve(Base_hexa, Basereel, Bessel, Besselr, Cal_comp, cgrad, clapl, E_comp, Eimag, Enervide, &
      gradvr, ianew, ianewL, iato, ibord, icheck, ie, igrph, ii, isbord, iso, ispinin, isrt, &
      isvois, ivois, Kar, Kari, Lato, lb1, lb2, lmaxso, lso, mato, MPI_host_num_for_mumps, mpirank0, mso, &
      natome, nbm, nbm_c, nbord, nbordf, nbtm, Neuman, Neumanr, new, newinv, &
      ngrph, nicm, nim, nligne, nligne_i, nligneso, nlmsam, nlmagm, nlmmax, nlmomax, nlmsa, nlmso, nlmso_i, &
      nphiato, nphiato_c, npoint, npsom, nso1, nsort, nsort_c, nsort_r, nsortf, nspin, nspino, nspinp, nspinr, &
      numia, nvois, Phiato, Phiato_c, poidsa, poidso, Relativiste, Repres_comp, rvol, smi, smr,  Rot_Spinor, Spin_axis_free, &
      Spinorbite, Time_fill, Time_tria, Vr, Ylm_comp, Ylmato, Ylms, Ylmso)

  if( icheck > 2 ) then
    write(3,120)
    do ii = 1,nligne
      if( Cal_comp ) then
        write(3,130) ii, ( smr(Lms,ii), smi(Lms,ii), Lms = 1,nlmso )
      else
        write(3,140) ii, smr(1:nlmso,ii)
      endif
    end do
  endif

  deallocate( Bessel, Besselr, Neuman, Neumanr )

  if( Spinorbite ) then
    is1 = 1
    is2 = 2
  else
    is1 = ispinin
    is2 = ispinin
  endif

  ia = iaabsi

  allocate( norm(nlmso_reg,nlmso_reg) )
  allocate( norm_t(nlmso_reg,nlmso_reg) )

  if( Basereel .and. nlmso_reg > 0  ) then
! When simulation is done using real basis, one must renormalize.

    do Lm = 1,nlmso_reg
      ii = nligneso + Lm
      if( Cal_comp ) then
        do Lm2 = 1,nlmso_reg
          norm(Lm2,Lm) = cmplx( -smi(Lm2,ii), smr(Lm2,ii), db )
        end do
      else
        do Lm2 = 1,nlmso_reg
          norm(Lm2,Lm) = cmplx( 0._db, smr(Lm2,ii), db )
        end do
      endif
    end do
    do Lm = 1,nlmso_reg
      norm(Lm,Lm) = 1._db + norm(Lm,Lm)
    end do
 
    if( icheck > 2 ) then
      write(3,'(/A)') '   Lm  norm(Lm,1..nlmso_reg)'
      do Lm = 1,nlmso_reg
        write(3,130) Lm, norm(Lm,:)
      end do
    endif

! Inversion
    Matrix_info = ' Norm matrix from Mat_FDM'
    call invcomp(Matrix_info,nlmso_reg,norm,nlmso_reg,nlmso_reg,0,Stop_job)
! for speedup access to norm
    norm_t = TRANSPOSE(norm)
  endif

  do ia = 1,natome

    if( ia /= iaabsi .and. .not. State_all_r .and. Cal_xanes ) cycle

    do Lm1 = 1,nlmsa(ia)
      
      L1 = Lato(Lm1,ia,igrph)
      m1 = mato(Lm1,ia,igrph)
      if( Spinorbite ) then
        is1 = iato(Lm1,ia,igrph)
      else
        is1 = ispinin
      endif
      Lm01 = L1**2 + L1 + 1 + m1

      ii = ianew(ia) + Lm1

      do Lm2 = 1,nlmsa(ia)

        L2 = Lato(Lm2,ia,igrph)
        m2 = mato(Lm2,ia,igrph)
        if( nspino == 2 ) then
          is2 = iato(Lm2,ia,igrph)
        else
          is2 = ispinin
        endif
        Lm02 = L2**2 + L2 + 1 + m2

        jj = ianew(ia) + Lm2

        cfac = ( 0._db, 0._db )

        do Lmf = 1,nlmso_reg

          if( Basereel ) then
            if( Cal_comp ) then
              ampl1 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,ii), smi(1:nlmso_reg,ii), db ))
              ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,jj), smi(1:nlmso_reg,jj), db ))
            else
              ampl1 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,ii) )
              ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,jj) )
            endif
          else
            ampl1 = cmplx(smr(Lmf,ii),smi(Lmf,ii),db)
            ampl2 = cmplx(smr(Lmf,jj),smi(Lmf,jj),db)
          endif

 ! One multiplies by -img in order -aimag(Taull) is the density of state
          cfac = cfac - img * ampl1 * conjg(  ampl2 )
        end do

        Taull(Lm01,is1,Lm02,is2,ia) = Taull(Lm01,is1,Lm02,is2,ia) + cfac

      end do
    end do

  end do

! Conjugate representation
  if( Repres_comp .and. .not. Spinorbite ) then
    do ia = 1,natome

      do Lm1 = 1,nlmsa(ia)
        L1 = Lato(Lm1,ia,igrph)
        m1 = mato(Lm1,ia,igrph)
        Lm01 = L1**2 + L1 + 1 + m1
        Lm01c = Lm01 - 2 * m1
        
        is1 = ispinin

        do Lm2 = 1,nlmsa(ia)

          L2 = Lato(Lm2,ia,igrph)
          m2 = mato(Lm2,ia,igrph)
          
          Lm02 = L2**2 + L2 + 1 + m2
          Lm02c = Lm02 - 2 * m2
          is2 = ispinin
          isg = (-1)**(m1+m2)
          Taull(Lm02c,is1,Lm01c,is2,ia) = Taull(Lm02c,is1,Lm01c,is2,ia) + isg * Taull(Lm01,is1,Lm02,is2,ia)
          
        end do
      end do
    end do

    if( icheck > 1 ) call Write_ampl(Basereel,Cal_comp,iaabsi,ianew,iato,igrph,irep_util,iso,ispinin,Lato,Lso, &
                    mato,mso,natome,newinv,ngrph,nim,nligne,nlmsa,nlmsam,nlmso_reg,norm,norm_t,npsom, &
                    nso1,nspino,rvol,smi,smr,Spinorbite,xyz)

  endif
  
! Atomic amplitude used for the RIXS formula
  if( RIXS_gen .and. Cal_xanes ) then
  
    Lmax_rixs = nint( sqrt( real( nlm_rixs ) ) ) - 1
    ia = iaabsi

    do Lm1 = 1,nlmsa(ia)
      ii = ianew(ia) + Lm1
      L1 = Lato(Lm1,ia,igrph)
      if( L1 > Lmax_rixs ) exit
      m1 = mato(Lm1,ia,igrph)
      if( Spinorbite ) then
        is1 = iato(Lm1,ia,igrph)
      else
        is1 = ispinin
      endif
      Lm01 = L1**2 + L1 + 1 + m1

      do Lmf = 1,nlmso_reg
        L2 = Lso(Lmf,igrph)
        m2 = mso(Lmf,igrph)
        Lm02 = L2**2 + L2 + 1 + mso(Lmf,igrph)
        if( Spinorbite ) then
          is2 = iso(Lmf,igrph)
        else
          is2 = ispinin
        endif
       
        if( Basereel ) then
          if( Cal_comp ) then
            Ampl_dft(Lm01,is1,Lm02,is2,mpirank0) = sum( norm_t(1:nlmso_reg,Lmf) &
                                                 * cmplx( smr(1:nlmso_reg,ii), smi(1:nlmso_reg,ii), db ))
          else
            Ampl_dft(Lm01,is1,Lm02,is2,mpirank0) = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,ii) )
          endif
        else
          Ampl_dft(Lm01,is1,Lm02,is2,mpirank0) = cmplx(smr(Lmf,ii),smi(Lmf,ii),db)
        endif

        if( Repres_comp .and. .not. Spinorbite ) &
          Ampl_dft(Lm01-2*m1,is1,Lm02-2*m2,is2,mpirank0) = (-1)**(L1+L2) * Ampl_dft(Lm01,is1,Lm02,is2,mpirank0)
!          Ampl_dft(Lm01-2*m1,is1,Lm02-2*m2,is2,mpirank0) = Ampl_dft(Lm01,is1,Lm02,is2,mpirank0)
!          Ampl_dft(Lm01-2*m1,is1,Lm02-2*m2,is2,mpirank0) = (-1)**(m1+m2) * Ampl_dft(Lm01,is1,Lm02,is2,mpirank0)

      end do
          
    end do
  endif
    
  if( nab_coop > 0 ) then
    iab = 0
    do ia = 1,natome

      ipra = iaprotoi(ia)

      do ib = ia+1,natome

        if( n_atom_coop > 0 ) then
          do j = 1,n_atom_coop
            if( ia_coop(j) == 0 ) cycle
            if( ia == ia_coop(j) .or. ib == ia_coop(j) ) exit
          end do
          if( j > n_atom_coop ) cycle
        endif

        iprb = iaprotoi(ib)
        Dist = sqrt( sum( ( Posi(:,ia) - Posi(:,ib) )**2 ) )
        if( Dist < Dist_coop(1) - eps10 ) cycle
        if( Dist_coop(2) > eps10 ) then
          if( Dist > Dist_coop(2) + eps10 ) cycle
        else
          if( Dist > Rmtg( ipra ) + Rmtg( iprb ) ) cycle
        endif

        iab = iab + 1
          
        do Lm1 = 1,nlmsa(ia)
          L1 = Lato(Lm1,ia,igrph)
          m1 = mato(Lm1,ia,igrph)
          Lm01 = L1**2 + L1 + 1 + m1
          if( Spinorbite ) then
            is1 = iato(Lm1,ia,igrph)
          else
            is1 = ispinin
          endif
          ii = ianew(ia) + Lm1

          do Lm2 = 1,nlmsa(ib)

            L2 = Lato(Lm2,ib,igrph)
            m2 = mato(Lm2,ib,igrph)
            Lm02 = L2**2 + L2 + 1 + m2
            if( nspino == 2 ) then
              is2 = iato(Lm2,ib,igrph)
            else
              is2 = ispinin
            endif
            jj = ianew(ib) + Lm2

            cfac = ( 0._db, 0._db )

            do Lmf = 1,nlmso_reg

              if( Basereel ) then
                if( Cal_comp ) then
                  ampl1 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,ii), smi(1:nlmso_reg,ii), db ))
                  ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,jj), smi(1:nlmso_reg,jj), db ))
                else
                  ampl1 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,ii) )
                  ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,jj) )
                endif
              else
                ampl1 = cmplx(smr(Lmf,ii),smi(Lmf,ii),db)
                ampl2 = cmplx(smr(Lmf,jj),smi(Lmf,jj),db)
              endif

              cfac = cfac + ampl1 * conjg( ampl2 )
            end do

! One multiplies by -img in order -aimag(Taull) is the density of state
! Conjugate to get the real part same sign than Green
            cfac = - img * cfac

            Tau_coop(Lm01,is1,Lm02,is2,iab,1) = Tau_coop(Lm01,is1,Lm02,is2,iab,1) + cfac
            if( Spinorbite ) Tau_coop(Lm02,is2,Lm01,is1,iab,2) = Tau_coop(Lm02,is2,Lm01,is1,iab,2) - conjg( cfac )

            if( Repres_comp .and. .not. Spinorbite ) then
              Lm01c = Lm01 - 2 * m1
              Lm02c = Lm02 - 2 * m2
              isg = (-1)**(m1+m2)
              irep = abs( irep_util(igrph,1) )
              if( abs(karact(49,irep)) > eps10 .or. abs(karact(53,irep)) > eps10 ) then   ! 3 fold axis case
                Tau_coop(Lm01c,is1,Lm02c,is2,iab,1) = Tau_coop(Lm01c,is1,Lm02c,is2,iab,1) + isg * cfac
              else  
! Empirical finding of the power of img in order sym 2/m with Ycomp and 4/m give same result
                Tau_coop(Lm01c,is1,Lm02c,is2,iab,1) = Tau_coop(Lm01c,is1,Lm02c,is2,iab,1) + isg * img**(-m2+1) * cfac
              endif
            endif
          
          end do          
        end do

      end do
    end do

  endif
  
  deallocate( norm, norm_t )
  deallocate( smi, smr )
  deallocate( lb1, lb2 )
  deallocate( newinv )

! Modification to be equivalent by symmetry with conventionnal 0.5(Tau - Tau^t*) given by multiple scattering theory
! With real harmonics, real part is zero
! with complex harmonics and no spin-orbit Tau(Lm,L'm') = -1^(m+m') Tau(L'-m',L-m)
! With spinorbit and no magnetism Tau(L,m,s,L',m',s') = -1^(m+s-m'-s') Tau(L',-m',-s',L,-m,-s)
  do ia = 1,natome

    if( Spinorbite .and. nspin == 2 ) exit
    
    if( ia /= iaabsi .and. .not. State_all_r .and. Cal_xanes ) cycle

    do Lm1 = 1,nlmsa(ia)
      L1 = Lato(Lm1,ia,igrph)
      m1 = mato(Lm1,ia,igrph)
      Lm01 = L1**2 + L1 + 1 + m1
      if( Spinorbite ) then
        is1 = iato(Lm1,ia,igrph)
      else
        is1 = ispinin
      endif

      do Lm2 = 1,nlmsa(ia)

        jj = ianew(ia) + Lm2

        L2 = Lato(Lm2,ia,igrph)
        m2 = mato(Lm2,ia,igrph)
        Lm02 = L2**2 + L2 + 1 + m2
        if( nspino == 2 ) then
          is2 = iato(Lm2,ia,igrph)
        else
          is2 = ispinin
        endif

        if( Spinorbite ) then
          isg = (-1)**(m1-m2-is1+is2)
          cfac = 0.5_db * ( Taull(Lm01,is1,Lm02,is2,ia) + isg * Taull(Lm02-2*m2,3-is2,Lm01-2*m1,3-is1,ia) )
          Taull(Lm01,is1,Lm02,is2,ia) = cfac
          Taull(Lm02-2*m2,3-is2,Lm01-2*m1,3-is1,ia) = isg * cfac
        elseif( Ylm_comp ) then
          isg = (-1)**(m1+m2)
          cfac = 0.5_db * ( Taull(Lm01,is1,Lm02,is2,ia) + isg * Taull(Lm02-2*m2,is2,Lm01-2*m1,is1,ia) )
          Taull(Lm01,is1,Lm02,is2,ia) = cfac
          Taull(Lm02-2*m2,is2,Lm01-2*m1,is1,ia) = isg * cfac
        else
          Taull(Lm01,is1,Lm02,is2,ia) = cmplx( 0._db, aimag(Taull(Lm01,is1,Lm02,is2,ia)), db )
        endif

      end do
    end do
  end do

! When State_all, recop is false.
  if( Recop .and. igrph == ngrph .and. ( Spinorbite .or. ispinin == nspin ) ) then
    allocate( taull_tem(nlmagm,nspinp,nlmagm,nspinp) )
    do ia = 1,natome
      if( ia /= iaabsi .and. .not. Atom_axe(ia) ) Cycle
      taull_tem(:,:,:,:) = Taull(:,:,:,:,ia)
      call recop_taull(nlmagm,nspinp,Sym_cubic,Taull_tem,Ylm_comp)
      Taull(:,:,:,:,ia) = taull_tem(:,:,:,:)
    end do
    deallocate( taull_tem )
    if( RIXS_gen .and. Cal_xanes ) then
  ! not used  (because recop set to false when RIXS_gen in fdm.90)
      allocate( taull_tem(nlm_rixs,nspinp,nlmomax,nspinp) )
      taull_tem(:,:,:,:) = Ampl_dft(:,:,:,:,mpirank0)
      call recop_ampl(nlm_rixs,nlmomax,nspinp,Sym_cubic,Taull_tem,Ylm_comp)
      Ampl_dft(:,:,:,:,mpirank0) = taull_tem(:,:,:,:) 
      deallocate( taull_tem )
    endif
  endif

  if( icheck > 1 .and. igrph == ngrph .and. ( Spinorbite .or. ispinin == nspin ) ) then
    do ia = 1,natome
      if( ia /= iaabsi .and. .not. State_all_r .and. Cal_xanes ) cycle
      if( ia == iaabsi ) then
        write(3,210)
      else
        write(3,220) ia
      endif
      Lmw = min(3,Lmaxa(ia))
      nlmw = min(nlmagm,(Lmw+1)**2)
      if( nspinp == 2 ) then
        write(3,230) ((( L, m, i, m = -L,L), L = 0,Lmw ), i = 1,2)
      else
        write(3,240) (( L, m, m = -L,L), L = 0,Lmw )
      endif
      do is1 = 1,nspinp
        Lm1 = 0
        do L1 = 0,Lmw
          do m1 = -L1,L1
            Lm1 = Lm1 + 1
            write(3,250) L1, m1, is1, (( Taull(Lm1,is1,Lm2,is2,ia), Lm2 = 1,nlmw), is2 = 1,nspinp )
          end do
        end do
      end do
    end do
  endif

  return
  105 format(/' nlmso_eg =',i4,' > nlmso_reg =',i4,// '   Increase Lmaxso !') 
  110 format(/' ---- Mat_FDM -------',100('-'),//' igrph =',i2)
  120 format(/'   ii   Solution(i,Lmf = 1,nlmso_reg)')
  130 format(i6,1p,250(1x,2e11.3))
  140 format(i6,1p,500e11.3)
  210 format(/' Multiple scattering amplitude')
  220 format(/' Atom ia =',i3,', multiple scattering amplitude')
  230 format('( l, m, s)',18(8x,3i3,6x))
  240 format('( l, m, s)',18(11x,2i3,7x))
  250 format(3i3,2x,1p,18(1x,2e11.3))

end

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

! Calculation of the Neuman and Bessel functions in the FDM grid of point at the border of the outer-sphere

subroutine Phiso(Adimp,Bessel,Besselr,E_comp,Ecinetic_out,Eclie_out,Eimag, &
       Eneg,Enervide,icheck,isp,isrt,lmaxso,mpirank0,Neuman,Neumanr,npsom,nsort,nsort_c,nsort_r,nspinr, &
       R_Rydb,Rsort,Rydberg,V0bd_out,xyz)

  use declarations
  implicit none

  integer:: i, ib, icheck, ipr, ir, isp, j, L, L2, lmaxso, mpirank0, npsom, nr, nsort, nsort_c, nsort_r, nspinr

  integer, dimension(nsort):: isrt

  complex(kind=db):: bs(0:lmaxso), fnormc, konde, nm(0:lmaxso), z
  complex(kind=db), dimension(nsort_c,0:lmaxso,nspinr):: Bessel, Neuman
  complex(kind=db), dimension(:), allocatable:: f12
  complex(kind=db), dimension(:,:), allocatable:: rbs, rnm, u

  logical:: E_comp, Eneg, Rydberg

  real(kind=db):: Adimp, cal_norm, clapl, clapl0, deltar, dr, Ecinetic_out, Ec, Eclie_out, ee, E_min, Eimag, Enervide, fnorm, &
                  konder, p1, p2, pp, R_rydb, Rmax, rr, Rsort, V0bd_out, zr

  real(kind=db), dimension(3):: p
  real(kind=db), dimension(0:lmaxso):: bsr, nmr
  real(kind=db), dimension(4,npsom):: xyz
  real(kind=db), dimension(nsort_r,0:lmaxso,nspinr):: Besselr, Neumanr
  real(kind=db), dimension(:), allocatable:: f1, f2, f12r, r, v
  real(kind=db), dimension(:,:), allocatable:: rbsr, rnmr, ur

  if( icheck > 1 ) write(3,110) Ecinetic_out * rydb, lmaxso

  if( E_comp ) then
    konde = sqrt( cmplx( Ecinetic_out, Eimag, db ) )
    if( abs( konde ) < eps10 ) konde = cmplx( eps10, 0._db, db )
    konder = real( konde, db )
!    fnorm = cal_norm(Ecinetic_out,eimag)
!    fnormc = cmplx( fnorm, 0._db, db)
    fnormc = sqrt( konde / pi )
    
    fnorm = real( fnormc, db )
  else
    konder = sqrt( Ecinetic_out )
    konde = cmplx( konder, 0._db, db )
    if( konder < eps10 .and. mpirank0 == 0 ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,120)
      end do
      stop
    endif
    fnorm = sqrt( konder / pi )
    fnormc = sqrt( konde / pi )
  endif

  if( Rydberg ) then
  
    E_min = 1.e-06_db / Rydb   
    if( abs( Enervide ) < E_min ) then
      ee = E_min
    else
      ee = Enervide
    endif
    if( ee < E_min ) then
      Rmax = - 2 / ee + R_rydb
    else
      Rmax = 100._db / bohr
    endif
  
  else
  
    Rmax = 0._db
  
  endif

  if( Rydberg .and. Rmax >= Rsort ) then
    deltar = Rmax - Rsort + adimp + eps10
    dr = 0.05_db / bohr
    nr = nint( deltar / dr ) + 1

    allocate( r(0:nr) )
    allocate( v(0:nr) )
    allocate( f1(0:nr) )
    allocate( f2(0:nr) )

    dr = deltar / ( nr - 1 )
    do ir = nr,0,-1
      r(ir) = rmax - ( ir - 1 ) * dr
    end do
    v(0:nr) = - 2 / r(0:nr)
    do ir = nr,0,-1
      v(ir) = max( v(ir), V0bd_out )
    end do

! To avoid discontinuity
    if( ee >= E_min ) then
      pp = - v(0) / ( r(0) - r(nr) )
      do ir = 0,nr
        v(ir) = v(ir) + ( r(ir) - r(nr) ) * pp
      end do
    endif

    if( E_comp ) then
      allocate( f12(0:nr) )
      allocate( u(-1:nr,2) )
      allocate( rnm(2,0:lmaxso) )
      allocate( rbs(2,0:lmaxso) )
    else
      allocate( f12r(0:nr) )
      allocate( ur(-1:nr,2) )
      allocate( rnmr(2,0:lmaxso) )
      allocate( rbsr(2,0:lmaxso) )
    endif

    Ec = Enervide - v(0)
    if( .not. Eneg ) Ec = max( Ec, Eclie_out )
    if( E_comp ) then
      konde = sqrt( cmplx( Ec, eimag, db ) )
      konder = real( konde, db )
    else
      konder = sqrt( ec )
      konde = cmplx( konder, 0._db, db )
    endif
    if( konder < eps10 .and. mpirank0 == 0 ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,120)
      end do
      stop
    endif
    if( E_comp ) then
      fnorm = cal_norm(Ecinetic_out,eimag)
      fnormc = cmplx( fnorm, 0._db, db)
    else
      fnorm = sqrt( konder / pi )
      fnormc = sqrt( konde / pi )
    endif

    do i = 1,2
      rr = rmax + i * dr
      if( E_comp ) then
        z = konde * rr
        call cbessneu(fnormc,z,lmaxso,lmaxso,bs,nm)
        rnm(i,0:lmaxso) = nm(0:lmaxso) * rr
        rbs(i,0:lmaxso) = bs(0:lmaxso) * rr
      else
        zr = konder * rr
        call cbessneur(fnorm,zr,lmaxso,lmaxso,bsr,nmr)
        rnmr(i,0:lmaxso) = nmr(0:lmaxso) * rr
        rbsr(i,0:lmaxso) = bsr(0:lmaxso) * rr
      endif
    end do

    clapl = 1 / dr**2
    clapl0 = - 2 * clapl
    f2(0:nr) = 1 / r(0:nr)**2
    f1(0:nr) = - clapl0 + v(0:nr) - Enervide

    do L = 0,lmaxso
      L2 = L**2 + L

      if( E_comp ) then
        f12(0:nr) = cmplx( f1(0:nr) + L2 * f2(0:nr), -eimag, db ) / clapl
        u(0,1) = rnm(1,L)
        u(-1,1) = rnm(2,L)
        u(0,2) = rbs(1,L)
        u(-1,2) = rbs(2,L)
        do j = 1,2
          do ir = 0,nr-1
            u(ir+1,j) = f12(ir) * u(ir,j) - u(ir-1,j)
          end do
          u(0:nr,j) = u(0:nr,j) / r(0:nr)
        end do
      else
        f12r(0:nr) = ( f1(0:nr) + L2 * f2(0:nr) ) / clapl
        ur(0,1) = rnmr(1,L)
        ur(-1,1) = rnmr(2,L)
        ur(0,2) = rbsr(1,L)
        ur(-1,2) = rbsr(2,L)
        do j = 1,2
          do ir = 0,nr-1
            ur(ir+1,j) = f12r(ir) * ur(ir,j) - ur(ir-1,j)
          end do
          ur(0:nr,j) = ur(0:nr,j) / r(0:nr)
        end do
      endif

      do ib = 1,nsort
        i = isrt(ib)
        p(1:3) = xyz(1:3,i)
        rr = sqrt( sum( p(:)**2 ) )
        do ir = nr-1,2,-1
          if( r(ir) > rr ) exit
        end do
        p1 = ( rr - r(ir) ) / ( r(ir+1) - r(ir) )
        p2 = 1 - p1

        if( E_comp ) then
          Neuman(ib,L,isp) = p1 * u(ir+1,1) + p2 * u(ir,1)
          Bessel(ib,L,isp) = p1 * u(ir+1,2) + p2 * u(ir,2)
        else
          Neumanr(ib,L,isp) = p1 * ur(ir+1,1) + p2 * ur(ir,1)
          Besselr(ib,L,isp) = p1 * ur(ir+1,2) + p2 * ur(ir,2)
        endif

      end do

      if( icheck > 2 ) then
        write(3,126) L
        do ir = nr,0,-1
          if( E_comp ) then
            write(3,127) r(ir)*bohr, v(ir)*rydb, u(ir,1:2)
          else
            write(3,127) r(ir)*bohr, v(ir)*rydb, ur(ir,1:2)
          endif
        end do
      endif

    end do

    deallocate( r )
    deallocate( v )
    deallocate( f1 )
    deallocate( f2 )
    if( E_comp ) then
      deallocate( f12 )
      deallocate( u )
      deallocate( rnm )
      deallocate( rbs )
    else
      deallocate( f12r )
      deallocate( ur )
      deallocate( rnmr )
      deallocate( rbsr )
    endif

  else

! Calculation of Bessel and Neuman function
    do ib = 1,nsort
      i = isrt(ib)
      p(1:3) = xyz(1:3,i)
      if( E_comp ) then
        z = konde * sqrt( sum( p(:)**2 ) )
        call cbessneu(fnormc,z,lmaxso,lmaxso,bs,nm)
        Neuman(ib,0:lmaxso,isp) = nm(0:lmaxso)
        Bessel(ib,0:lmaxso,isp) = bs(0:lmaxso)
      else
        zr = konder * sqrt( sum( p(:)**2 ) )
        call cbessneur(fnorm,zr,lmaxso,lmaxso,bsr,nmr)
        Neumanr(ib,0:lmaxso,isp) = nmr(0:lmaxso)
        Besselr(ib,0:lmaxso,isp) = bsr(0:lmaxso)
      endif
    end do

  endif

  if( icheck > 2 ) then
    do L = 0,lmaxso
      if( E_comp ) then
        write(3,130) L, isp, konde
        do ib = 1,nsort
          i = isrt(ib)
          p(1:3) = xyz(1:3,i)
          rr = sqrt( sum( p(:)**2 ) )
          z = konde * rr
          write(3,140) isrt(ib), rr*bohr, Bessel(ib,L,isp), Neuman(ib,L,isp)
        end do
      else
        write(3,150) L, isp, konder
        do ib = 1,nsort
          i = isrt(ib)
          p(1:3) = xyz(1:3,i)
          rr = sqrt( sum( p(:)**2 ) )
          zr = konder * rr
          write(3,140) isrt(ib), rr*bohr, Besselr(ib,L,isp), Neumanr(ib,L,isp)
        end do
      endif
    end do
  endif

  return
  110 format(/' ---- Phiso --------',100('-')//, 20x,' Ecinetic_out =',f10.5,' eV,  lmaxso =',i3)
  120 format(//' The wave vector is zero, what is forbidden !'/)
  126 format(/'    radius      v            bess-coul             neuman-coul     l =',i3)
  127 format(2f10.3,1p,4e12.4)
  130 format(/' isrt     r',13x,'bessel',18x,'neuman',9x, 'l =',i3,', isp =',i2,', konde =',1p,2e12.4)
  140 format(i5,f8.3,1x,1p,4e12.4)
  150 format(/' isrt     r',7x,'bessel      neuman',5x, 'l =',i3,', isp =',i2,', konder =',1p,e12.4)
end

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

! Calculate one row abvr, abvi and smr, smi

subroutine CalcMatRow( abvr, abvi, Base_hexa, Basereel, Bessel, Besselr, Cal_comp, cgrad, clapl, E_comp, Eimag, &
    Enervide, gradvr, ianew, ianewL, iato, ibord, icheck, igrph, ii, isbord, iso, ispin0, isrt, &
    isvois, ivois, Kar, Kari, Lato, Lb1i, Lb1r, Lb2i, Lb2r, lmaxso, Lso, mato, mletl, mso, &
    natome, nbm, nbm_c, nbord, nbordf, nbtm, Neuman, Neumanr, new, newinv, ngrph, nicm, nim, nligne, &
    nligne_i, nligneso, nlmagm, nlmmax, nlmomax, nlmsa, nlmsam, nlmso, nlmso_i, nphiato, nphiato_c, npoint, npsom, &
    nso1, nsort, nsort_c, nsort_r, nsortf, nspin, nspino, nspinp, nspinr, &
    numia, nvois, Phiato, Phiato_c, poidsa, poidso, Relativiste, Repres_comp, Rot_Spinor, rvol, Spin_axis_free, Spinorbite, &
    smi, smr, Vr, Ylm_comp, Ylmato, Ylms, Ylmso, MUMPS )

  use declarations
  implicit none
  
  integer, parameter:: nletm = 52

  integer:: i, ii1, ia, ib, icheck, ii, igrph, index_Rotational, indg, isoo, isop, isp, ispin, ispin_2, ispin0, ispint, &
    ispp, ispq, ispt, iss, isss, isym, iv, j, jj, jjj, k, &
    L, Lb1i, Lb1r, Lb2i, Lb2r, Ljj, Lm, Lm0, lmaxso, Lmf, lmf0, Lmp, Lmp0, Lms, Lmsp, Ln, lp, m, mf, mp, &
    n, natome, nbm, nbm_c, nbtm, ngrph, nicm, nim, nligne, nligne_i, nligneso, & 
    nlmagm, nlmso, nlmso_i, nphiato, nphiato_c, npoint, npsom, nsort, nsort_c, nsortf, nsort_r, nspino, nspinp, nso1, &
    nlmsam, nspin, nspinr, nlmmax, nlmomax, nvois
 
  integer, dimension(0:npoint):: new
  integer, dimension(natome):: ianew, ianewL, nbord, nbordf, nlmsa
  integer, dimension(nsort):: isrt
  integer, dimension(npsom):: numia
  integer, dimension(nso1,ngrph):: Lso, mso, iso
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato
  integer, dimension(nbtm,natome):: ibord, isbord
  integer, dimension(npsom,nvois):: ivois, isvois
  integer, dimension(nligne):: newinv
  integer, dimension(:), allocatable:: jnzero

  integer, save:: letd, letu, letc, nlet

  character(len=1):: let(nletm)
  character(len=3):: mlet(nletm**3)
  character(len=3), dimension(nligne):: mletl

  complex(kind=db):: cfac, cfad, charmc, Phi_out, Yc, Ycomp
  complex(kind=db), dimension(2,2):: Mat_Pauli_Rot
  complex(kind=db), dimension(nspino):: cfv_comp, cfvs_c, cso
  complex(kind=db), dimension(nsort_c,0:lmaxso,nspino):: Bessel, Neuman
  complex(kind=db), dimension(nphiato_c,nlmagm,nspinp,nspino,natome):: Phiato_c
  complex(kind=db), dimension(nbm_c,nlmmax,nspinp,nspino,natome):: Ylms

  logical:: Base_hexa, Basereel, Cal_comp, E_comp, MUMPS, Relativiste, Repres_comp, Spin_axis_free, Spinorbite, &
            Ylm_comp

  real(kind=db):: a2s4, bder, charmr, charmi, crelat, Enervide, Eimag, f, fac, fad, fr, VmE
  
  real(kind=db), dimension(2):: Charac_i, Charac_r, csor, csoi, cfvs_r, cfvs_i, Fac_Rotational, v
  real(kind=db), dimension(3,3):: Rot_Spinor
  real(kind=db), dimension(nopsm,nspino):: Kar, Kari
  real(kind=db), dimension(nvois):: cgrad
  real(kind=db), dimension(0:nvois):: clapl
  real(kind=db), dimension(npoint,nspin):: Vr
  real(kind=db), dimension(nbm,natome):: poidsa
  real(kind=db), dimension(nim):: rvol
  real(kind=db), dimension(nsortf):: poidso
  real(kind=db), dimension(nsort,nlmomax):: Ylmso
  real(kind=db), dimension(nicm,3,nspin):: gradvr
  real(kind=db), dimension(nbtm,nlmmax,natome):: Ylmato
  real(kind=db), dimension(nphiato,nlmagm,nspinp,nspino,natome):: Phiato
  real(kind=db), dimension(Lb1i:Lb2i):: abvi
  real(kind=db), dimension(Lb1r:Lb2r):: abvr
  real(kind=db), dimension(nlmso,nligne):: smr
  real(kind=db), dimension(nlmso_i,nligne_i):: smi
  real(kind=db), dimension(nsort_r,0:lmaxso,nspino):: Besselr, Neumanr
  
  real(kind=db), save, dimension(nletm**3):: vletr, vleti

  data let/ 'a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z', &
            'A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/

  ispin = ispin0
  if( icheck > 2 .and. ( ( ii == nligne .and. .not. MUMPS ) .or. ( ii == 1 .and. MUMPS ) ) ) then
    write(3,120)
    nlet = 0
    letu = 0
    letd = 1
    letc = 1
  endif
  
! alfa_sf is the fine structure constant
  a2s4 = 0.25_db * alfa_sf**2

! ii is the index of line in the FDM matrix

  i = newinv(ii)
! i > 0 is the index in the FDM grid of point
! i = 0 means in the outer sphere
! i < 0 means abs(i) = ia is the index of the atom where is the point

  if( Spinorbite .and. i > 0 ) then
    do ia = natome,1,-1
      if( ii > ianewL(ia) ) exit
    end do
    if( ia == 0 ) then
      ispin = 2 - mod(ii,2)
    else
      ii1 = ii - ( ianewL(ia) + nlmsa(ia) )
      ispin = 2 - mod(ii1,2)
    endif
  endif 

! Points of the FDM grid
  if( i > 0 ) then

    ispp = min( ispin, nspin )
    indg = 3 - ispin
    ispq = min( indg, nspin )

    isp = min( ispin, nspino )

! Diagonal term:
    VmE = Vr( i, ispp ) - Enervide
    abvr(ii) = VmE - clapl(0)
    if( E_comp ) abvi(ii) = - Eimag

    if( Relativiste ) abvr(ii) = abvr(ii) - a2s4 * VmE**2
    if( Relativiste .or. Spinorbite ) bder = - a2s4 / ( 1 - a2s4 * VmE )

! iv is the number of the neighboring points
    do iv = 1,nvois

! j is the index of the neighboring point
      j = ivois(i,iv)
      ia = numia(j)

! Coefficient connecting the central point to the neighboring point
      cfvs_r(isp) = - clapl(iv)
      if( Cal_comp ) cfvs_i(isp) = 0._db 

      Charac_r(1:nspino) = Kar( abs( isvois(i,iv) ), 1:nspino )
      Charac_i(1:nspino) = Kari( abs( isvois(i,iv) ), 1:nspino )

      if( Spinorbite .or. Relativiste ) then
        if( Base_hexa ) then
          k = 1 + mod((iv-1)/2,4)
        else
          k = 1 + mod((iv-1)/2,3)
        endif
      endif

      if( Relativiste ) then
        if( k == 4 ) then
          crelat = bder * cgrad(iv) * gradvr(i,2,ispp)
        else
          crelat = bder * cgrad(iv) * gradvr(i,k,ispp)
        endif
        cfvs_r(isp) = cfvs_r(isp) + crelat
      endif

      if( Spinorbite ) then
! rotationnal   iS.( gradV x grad ), S = Pauli matrices
        cso(:) = ( 0._db, 0._db )
        
        do index_Rotational = 1,3
! There is a rotation of the Pauli matrices vector when the spin axis is not parallel to the z axis of the FDM grid of points
! When parallel, Mat_Pauli_Rot = Pauli-matrix(index_Rotational) 
          Mat_Pauli_Rot(1,1) = cmplx( Rot_Spinor(index_Rotational,3), 0._db, db )
          Mat_Pauli_Rot(1,2) = cmplx( Rot_Spinor(index_Rotational,1), -1._db * Rot_Spinor(index_Rotational,2), db ) 
          Mat_Pauli_Rot(2,1) = conjg( Mat_Pauli_Rot(1,2) )
          Mat_Pauli_Rot(2,2) = - Mat_Pauli_Rot(1,1)
                
          Mat_Pauli_Rot(:,:) = img * Mat_Pauli_Rot(:,:)
!if( iv == 1 ) then
!  write(6,'(a18,i3)') 'index_Rotational =', index_Rotational          
!  write(6,'(a10,1p,2(1x,2e13.5))') 'Mat_Pauli=', Mat_Pauli_Rot(1,:)          
!  write(6,'(10x,1p,2(1x,2e13.5))') Mat_Pauli_Rot(2,:)
!endif
          Fac_Rotational(:) = 0._db
          
          do iss = 1,nspino
            isss = min( iss, nspin )
                      
            select case(index_Rotational)
            
              case(1)
              
                if( k == 2 .or. k == 4 ) then
                  Fac_Rotational(iss) = - GradVr(i,3,isss)
                elseif( k == 3 ) then
                  Fac_Rotational(iss) = GradVr(i,2,isss)
                endif
                
              case(2)
              
                if( k == 1 ) then
                  Fac_Rotational(iss) = GradVr(i,3,isss)
                elseif( k == 3 ) then
                  Fac_Rotational(iss) = - GradVr(i,1,isss)
                endif
                   
              case(3)
              
                if( k == 1 ) then
                  Fac_Rotational(iss) = - GradVr(i,2,isss)
                elseif( k == 2 ) then
                  Fac_Rotational(iss) = GradVr(i,1,isss)
                endif
                   
            end select
          
          end do

          cso(isp) = cso(isp) + Mat_Pauli_Rot(isp,isp) * Fac_Rotational(isp)   
          cso(indg) = cso(indg) + Mat_Pauli_Rot(isp,indg) * Fac_Rotational(indg)   

        end do

        csor(:) = real( cso(:), db )
        csoi(:) = aimag( cso(:) )

! Below is equivalent to above when no rotation of Spinor
!        csor(:) = 0._db
!        csoi(:) = 0._db
!        select case(k)
!          case(1)
!            if( isp == 1 ) then
!              csor(indg) = gradvr(i,3,ispq)   ! Sy
!              csoi(isp) = - gradvr(i,2,ispp)  ! Sz
!            else
!              csor(indg) = - gradvr(i,3,ispq) ! Sy
!              csoi(isp) = gradvr(i,2,ispp)    ! Sz
!            endif
!          case(2,4)
!            csoi(indg) = - gradvr(i,3,ispq)   ! Sx
!            if( isp == 1 ) then
!              csoi(isp) = gradvr(i,1,ispp)    ! Sz
!            else
!              csoi(isp) = - gradvr(i,1,ispp)  ! Sz
!            endif  
!          case(3)
!            csoi(indg) = gradvr(i,2,ispq)     ! Sy
!            if( isp == 1 ) then
!              csor(indg) = - gradvr(i,1,ispq) ! Sz 
!            else
!              csor(indg) = gradvr(i,1,ispq)   ! Sz
!            endif
!        end select

! Below is not equivalent to above. It was in the code up to June 2025
! Sign ok only when isp = 1 and case(3) would have to contain csoi. See just above to see the good equations.
!        csor(:) = 0._db
!        csoi(:) = 0._db
!        select case(k)
!          case(1)
!            csoi(isp) = - gradvr(i,2,ispp)
!            csor(indg) = gradvr(i,3,ispq)
!          case(2,4)
!            csoi(isp) = gradvr(i,1,ispp)
!            csoi(indg) = - gradvr(i,3,ispq)
!          case(3)
!            csor(indg) = - gradvr(i,1,ispq) + gradvr(i,2,ispq)
!        end select

        csor(:) = bder * cgrad(iv) * csor(:)
        csoi(:) = bder * cgrad(iv) * csoi(:)

        cfvs_r(isp) = cfvs_r(isp) + csor(isp)
        cfvs_r(indg) = csor(indg)
        cfvs_i(:) = csoi(:)
        if( Repres_comp ) then
          v(:) = Charac_r(:) * cfvs_r(:) - Charac_i(:) * cfvs_i(:)
          cfvs_i(:) = Charac_i(:) * cfvs_r(:) + Charac_r(:) * cfvs_i(:)
          cfvs_r(:) = v(:)
        else
          cfvs_i(:) = Charac_r(:) * cfvs_i(:)
          cfvs_r(:) = Charac_r(:) * cfvs_r(:)
        endif
      else
        if( Repres_comp ) cfvs_i(isp) = Charac_i(isp) * cfvs_r(isp)
        cfvs_r(isp) = Charac_r(isp) * cfvs_r(isp)
      endif

      cfvs_r(1:nspino) = cfvs_r(1:nspino) * rvol(i)
      if( Cal_comp ) then
        cfvs_i(1:nspino) = cfvs_i(1:nspino) * rvol(i)
        cfvs_c(1:nspino) = cmplx( cfvs_r(1:nspino), cfvs_i(1:nspino), db )
      endif

! Neighbor is in the interstitial FDM grid :
      if( ia == 0 ) then

        if( Spinorbite ) then
          jj = new(j) - nspino
          do ispin_2 = 1,nspino
            jjj = jj + ispin_2
            abvr( jjj ) = abvr( jjj ) + cfvs_r(ispin_2) / rvol(j)
            if( Cal_comp ) abvi( jjj ) = abvi( jjj ) + cfvs_i(ispin_2) / rvol(j)
          end do
        else
          jjj = new(j)
          abvr( jjj ) = abvr( jjj ) + cfvs_r(1) / rvol(j)
          if( Cal_comp ) abvi( jjj ) = abvi( jjj ) + cfvs_i(1) / rvol(j)
        endif

      elseif( ia > 0 ) then

! Neighbor is in the atom ia :

        do ib = nbordf(ia)+1,nbord(ia)
          if( j == ibord(ib,ia) ) exit
        end do
        isym = abs( isbord(ib,ia) )
        if( isym /= 1 ) then
          do ispin_2 = 1,nspino
            if( Repres_comp ) then
              cfvs_c(ispin_2) = cmplx( Kar(isym,ispin_2), Kari(isym,ispin_2), db ) * cfvs_c(ispin_2)
              cfvs_r(ispin_2) = real( cfvs_c(ispin_2), db)
              cfvs_i(ispin_2) = aimag( cfvs_c(ispin_2) )
            elseif( Cal_comp ) then
              cfvs_c(ispin_2) = Kar(isym,ispin_2) * cfvs_c(ispin_2)
              cfvs_r(ispin_2) = real( cfvs_c(ispin_2), db)
              cfvs_i(ispin_2) = aimag( cfvs_c(ispin_2) )
            else
              cfvs_r(ispin_2) = Kar(isym,ispin_2) * cfvs_r(ispin_2)
            endif
          end do
        endif
        
        do Lm = 1,nlmsa(ia)

          jjj = ianew(ia) + Lm
          L = Lato(Lm,ia,igrph)
          Lms = L**2 + L + 1 + mato(Lm,ia,igrph) 
          isoo = iato(Lm,ia,igrph)

          do ispin_2 = 1,nspino       ! it is the spin. When ispin_2 /= ispin, this is the spin-flip process
            if( Spin_axis_free ) then
              m = mato(Lm,ia,igrph)
            else
              m = mato(Lm,ia,igrph) + ispin_2 - isoo  ! the 'real' m
              if( m > L .or. m < -L ) cycle
            endif
            Lm0 = L**2 + L + 1 + m
            
            if( Cal_comp ) then
 
              if( nphiato_c > 0 .and. Spinorbite ) then
                cfac = cfvs_c(ispin_2) * Phiato_c(ib,Lms,ispin_2,isoo,ia)
              elseif( nphiato_c > 0 ) then 
                cfac = cfvs_c(isp) * Phiato_c(ib,Lms,ispin,1,ia)
              elseif( Spinorbite ) then
                cfac = cfvs_c(ispin_2) * Phiato(ib,Lms,ispin_2,isoo,ia)
              else
                cfac = cfvs_c(isp) * Phiato(ib,Lms,ispin,1,ia)
              endif
            
              abvr( jjj ) = abvr( jjj ) + Real( cfac, db )
              abvi( jjj ) = abvi( jjj ) + aimag( cfac )
              
            else
            
              abvr( jjj ) = abvr( jjj ) + cfvs_r(isp) * Phiato(ib,Lms,ispp,isoo,ia)
 
            endif
 
          end do
      
        end do

! Neighbor is in the outer-sphere :
      elseif( ia == -2 ) then
!if( ii == 115) stop
        do ib = nsortf+1,nsort
          if( j == isrt(ib) ) exit
        end do
     
        do Lm = 1,nlmso
          L = Lso(Lm,igrph)
          m = mso(Lm,igrph)
          Lm0 = L**2 + L + 1 + m
          ispq = min( iso(Lm,igrph), nspinr )
          jjj = nligneso + Lm
      
          if( Cal_comp ) then   
            if( Ylm_comp .and. m /= 0 ) then
              Ycomp = Yc(m,Ylmso(ib,Lm0),Ylmso(ib,Lm0-2*m))
            else
              Ycomp = cmplx( Ylmso(ib,Lm0), 0._db, db )
            endif
            if( Basereel ) then
              if( E_comp ) then
                Phi_out = Neuman(ib,L,ispq)
              else
                Phi_out = cmplx( Neumanr(ib,L,ispq), 0._db, db ) 
              endif
            else
              if( E_comp ) then
                Phi_out = Bessel(ib,L,ispq) + img * Neuman(ib,L,ispq)
              else
                Phi_out = cmplx( Besselr(ib,L,ispq), Neumanr(ib,L,ispq), db )
              endif
            endif
            cfv_comp(1:nspino) = Phi_out * Ycomp * cmplx( cfvs_r(1:nspino), cfvs_i(1:nspino), db )
            if( isp == iso(Lm,igrph) ) then
              abvr( jjj ) = abvr( jjj ) + real( cfv_comp(isp), db)
              abvi( jjj ) = abvi( jjj ) + aimag( cfv_comp(isp) )
            else
      ! here only if spinorbit
              abvr( jjj ) = abvr( jjj ) + real( cfv_comp(indg),db)
              abvi( jjj ) = abvi( jjj ) + aimag( cfv_comp(indg) )
            endif
          else
            fr = Ylmso(ib,Lm0) * Neumanr(ib,L,ispq)
            abvr( jjj ) = abvr( jjj ) + cfvs_r(isp) * fr
          endif

        end do

        do Lm = 1,nlmso
          L = Lso(Lm,igrph)
          m = mso(Lm,igrph)
          Lm0 = L**2 + L + 1 + m
          ispq = min( iso(Lm,igrph), nspinr )

          if( Cal_comp ) then
            if( Ylm_comp .and. m /= 0 ) then
              Ycomp = Yc(m,Ylmso(ib,Lm0),Ylmso(ib,Lm0-2*m))
            else
              Ycomp = cmplx( Ylmso(ib,Lm0), 0._db, db )
            endif
            if( E_comp ) then
              cfv_comp(1:nspino) = Ycomp * Bessel(ib,L,ispq) * cmplx( cfvs_r(1:nspino), cfvs_i(1:nspino), db )
            else
              cfv_comp(1:nspino) = Ycomp * Besselr(ib,L,ispq) * cmplx( cfvs_r(1:nspino), cfvs_i(1:nspino), db )
            endif
            if( isp == iso(Lm,igrph) ) then
              smr(Lm,ii) = smr(Lm,ii) - real( cfv_comp(isp), db )
              smi(Lm,ii) = smi(Lm,ii) - aimag( cfv_comp(isp) )
            else
              smr(Lm,ii) = smr(Lm,ii) - real( cfv_comp(indg), db)
              smi(Lm,ii) = smi(Lm,ii) - aimag( cfv_comp(indg) )
            endif
          else
            fr = Ylmso(ib,Lm0) * Besselr(ib,L,ispq)
            smr(Lm,ii) = smr(Lm,ii) - cfvs_r(isp) * fr
          endif
        end do

      endif

    end do ! end of the loop on the neighbor points

! Expansion in the outer sphere
  elseif( i == 0 ) then

    Lm = ii - nligneso
    m = mso(Lm,igrph)
    Lm0 = Lso(Lm,igrph)**2 + Lso(Lm,igrph) + 1 + m
    isp = iso(Lm,igrph)

! coefficients of scattering amplitude
    do Lmp = 1,nlmso 

      if( iso(Lmp,igrph) /= isp ) cycle
      L = Lso(Lmp,igrph)
      mp = mso(Lmp,igrph)
      Lmp0 = L**2 + L + 1 + mp
      ispq = min( iso(Lmp,igrph), nspinr )
      Ljj = nligneso + Lmp
      if( Cal_comp ) then
        charmc = ( 0._db, 0._db )
      else
        charmr = 0._db
      endif

      do ib = 1,nsortf
        if( Cal_comp ) then
          if( Ylm_comp .and.  m /= 0 .and. mp /= 0 ) then
            cfac = conjg( Yc(m,Ylmso(ib,Lm0),Ylmso(ib,Lm0-2*m)) ) * Yc(mp,Ylmso(ib,Lmp0), Ylmso(ib,Lmp0-2*mp) )
          elseif( Ylm_comp .and.  m /= 0 ) then
            cfac = conjg( Yc(m,Ylmso(ib,Lm0),Ylmso(ib,Lm0-2*m)) ) * Ylmso(ib,Lmp0)
          elseif( Ylm_comp .and.  mp /= 0 ) then
            cfac = Ylmso(ib,Lm0) * Yc(mp,Ylmso(ib,Lmp0), Ylmso(ib,Lmp0-2*mp) )
          else
            cfac = cmplx(Ylmso(ib,Lm0) * Ylmso(ib,Lmp0), 0._db, db)
          endif
          cfac = poidso(ib) * cfac
 
          if( Basereel ) then
            if( E_comp ) then
              charmc = charmc + cfac * Neuman(ib,L,ispq)
            else
              charmc = charmc + cfac * Neumanr(ib,L,ispq)
            endif
          else
            if( E_comp ) then
              charmc = charmc + cfac * ( Bessel(ib,L,ispq) + img * Neuman(ib,L,ispq) )
            else
              charmc = charmc + cfac * cmplx( Besselr(ib,L,ispq), Neumanr(ib,L,ispq), db )
            endif
          endif
        else
          fac = poidso(ib) * Ylmso(ib,Lm0) * Ylmso(ib,Lmp0)
          charmr = charmr + fac * Neumanr(ib,L,ispq)
        endif
      end do

      if( Cal_comp ) then
        abvr(Ljj) = abvr(Ljj) + real( charmc, db )
        abvi(Ljj) = abvi(Ljj) + aimag( charmc )
      else
        abvr(Ljj) = abvr(Ljj) + charmr
      endif
      
    end do

! Coefficients of FDM grid amplitudes and of outer-sphere (sm)
    do ib = 1,nsortf

      j = isrt(ib)
      jj = new(j) - nspino + isp

      if( Ylm_comp .and. m /= 0 ) then
        cfac = poidso(ib) * conjg( Yc(m,Ylmso(ib,Lm0),Ylmso(ib,Lm0-2*m)) )
        abvr( jj ) = abvr( jj ) - real( cfac, db ) / rvol(j)
        abvi( jj ) = abvi( jj ) - aimag( cfac ) / rvol(j)
      else
        fac = poidso(ib) * Ylmso(ib,Lm0)
        if( Cal_comp ) cfac = cmplx( fac, 0._db, db )
        abvr( jj ) = abvr( jj ) - fac / rvol(j)
      endif

      do Lmf = 1,nlmso
        ispq = min( iso(Lmf,igrph), nspinr )
        if( isp /= iso(Lmf,igrph) ) cycle
        L = Lso(Lmf,igrph)
        mf = mso(Lmf,igrph)
        lmf0 = L**2 + L + 1 + mf
        if( Cal_comp ) then
          if( Ylm_comp .and. mf /= 0 ) then
            cfad = cfac * Yc(mf,Ylmso(ib,lmf0),Ylmso(ib,lmf0-2*mf))
          else
            cfad = cfac * cmplx( Ylmso(ib,lmf0), 0._db, db )
          endif
          if( E_comp ) then
            cfad = cfad * Bessel(ib,L,ispq)
          else
            cfad = cfad * Besselr(ib,L,ispq)
          endif
          smr(Lmf,ii) = smr(Lmf,ii) - real( cfad, db )
          smi(Lmf,ii) = smi(Lmf,ii) - aimag( cfad )
        else
          fad = fac * Ylmso(ib,lmf0)
          smr(Lmf,ii) = smr(Lmf,ii) - fad * Besselr(ib,L,ispq)
        endif

      end do
    end do

  else

! Expansion in the atom
    ia = -i

    Ln = ii - ianewL(ia)
    Lm = Ln
     
    L = Lato(Lm,ia,igrph)
    isoo = iato(Lm,ia,igrph)
    Lms = L**2 + L + 1 + mato(Lm,ia,igrph) 
    
    do ispt = 1,nspino
      ! Without Spin_axis_free, one expands on spin = isoo
      if( .not. Spin_axis_free .and. isoo /= ispt ) cycle 
!      if( isoo /= ispt ) cycle 
      
      if( Spinorbite ) then
        ispint = min( ispt, nspinp )  ! nspinp = 2 in this case
      else
        ispint = min( ispin, nspinp )
      endif

      if( .not. Spin_axis_free ) then
        m = mato(Lm,ia,igrph) + ispt - isoo    ! it is the m of the harmonics
        if( m > L .or. m < -L ) cycle ! in fact does not happen because without spin_axis_free, ispt = isoo
        Lm0 = L**2 + L + 1 + m    
      endif

      do ib = 1,nbordf(ia)

        j = ibord(ib,ia)
        isym = abs( isbord(ib,ia) )
        jj = new(j) - nspino + ispt

        if( Cal_comp ) then
      
          if( Spin_axis_free ) then
            Ycomp = Ylms(ib,Lms,ispt,isoo,ia)
          elseif( Ylm_comp ) then
            Ycomp = Yc(m,Ylmato(ib,Lm0,ia),Ylmato(ib,Lm0-2*m,ia))
          else
            Ycomp = cmplx( Ylmato(ib,Lm0,ia), 0._db, db )
          endif
        
! complex conjugate on Character because it is the inverse symmetry to go from FDM grid to the atom
          cfac = poidsa(ib,ia) * conjg( Ycomp ) * cmplx( Kar(isym,ispt), -Kari(isym,ispt), db ) / rvol(j)

          abvr( jj ) = abvr( jj ) - real( cfac, db )
          abvi( jj ) = abvi( jj ) - aimag( cfac )

        else
      
          f = poidsa(ib,ia) * Kar(isym,ispt) / rvol(j)
          abvr( jj ) = abvr( jj ) - f * Ylmato(ib,Lm0,ia)
  
        endif

      end do

    end do
    
    do ispt = 1,nspino
      if( .not. Spin_axis_free .and. isoo /= ispt ) cycle 
!      if( isoo /= ispt ) cycle 
      
      if( Spinorbite ) then
        ispint = min( ispt, nspinp )  ! nspinp = 2 in this case
      else
        ispint = min( ispin, nspinp )
      endif

      if( .not. Spin_axis_free ) then
        m = mato(Lm,ia,igrph) + ispt - isoo    ! it is the m of the harmonics, in fact ispt = isoo here 
        if( m > L .or. m < -L ) cycle ! does not happen
        Lm0 = L**2 + L + 1 + m
      endif

      do Lmp = 1,nlmsa(ia)
        Lp = Lato(Lmp,ia,igrph)
        isop = iato(Lmp,ia,igrph)
        Lmsp = Lp**2 + Lp + 1 + mato(Lmp,ia,igrph)
        if( .not. Spin_axis_free ) then
          mp = mato(Lmp,ia,igrph) + ispt - isop
          if( mp > Lp .or. mp < -Lp ) cycle
        endif
  
        charmr = 0._db
        if( Cal_comp ) charmi = 0._db

        do ib = 1,nbordf(ia)
 
          if( Spin_axis_free ) then
            cfac = poidsa(ib,ia) * conjg( Ylms(ib,Lms,ispt,isoo,ia) ) * Phiato_c(ib,Lmsp,ispint,isop,ia)
          elseif( Ylm_comp .and. nphiato_c > 0  ) then
            cfac = poidsa(ib,ia) * conjg( Yc(m, Ylmato(ib,Lm0,ia), Ylmato(ib,Lm0-2*m,ia) ) ) &
                                 * Phiato_c(ib,Lmsp,ispint,isop,ia)
          elseif( Ylm_comp ) then
            cfac = poidsa(ib,ia) * conjg( Yc(m, Ylmato(ib,Lm0,ia), Ylmato(ib,Lm0-2*m,ia) ) ) &
                                 * Phiato(ib,Lmsp,ispint,isop,ia)
          elseif( nphiato_c > 0 ) then 
            cfac = poidsa(ib,ia) * Ylmato(ib,Lm0,ia) * Phiato_c(ib,Lmsp,ispint,isop,ia)                      
          else
            fac = poidsa(ib,ia) * Ylmato(ib,Lm0,ia) * Phiato(ib,Lmsp,ispint,isop,ia)
          endif

          if( Cal_comp ) then
            charmr = charmr + real( cfac, db )
            charmi = charmi + aimag( cfac )
          else
            charmr = charmr + fac
          endif
  
        end do

        Ljj = ianew(ia) + Lmp
        abvr(Ljj) = abvr(Ljj) + charmr
        if( Cal_comp ) abvi(Ljj) = abvi(Ljj) + charmi

      end do  ! end of loop on Lmp
  
    end do
      
  endif

!---- Writing of the matrix -------------------------

  if( icheck > 2 ) then
    boucle_jj: do jj = Lb1r,Lb2r
      j = jj - Lb1r + 1
      if( Cal_comp ) then
        if( ( abs(abvr(jj)) < eps6 ) .and. ( abs(abvi(jj)) < eps6 ) ) then
          mletl(j) = ' .'
          cycle
        endif
      else
        if( abs(abvr(jj))  <  eps6 ) then
          mletl(j) = ' .'
          cycle
        endif
      endif
      if( abs( abvr(jj) - 1 ) < eps6 ) then
        mletl(j) = ' 1'
      elseif( abs(abvr(jj) + 1) < eps6 ) then
        mletl(j) = '-1'
      else
        do k = 1,nlet
          if( Cal_comp ) then
            if( abs(abvr(jj) - vletr(k)) > eps6 .or. abs(abvi(jj) - vleti(k)) > eps6 ) cycle
          else
            if( abs(abvr(jj) - vletr(k)) > eps6 ) cycle
          endif
          mletl(j) = mlet(k)
          cycle boucle_jj
        end do
        if( nlet == nletm**3 ) then
          mletl(j) = ' ? '
        else
          nlet = nlet + 1
          letu = letu + 1
          if( letu > nletm ) then
            letu = 1
            letd = letd + 1
            if( letd > nletm ) then
              letd = 1
              letc = letc + 1
            endif
          endif
          
          mlet(nlet) = let(letc) // let(letd) // let(letu)
          mletl(j) = mlet(nlet)
          vletr(nlet) = abvr(jj)
          if( Cal_comp ) vleti(nlet) = abvi(jj)
        endif
      endif
    end do boucle_jj

    write(3,130) ii, i, Lb1r, Lb2r
    write(3,140) ( mletl(j), j = 1,Lb2r - Lb1r + 1 )

    if( icheck > 3 ) then
      allocate( jnzero(Lb2r - Lb1r + 1) )
      n = 0
      do j = Lb1r,Lb2r
        if( Cal_comp ) then
          if( abs( abvr(j) ) +  abs( abvi(j) ) < eps10 ) cycle
         else
          if( abs( abvr(j) ) < eps10 ) cycle
        endif
         n = n + 1
         jnzero(n) = j 
      end do        
      if( Cal_comp ) then
        write(3,'(300(i5,i4,2e12.4))') ( jnzero(j), newinv(jnzero(j)), abvr(jnzero(j)), abvi(jnzero(j)), j = 1,n )
      else
        write(3,'(300(i5,i4,e12.4))') ( jnzero(j), newinv(jnzero(j)), abvr(jnzero(j)), j = 1,n )
      endif
      deallocate( jnzero )
    endif
    if( sum( abs(smr(:,ii)) ) > eps10 ) then
      if( Cal_comp ) then
        write(3,150) ( smr(Lms,ii), smi(Lms,ii), Lms = 1,nlmso )
      else
        write(3,155) ( smr(Lms,ii), Lms = 1,nlmso )
      endif
    endif
  endif

  if( icheck > 2 .and. ( ( ii == 1 .and. .not. MUMPS ) .or. ( ii == nligne .and. MUMPS ) ) ) then
    write(3,*) 
    if( .not. Cal_comp ) then
      write(3,160) ( mlet(k), vletr(k), k = 1,nlet )
    else
      write(3,165) ( mlet(k), vletr(k), vleti(k), k = 1,nlet )
    endif
  endif
    
  return
  120 format(/'  FDM Matrix:')
  130 format(/' Line =',i5,', i =',i5,', Lb1, Lb2 =',2i5)
  140 format(1000(1x,a3))
  150 format('   sm =',1p,250(1x,2e11.3))
  155 format('   sm =',1p,500e11.3)
  160 format(8(a3,' =',1pe10.3,1x))
  165 format(1p,4(a3,' =',2e10.3,1x))
  end

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

! Function giving the complex harmonics from the real harmonics

function Yc(m,Ylm1,Ylm2)

  use declarations
  implicit none

  integer:: m
  
  complex(kind=db):: Yc

  real(kind=db):: Ylm1, Ylm2

  if( m == 0 ) then
    Yc = cmplx( Ylm1, 0._db, db)
  else
    if( m < 0 ) then
      Yc = cmplx( Ylm2,-Ylm1, db ) * sqrt_1o2
    else
      Yc = cmplx( Ylm1, Ylm2, db ) * ( (-1)**m ) * sqrt_1o2
    endif
  endif

  return
end

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

! Calculation of the new indexation of the FDM points, that is their line number in the general FDM matrix.

subroutine Newind(distai,ianew,ianewL,ibord,icheck,isrt,ivois,lb1,lb2,mpirank0,natome, &
                nbordf,nbtm,new,newinv,nligne,nlmsa,nlmso_eq,npoint,nsort,nsortf,nspino,npsom,numia, &
                nvois,xyz)

  use declarations
  implicit none
  
  integer:: i, ia, ib, icheck, ii, ipr, iv, j, L, lb11, lb22, Lm, mpirank0, n, natome, &
            nbtm, nligne, nlmso_eq, npf, npf1, npoint, nsort, nsortf, nspino, npsom, nvois 

  integer, dimension(natome):: ianew, ianewL, nlmsa
  integer, dimension(nligne):: lb1, lb2, newinv
  integer, dimension(nsort):: isrt
  integer, dimension(natome):: nbordf
  integer, dimension(npsom):: numia
  integer, dimension(0:npoint):: new, newL
  integer, dimension(nbtm,natome):: ibord
  integer, dimension(npsom,nvois):: ivois
  
  real(kind=db), dimension(natome):: distai
  real(kind=db), dimension(4,npsom):: xyz

  ii = 0
  do i = 1,npoint
    if( numia(i) == 0 ) then
      ii = ii + nspino
      new(i) = ii
    else
      new(i) = 0
    endif
  end do

  newL(:) = new(:)
  ianewL(:) = ianew(:)
  
  npf = 1
  npf1 = 0
  do ia = 1,natome
    do i = npf,npoint
      if( numia(i) == 0 .and. xyz(4,i) >= distai(ia) ) exit
    end do
    npf = i
    if( npf1 == npf ) then
      ianew(ia) = ianew(ia-1) + nlmsa(ia-1)
      ianewL(ia) = ianewL(ia-1) + nlmsa(ia-1)
    else
      ianew(ia) = new(npf) - nspino
      ianewL(ia) = newL(npf) - nspino
    endif
    npf1 = npf
    do i = npf,npoint
      if( numia(i) == 0 ) new(i) = new(i) + nlmsa(ia)
      if( numia(i) == 0 ) newL(i) = newL(i) + nlmsa(ia)
    end do
  end do

! Calculation of the positions of the coefficients in the matrix.
  do ii = 1,nligne
    lb1(ii) = ii
    lb2(ii) = ii
  end do

  n = nspino - 1
  do i = 1,npoint
    if( numia(i) /= 0 ) cycle
    ii = newL(i)
    do iv = 1,nvois
      j = ivois(i,iv)
      if( j == 0 .and. mpirank0 == 0 ) then
        call write_error
        do ipr = 6,9,3
          write(ipr,105) i, iv
        end do
        stop
      endif
      ia = numia(j)
      if( ia == 0 ) then
        lb1(ii) = min(lb1(ii),new(j)-n)
        lb2(ii) = max(lb2(ii),new(j))
      elseif( ia > 0 ) then
        lb1(ii) = min(lb1(ii),ianew(ia) + 1)
        lb2(ii) = max(lb2(ii),ianew(ia) + nlmsa(ia))
      elseif( ia == -2 ) then
        lb2(ii) = nligne
      endif
    end do
    if( nspino == 2 ) then
      lb2(ii-1) = lb2(ii)
      lb1(ii-1) = lb1(ii)
    endif
  end do

  do ia = 1,natome
    lb11 = ianew(ia) + 1
    lb22 = ianew(ia) + nlmsa(ia)
    
    do ib = 1,nbordf(ia)
      lb11 = min( lb11, new( ibord(ib,ia) ) - n )
      lb22 = max( lb22, new( ibord(ib,ia) ) )
    end do

    do L = 1,nlmsa(ia)
      ii = ianewL(ia) + L
      lb1(ii) = lb11
      lb2(ii) = lb22
    end do

  end do

  if( nlmso_eq > 0 ) then
    lb11 = nligne - nlmso_eq + 1
    do ib = 1,nsortf
      lb11 = min(lb11,new(isrt(ib))-n)
    end do
    do ii = nligne-nlmso_eq+1,nligne
      lb1(ii) = lb11
      lb2(ii) = nligne
    end do
  endif

  do ii = nligne-1,1,-1
    lb1(ii) = min(lb1(ii),lb1(ii+1))
  end do

  newinv = 0
  do i = 1,npoint
    ii = newL(i)
    if( ii /= 0 ) then
      newinv(ii) = i
      if( nspino == 2 ) newinv(ii-1) = i
    endif
  end do
  do ia = 1,natome
    do Lm = 1,nlmsa(ia)
      newinv(ianewL(ia)+Lm) = - ia
    end do
  end do

  if( icheck > 1 ) then
    write(3,110)
    write(3,120) nligne
    write(3,'(a13,i3)') '   nlmso_eq =', nlmso_eq
    write(3,'(/a13,12i6)') ' ianew(ia)  =', ianew(1:natome)
    write(3,'(a13,12i6)') ' ianewL(ia) =', ianewL(1:natome)
    write(3,'(a13,12i6)') ' nlmsa(ia)  =', nlmsa(1:natome)
  endif
  if( icheck > 2 ) then
    write(3,130)
    write(3,140) (i, new(i), i = 1,npoint)
    write(3,150)
    write(3,140) (i, newL(i), i = 1,npoint)
    write(3,160)
    write(3,170) (ii, lb1(ii), lb2(ii), newinv(ii), ii = 1,nligne)
  endif

  return
  105 format(//' Error in newind for i, iv =',2i6)
  110 format(/' ---- Newind -------',100('-'))
  120 format(/' Number of line =',i6,//)
  130 format(/1x,10('    i  new(i)'))
  140 format(10(i7,i6))
  150 format(/1x,10('    i newL(i)'))
  160 format(/2(2x,'    ii   lb1   lb2   newinv'))
  170 format(2(2x,3i6,i9))
end

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

! Subroutine for inversion of complex general matrix

subroutine invcomp(Matrix_info,n,mat,nm,lwork,is,Stop_job)

  use declarations
  implicit none

  integer:: i, info, ipr, is, j, lwork, n, nm

  integer, dimension(nm):: ipiv

  character(len=Length_name):: Matrix_info

  complex(kind=db), dimension(lwork):: work
  complex(kind=db), dimension(nm,nm):: Mat
  
  complex(kind=db), dimension(:,:), allocatable:: Mat_o

  logical:: Stop_job
  
  allocate( Mat_o(nm,nm) )

  Mat_o(:,:) = Mat(:,:)
  
  if( is == 0 ) then
!cray      call cgetrf(n,n,mat,nm,ipiv,info)
    call zgetrf(n,n,mat,nm,ipiv,info)
  elseif( is == 1 ) then
!cray      call csytrf('u',n,mat,nm,ipiv,work,lwork,info)
    call zsytrf('u',n,mat,nm,ipiv,work,lwork,info)
  else
    call zhetrf('u',n,mat,nm,ipiv,work,lwork,info)
  endif

  if( info /= 0 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,110) info
    end do

    write(9,'(/A)') Matrix_info    
    write(9,'(/a6,i5)') ' dim =', nm
    do i = 1,nm
      write(9,'(1p,1000(1x,2e13.5))') Mat_o(i,:)
    end do
    
    if( Stop_job ) stop
    Stop_job = .true.
  endif

  if( is == 0 ) then
!cray        call cgetri(n,mat,nm,ipiv,work,lwork,info)
    call zgetri(n,mat,nm,ipiv,work,lwork,info)
  elseif( is == 1 ) then
!cray      call csytri('u',n,mat,nm,ipiv,work,info)
    call zsytri('u',n,mat,nm,ipiv,work,info)
    do i = 1,nm
      do j = i + 1,nm
        mat(j,i) = mat(i,j)
     end do
    end do
  else
    call zhetri('u',n,mat,nm,ipiv,work,info)
    do i = 1,nm
      do j = i + 1,nm
        mat(j,i) = conjg( mat(i,j) )
     end do
    end do
  endif

  if( info /= 0 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,120) info
    end do

    write(9,'(/A)') Matrix_info    
    write(9,'(/a6,i5)') ' dim =', nm
    do i = 1,nm
      write(9,'(1p,1000(1x,2e13.5))') Mat_o(i,:)
    end do
    
    if( Stop_job ) stop
    Stop_job = .true.
  endif

  deallocate( Mat_o )

  return
  110 format(/' CGETRF : info =',i5)
  120 format(/' CGETRI : info =',i5)
end

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

! Sub routine for inversion of a symetric real matrix

subroutine invreel(n,mat,nm,lwork)

  use declarations
  implicit none

  integer:: i, info, ipr, j, lwork, n, nm

  integer, dimension(nm):: ipiv

  real(kind=db), dimension(lwork):: work
  real(kind=db), dimension(nm,nm):: mat

  call dsytrf('u',n,mat,nm,ipiv,work,lwork,info)

  if( info /= 0 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,110) info
    end do
    stop
  endif

  call dsytri('u',n,mat,nm,ipiv,work,info)
  do i = 1,nm
    do j = i + 1,nm
      mat(j,i) = mat(i,j)
   end do
  end do

  if( info /= 0 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,120) info
    end do
    stop
  endif

  return
  110 format(/' DSYTRF : info =',i5)
  120 format(/' DSYTRI : info =',i5)
end

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

subroutine Write_ampl(Basereel,Cal_comp,iaabsi,ianew,iato,igrph,irep_util,iso,ispinin,Lato,Lso, &
                      mato,mso,natome,newinv,ngrph,nim,nligne,nlmsa,nlmsam,nlmso_reg,norm,norm_t,npsom, &
                      nso1,nspino,rvol,smi,smr,Spinorbite,xyz)

  use declarations
  implicit none
  
  integer:: i, i1, i2, ia, iaabsi, idir1, idir2, idir3, igrph, ii, ipas, isens, is1, isp, ispinin, jj, &
            jsp, L1, Lm, lm1, lm2, Lmf, Lms, m1, natome, ngrph, nim, nligne, nlmsam, nlmso_reg, nlmw, npsom, &
            nso1, nspino

  integer, dimension(natome):: ianew, nlmsa
  integer, dimension(nligne):: newinv  
  integer, dimension(nrepm,2):: irep_util
  integer, dimension(nso1,ngrph):: Lso, mso, iso
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato
  
  character(len=1):: mot1

  complex(kind=db):: ampl1, ampl2, cfac
  complex(kind=db), dimension(nlmso_reg):: smc
  complex(kind=db), dimension(:), allocatable:: Tau_out
  complex(kind=db), dimension(nlmso_reg,nlmso_reg):: norm, norm_t
  
  logical:: Basereel, Cal_comp, Spinorbite

  real(kind=db), dimension(nim):: rvol
  real(kind=db), dimension(4,npsom):: xyz
  real(kind=db), dimension(nlmso_reg,nligne):: smi, smr

  if( nlmso_reg == 0 ) return
  
  if( Spinorbite ) then
    write(3,'(/A)') ' Wave function, versus (l_out,m_out,s_out):'
  else
    write(3,'(/A)') ' Wave function, versus (l_out,m_out):'
  endif
  select case(mso(1,igrph))
    case(-1,-2,-3)
      idir1 = 1; idir2 = 3; idir3 = 2
      mot1 = 'y'
    case(0)
      idir1 = 1; idir2 = 2; idir3 = 3
      mot1 = 'z'
    case(1,2,3)
      idir1 = 2; idir2 = 3; idir3 = 1
      mot1 = 'x'
  end select
  nlmw = min(nlmso_reg,9)
  if( Cal_comp .or. Basereel ) then
    if( Spinorbite ) then
      write(3,150) mot1, (Lso(Lm,igrph), mso(Lm,igrph), iso(Lm,igrph), Lm = 1,nlmw )
    else
      write(3,160) mot1, (Lso(Lm,igrph), mso(Lm,igrph), Lm = 1,nlmw )
    endif
  else
    if( Spinorbite ) then
      write(3,170) mot1, (Lso(Lm,igrph), mso(Lm,igrph), iso(Lm,igrph), Lm = 1,nlmw )
    else
      write(3,180) mot1, (Lso(Lm,igrph), mso(Lm,igrph), Lm = 1,nlmw )
    endif
  endif
 
  do jsp = 1,nspino
    if( Spinorbite ) write(3,'(a6,i2)') ' isp =', jsp
    do isens = 1,2
      if( isens == 1 ) then
        i1 = nligne  
        i2 = 1
        ipas = - 1
      else
        i1 = 1
        i2 = nligne
        ipas = 1
      endif
      do ii = i1,i2,ipas
        i = newinv(ii)
        if( i <= 0 ) cycle
        if( Spinorbite ) then
          isp = mod(ii+1,2) + 1
          if( isp /= jsp ) cycle
          if( irep_util(igrph,isp) == 0 ) cycle
        endif
        if( abs(xyz(idir1,i)) > eps10 .or. abs(xyz(idir2,i)) > eps10 ) cycle
        if( ( isens == 1 .and. xyz(idir3,i) > - eps10 ) .or. ( isens == 2 .and. xyz(idir3,i) < - eps10 ) ) cycle
        if( Cal_comp ) then
          if( Basereel ) then
            do Lmf = 1,nlmw
              smc(Lmf) = sum( norm(Lmf,1:nlmso_reg) * cmplx( smr(1:nlmso_reg,ii), smi(1:nlmso_reg,ii), db ))
            end do
            write(3,190) xyz(idir3,i)*bohr, smc(1:nlmw) / rvol(i)
          else
            write(3,200) xyz(idir3,i)*bohr, ( smr(Lms,ii) / rvol(i), smi(Lms,ii) / rvol(i), Lms = 1,nlmw )
          endif
        else
          if( Basereel ) then
            do Lmf = 1,nlmw
              smc(Lmf) = sum( norm(Lmf,1:nlmso_reg) * smr(1:nlmso_reg,ii) )
            end do
            write(3,190) xyz(idir3,i)*bohr, smc(1:nlmw) / rvol(i)
          else
            write(3,200) xyz(idir3,i)*bohr, smr(1:nlmw,ii) / rvol(i)
          endif
        endif
      end do
    end do
    
  end do

  write(3,'(/A)') ' Absorbing atom amplitudes, diagonal term versus outer sphere basis'
  ia = iaabsi
  if( Spinorbite ) then
    write(3,'(18x,18(3i3,12x))') (Lato(Lm,ia,igrph), mato(Lm,ia,igrph), iato(Lm,ia,igrph), Lm = 1,nlmsa(ia) )
  else
    write(3,'(18x,18(2x,2i3,15x))') (Lato(Lm,ia,igrph), mato(Lm,ia,igrph), Lm = 1,nlmsa(ia) )
  endif

  allocate( Tau_out(nlmsa(ia)) )
  Tau_out(:) = (0._db,0._db)

  do lm1 = 1,nlmsa(ia)
    L1 = Lato(lm1,ia,igrph)
    m1 = mato(lm1,ia,igrph)
    if( Spinorbite ) then
      is1 = iato(lm1,ia,igrph)
    else
      is1 = ispinin
    endif

    Tau_out(:) = (0._db,0._db)

    do lm2 = 1,nlmsa(ia)
      jj = ianew(ia) + lm2

      do Lmf = 1,nlmso_reg
        if( Lso(Lmf,igrph) == L1 .and. mso(Lmf,igrph) == m1 .and. iso(Lmf,igrph) == is1 ) exit
      end do
      if( Lmf > nlmso_reg ) cycle    
      
      if( Basereel ) then
        if( Cal_comp ) then
          ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,jj), smi(1:nlmso_reg,jj), db ))
        else
          ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,jj) )
        endif
      else
        ampl2 = cmplx(smr(Lmf,jj),smi(Lmf,jj),db)
      endif

      Tau_out(lm2) = - ampl2 

    end do
    write(3,250) Lato(lm1,ia,igrph), mato(lm1,ia,igrph), iato(lm1,ia,igrph), Tau_out(:)
  end do

  write(3,'(/A)') ' Absorbing atom multiple scattering amplitudes'
  ia = iaabsi
  if( Spinorbite ) then
    write(3,'(18x,18(3i3,12x))') (Lato(Lm,ia,igrph), mato(Lm,ia,igrph), iato(Lm,ia,igrph), Lm = 1,nlmsa(ia) )
  else
    write(3,'(18x,18(2x,2i3,15x))') (Lato(Lm,ia,igrph), mato(Lm,ia,igrph), Lm = 1,nlmsa(ia) )
  endif

  do lm1 = 1,nlmsa(ia)
    ii = ianew(ia) + lm1

    Tau_out(:) = (0._db,0._db)

    do lm2 = 1,nlmsa(ia)
      jj = ianew(ia) + lm2
      
      cfac = (0._db, 0._db )

      do Lmf = 1,nlmso_reg

        if( Basereel ) then
          if( Cal_comp ) then
            ampl1 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,ii), smi(1:nlmso_reg,ii), db ))
            ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * cmplx( smr(1:nlmso_reg,jj), smi(1:nlmso_reg,jj), db ))
          else
            ampl1 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,ii) )
            ampl2 = sum( norm_t(1:nlmso_reg,Lmf) * smr(1:nlmso_reg,jj) )
          endif
        else
          ampl1 = cmplx(smr(Lmf,ii),smi(Lmf,ii),db)
          ampl2 = cmplx(smr(Lmf,jj),smi(Lmf,jj),db)
        endif

        cfac = cfac + conjg( ampl1 ) * ampl2
     
      end do

! One multiplies by -img in order -aimag(Taull) is the density of state
! Conjugate to get the real part same sign that Green
      cfac = - img * conjg( cfac )
 
      Tau_out(lm2) = cfac

    end do
    write(3,250) Lato(lm1,ia,igrph), mato(lm1,ia,igrph), iato(lm1,ia,igrph), Tau_out(:)

  end do

  write(3,'(/A)') ' Cluster outgoing scattering amplitude'
  if( Spinorbite ) then
    write(3,'(18x,18(3i3,12x))') (Lso(Lm,igrph), mso(Lm,igrph), iso(Lm,igrph), Lm = 1,nlmso_reg )
  else
    write(3,'(18x,18(2x,2i3,15x))') (Lso(Lm,igrph), mso(Lm,igrph), Lm = 1,nlmso_reg )
  endif

  deallocate( Tau_out )

  do lm1 = 1,nlmso_reg
    smc(:) = (0._db,0._db)
    do lm2 = 1,nlmso_reg
    
      jj = nligne - nlmso_reg + lm2
       
      if( Basereel ) then
        if( Cal_comp ) then
          ampl1 = sum( norm_t(1:nlmso_reg,lm1) * cmplx( smr(1:nlmso_reg,jj), smi(1:nlmso_reg,jj), db ))
        else
          ampl1 = sum( norm_t(1:nlmso_reg,lm1) * smr(1:nlmso_reg,jj) )
        endif
      else
        ampl1 = img * cmplx(smr(lm1,jj),smi(lm1,jj),db)
      endif
      smc(lm2) = ampl1      
    end do
    write(3,250) Lso(lm1,igrph), mso(lm1,igrph), iso(lm1,igrph), smc(:)
  end do
 
  return
  150 format(4x,a1,2x,18(7x,3i3,7x))
  160 format(4x,a1,2x,18(9x,2i3,8x))
  170 format(4x,a1,2x,18(3x,3i3,3x))
  180 format(4x,a1,2x,18(2x,2i3,3x))
  190 format(f7.3,1p,18(1x,2e11.3))
  200 format(f7.3,1p,36e11.3)
  250 format(3i3,2x,1p,18(1x,2e11.3))
end

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

! Recopie des amplitudes pour les (L,m) non calculees

subroutine recop_taull(nlmagm,nspinp,Sym_cubic,Taull,Ylm_comp)

  use declarations
  implicit none

  integer:: i2, i3, i4, i5, i6, i7, i8, i9, j2, j3, j4, j5, j6, j7, j8, j9, nlmagm, nspinp

  logical:: Sym_cubic, Ylm_comp

  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp):: Taull

  i2 = 2;  j2 = 2
  i3 = 3;  j3 = 3
  i4 = 4;  j4 = 4
  i5 = 5;  j5 = 5
  i6 = 6;  j6 = 6
  i7 = 7;  j7 = 7
  i8 = 8;  j8 = 8
  i9 = 9;  j9 = 9
     
  if( Sym_cubic ) then
! L == 1
    if( nlmagm > 1 ) then
      Taull(i2,:,j2,:) = Taull(i3,:,j3,:)
      Taull(i4,:,j4,:) = Taull(i3,:,j3,:)
    endif

! L == 2
    if( nlmagm > 4 ) then
      if( Ylm_comp ) then
        if( abs( Taull(i8,1,j8,1) ) > eps10 ) then
          Taull(i5,:,j5,:) = 0.5_db*( Taull(i7,:,j7,:) + Taull(i8,:,j8,:))
          Taull(i9,:,j9,:) = Taull(i5,:,j5,:)
          Taull(i5,:,j9,:) = 0.5_db*( Taull(i7,:,j7,:) - Taull(i8,:,j8,:))
          Taull(i9,:,j5,:) = Taull(i5,:,j9,:)
        else
          Taull(i8,:,j8,:) = 2 * Taull(i5,:,j5,:) - Taull(i7,:,j7,:)
          Taull(i6,:,j6,:) = Taull(i8,:,j8,:)
        endif
      else
        Taull(i9,:,j9,:) = Taull(i7,:,j7,:)
        Taull(i6,:,j6,:) = Taull(i5,:,j5,:)
        Taull(i8,:,j8,:) = Taull(i5,:,j5,:)
      endif
    endif

  elseif( .not. Ylm_comp ) then
! sym_4
    if( nlmagm > 1 ) Taull(i4,:,j4,:) = Taull(i2,:,j2,:)
    if( nlmagm > 4 ) Taull(i6,:,j6,:) = Taull(i8,:,j8,:)

  endif
    
  return
end

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

! Recopie des amplitudes pour les (L,m) non calculees

! not used !

subroutine recop_ampl(nlm_rixs,nlmomax,nspinp,Sym_cubic,Ampl,Ylm_comp)

  use declarations
  implicit none

  integer:: nlm_rixs, nlmomax, nspinp

  logical:: Sym_cubic, Ylm_comp

  complex(kind=db), dimension(nlm_rixs,nspinp,nlmomax,nspinp):: Ampl

  if( Sym_cubic ) then
! L == 1
    if( nlm_rixs > 1 ) then
      Ampl(2,:,2,:) = Ampl(3,:,3,:)
      Ampl(4,:,4,:) = Ampl(3,:,3,:)
    endif

! L == 2
    if( nlm_rixs > 4 ) then
      if( Ylm_comp ) then
        if( abs( Ampl(8,1,8,1) ) > eps10 ) then
          Ampl(5,:,5,:) = 0.5_db*( Ampl(7,:,7,:) + Ampl(8,:,8,:))
          Ampl(9,:,9,:) = Ampl(5,:,5,:)
          Ampl(5,:,9,:) = 0.5_db*( Ampl(7,:,7,:) - Ampl(8,:,8,:))
          Ampl(9,:,5,:) = Ampl(5,:,9,:)
        else
          Ampl(8,:,8,:) = 2 * Ampl(5,:,5,:) - Ampl(7,:,7,:)
          Ampl(6,:,6,:) = Ampl(8,:,8,:)
        endif
      else
        Ampl(9,:,9,:) = Ampl(7,:,7,:)
        Ampl(6,:,6,:) = Ampl(5,:,5,:)
        Ampl(8,:,8,:) = Ampl(5,:,5,:)
      endif
    endif

  elseif( .not. Ylm_comp ) then
! sym_4
    if( nlm_rixs > 1 ) Ampl(4,:,4,:) = Ampl(2,:,2,:)
    if( nlm_rixs > 4 ) Ampl(6,:,6,:) = Ampl(8,:,8,:)

  endif

  return
end

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

! Calculation of multiple scattering amplitude using the multiple scattering theory

subroutine Mat_MST(Axe_atom_gr,Cal_xanes,Dist_coop,Ecinetic,Eimag,Full_atom,ia_coop,ia_eq,ia_rep, &
                    iaabsi,iaprotoi,iato,icheck,igroupi,igrph,iopsymr,irep_util,is_eq,ispin,karact,Lato,Lmaxa,Lmax_g,mato, &
                    n_atom_0,n_atom_coop,n_atom_ind,n_atom_proto,nab_coop,natome,natomp,nb_eq,nb_rpr,nb_rep_t,nb_sym_op,nb_test, &
                    ngroup_m,ngrph,nlmagm,nlmsa,nlmsam,nlmsamax,Normaltau,nspin,nspino,nspinp,pos,posi,recop,Repres_comp, &
                    Rmtg,rot_atom,Solsing,Spin_axis_free_L,Spinorbite,State_all_r,Sym_cubic,Tau_ato,Tau_nondiag,Tau_coop,Taull, &
                    Test_case,Time_fill,Time_tria,Ylm_comp)

  use declarations
  implicit none

  integer, parameter:: nletm = 52

  integer:: i, I_rep_dim, i0, i1, i2, ia, iaabsi, iapr, icheck, ib, ib1, ibb, icolone, icolone0, iga, igr, igrph, &
    iligne, iligne0, ind, irep, is1, is2, isg, iso_p, iso1, iso2, isp, ispin, its, j0, jligne, jligne0, jsp, &
    L, L1, L2, LL, Lm, Lm_p, Lm0, Lm1, Lm10, Lm1p, Lm01, Lm01c, Lm02, Lm02c, Lm2, Lm20, Lm2p, Lma, Lmax, Lmax_g, &
    Lmb, Lmc, Lmd1, Lmd2, m, m_p, m1, m2, ms1, ms2, &
    n_atom_0, n_atom_ind, n_atom_proto, n1, n_atom_coop, nab_coop, natome, natomp, nb_sym_op, nb_test, ndim, &
    ngroup_m, ngrph, nlm, nlmabs, nlmagm, nlmc, nlmch, nlmr, nlmsam, nlmsamax, nlmsmax, nspin, nspino, nspinp, nspinr

  integer, dimension(n_atom_coop):: ia_coop
  integer, dimension(natome):: iaprotoi, igroupi, Lmaxa, nb_eq, nlmsa
  integer, dimension(nopsm):: iopsymr
  integer, dimension(natome,natome):: nb_rpr
  integer, dimension(nb_sym_op,natome,natome):: ia_rep, nb_rep_t
  integer, dimension(nb_sym_op,natome):: ia_eq, is_eq
  integer, dimension(nrepm,2):: irep_util
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato

  character(len=Length_name):: Matrix_info
   
  complex(kind=db):: fnormc, Gmat, Gmatc, Gmats, kr
  complex(kind=db), dimension(nspin):: konde
  complex(kind=db), dimension(nopsm,nrepm):: karact
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,n_atom_0:n_atom_ind)::tau_ato
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,natome):: Taull
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,nab_coop,nspino):: Tau_coop
  complex(kind=db), dimension(:), allocatable:: bessel, neuman, ylmc
  complex(kind=db), dimension(:,:), allocatable:: hankelm, Mat, MAT_AC, Mat_CA, Mat_tau, Matp, Taullp, Taullq
  complex(kind=db), dimension(:,:,:), allocatable:: Cmatr, Dlmm_ia, Taup, Taup_inv
  complex(kind=db), dimension(:,:,:,:), allocatable :: Taull_tem
  complex(kind=db), dimension(:,:,:,:,:), allocatable:: DLmsLms_sym_AC, DLmsLms_sym_CA, Taull_p
  complex(kind=db), dimension(:,:,:,:,:), allocatable:: Cmat
  
  logical::Brouder, Cal_xanes, Ereel, Full_atom, Normaltau, Solsing, Recop, &
    Repres_comp, Spin_axis_free, Spin_axis_free_L, Spinorbite, State_all_r, Stop_job, Sym_cubic, Sym_3, Tau_nondiag, Ylm_comp
  logical, dimension(nb_test):: Test_case

  real(kind=db):: cosang, Eimag, fac, fnorm, g, Gaunt_r, Gauntcp, Gmatr, r, rkr, tp1, tp2, tp3, Time_fill, Time_tria

  real(kind=db), dimension(2):: Dist_coop 
  real(kind=db), dimension(3):: w
  real(kind=db), dimension(nspin):: Ecinetic, konder
  real(kind=db), dimension(3,3):: Mat_rot
  real(kind=db), dimension(0:n_atom_proto):: Rmtg
  real(kind=db), dimension(3,ngroup_m):: Axe_atom_gr
  real(kind=db), dimension(3,natome):: posi
  real(kind=db), dimension(3,natomp):: pos
  real(kind=db), dimension(3,3,natome):: rot_atom
  real(kind=db), dimension(:), allocatable:: besselr, neumanr, Ylm

  real(kind=sg):: time

  Stop_job = .true. 
  
  Spin_axis_free = Spin_axis_free_L .and. .not. Test_case(5)

! Actually all representation are 1 dimension, sorry.
  I_rep_dim = 1

! Corresponds to Ch. Brouder et al. PRB 54, 7334 (1996) to get a more stable multiple scattering matrix.
! Not coded with spin-orbit
  Brouder = .false.

! Filling of the MST matrix
  call CPU_TIME(time)
  tp1 = real(time,db)

  ndim = sum( nlmsa(1:natome) )
  nlmabs = nlmsa(iaabsi)
  nlmsmax = 0
  do ia = 1,natome
    nlmsmax = max( nlmsmax, nlmsa(ia) )
  end do

  konde(:) = sqrt( cmplx(Ecinetic(:), Eimag, db) )

  if( icheck > 1 ) then
    write(3,110)
    write(3,115) igrph, ndim
  endif

  if( ndim /= 0 ) then

  allocate( taup(nlmsmax,nlmsmax,natome) )
  allocate( taullp(nlmsmax,nlmsmax) )
  allocate( taullq(nlmsmax,nlmsmax) )

  nlmch = ndim
  Taup(:,:,:) = (0._db, 0._db)

  if( Tau_nondiag ) allocate( Taup_inv(nlmsmax,nlmsmax,natome) )

  if( Spin_axis_free .and. natome == 1 ) then

    allocate( Taull_p(nlmagm,nspinp,nlmagm,nspinp,natome) )
    Taull_p(:,:,:,:,:) = (0._db, 0._db)
    do ia = 1,natome
      nlm = ( Lmaxa(ia) + 1 )**2
      if( Full_atom ) then
        iapr = ia 
      else
        iapr = iaprotoi( ia )
      endif
      Taull_p(1:nlm,:,1:nlm,:,ia) = Tau_ato(1:nlm,:,1:nlm,:,iapr) 
    end do
    
    call Rotate_TauLL(icheck,Lmaxa,natome,nlmagm,nspinp,Rot_atom,TauLL_p,.false.)
  endif

  boucle_at: do ia = 1,natome

    if( nspin == 2 .and. .not. Spin_axis_free ) then
      igr = igroupi(ia)
      iga = igroupi(iaabsi)

      cosang = sum( Axe_Atom_gr(:,iga) * Axe_Atom_gr(:,igr) )
      if( abs( cosang - 1 ) < eps4 ) then
        its = 1
      elseif( abs( cosang + 1 ) < eps4 ) then
        its = -1
      else
        its = 0
      endif
    else
      its = 1
    endif

    if( Full_atom ) then
      iapr = ia
    else
      iapr = iaprotoi(ia)
    endif

    do Lm1 = 1,nlmsa(ia)

      L1 = Lato(Lm1,ia,igrph)
      m1 = mato(Lm1,ia,igrph)
      if( nspino == 2 ) then
        i1 = iato(Lm1,ia,igrph)
      else
        i1 = ispin
      endif
      if( its == - 1 ) then
        i1 = 3 - i1
        m1 = - m1
      endif
      Lm10 = L1**2 + L1 + 1 + m1
      do Lm2 = 1,nlmsa(ia)
        L2 = Lato(Lm2,ia,igrph)
        m2 = mato(Lm2,ia,igrph)
        if( nspino == 2 ) then
          i2 = iato(Lm2,ia,igrph)
        else
          i2 = ispin
        endif
        if( its == - 1 ) then
          i2 = 3 - i2
          m2 = - m2
        endif
        Lm20 = L2**2 + L2 + 1 + m2
        if( Spin_axis_free .and. natome == 1 ) then
          Taup(Lm1,Lm2,ia) = Taull_p(Lm10,i1,Lm20,i2,ia)
        elseif( its == 0 ) then
          if( i1 == i2 ) Taup(Lm1,Lm2,ia) = 0.5_db * ( Tau_ato(Lm10,1,Lm20,1,iapr) + Tau_ato(Lm10,nspinp,Lm20,nspinp,iapr) )
        else
          Taup(Lm1,Lm2,ia) = Tau_ato(Lm10,i1,Lm20,i2,iapr)
        endif
      end do
    end do
    
    if( Tau_nondiag .and. nlmsa(ia) > 0 ) then
      nlm = nlmsa(ia)
      allocate( Mat(nlm,nlm) )
      Mat(1:nlm,1:nlm) = Taup(1:nlm,1:nlm,ia)
!      Mat = Transpose( Mat )
      if( Ylm_comp ) then
        Matrix_info = '  Taup matrix  with Ylm_com true from Mat_MST'
        call invcomp(Matrix_info,nlm,mat,nlm,nlm,0,Stop_job)
      else
        Matrix_info = '  Taup matrix  with Ylm_com false from Mat_MST'
        call invcomp(Matrix_info,nlm,mat,nlm,nlm,1,Stop_job)
      endif
      Taup_inv(1:nlm,1:nlm,ia) = Mat(1:nlm,1:nlm)
      deallocate( Mat )
    endif
    
  end do boucle_at

  if( Spin_axis_free .and. natome == 1 ) deallocate( Taull_p )

  allocate( mat(ndim,nlmch) )
  mat(:,:) = (0._db, 0._db)

  konder(1:nspin) = real( konde(1:nspin), db )
  if( abs(Eimag) < eps10 .and. Ecinetic(1) > eps10 .and. Ecinetic(nspin) > eps10 ) then
    Ereel = .true.
  else
    Ereel = .false.
  endif

  if( natome > 1 .or. nb_sym_op > 1 ) then
    if( Spin_axis_free ) then
      allocate( DLmsLms_sym_AC(nlmsamax,-Lmax_g:Lmax_g,nspino,nb_sym_op,natome) )
      allocate( DLmsLms_sym_CA(nlmsamax,-Lmax_g:Lmax_g,nspino,nb_sym_op,natome) )
      DLmsLms_sym_AC(:,:,:,:,:) = ( 0._db, 0._db )
      DLmsLms_sym_CA(:,:,:,:,:) = ( 0._db, 0._db )
      do ia = 1,natome
        Lmax = Lmaxa(ia)
        call DLmsLms_sym_cal(DLmsLms_sym_AC,DLmsLms_sym_CA,ia,iato,icheck,igrph,iopsymr,is_eq,Lato,Lmax,Lmax_g,mato,natome, &
                       nb_eq,nb_sym_op,ngrph,nlmsa,nlmsam,nlmsamax,nspino,Rot_atom)
      end do
    else
      allocate( Cmat(natome,nlmsamax,nb_sym_op,-Lmax_g:Lmax_g,nspino) )
      call Cmat_cal(Cmat,iato,icheck,igrph,iopsymr,irep_util,is_eq,karact,Lato,Lmax_g,mato,natome,nb_eq,nb_sym_op,ngrph, &
                nlmsamax,nlmsa,nlmsam,nspino,Rot_atom,Ylm_comp)
    endif
  endif

  fnorm = 1._db
  fnormc = ( 1._db, 0._db )

  iligne = 0
  icolone = 0

  do ia = 1,natome   ! Loop over line
    iligne0 = sum( nlmsa(1:ia-1) )

    if( nb_sym_op > 1 .or. Normaltau .or. Ylm_comp ) then
      ib1 = 1
    else
      ib1 = ia + 1
    endif

    do ib = ib1,natome   ! Loop over column
      icolone0 = sum( nlmsa(1:ib-1) )
      Lmax = Lmaxa(ia) + Lmaxa(ib)
      if( Ereel ) then
        allocate( neumanr(0:Lmax) )
        allocate( besselr(0:Lmax) )
      else
        allocate( neuman(0:Lmax) )
        allocate( bessel(0:Lmax) )
      endif
      nspinr = min( nspin, nspino )
      allocate( hankelm(0:Lmax,nspinr) )
      nlmr = (Lmax + 1)**2
      nlmc = ( (Lmax + 1) * (Lmax + 2) ) / 2
      allocate( ylmc(nlmc) )
      if( .not. Ylm_comp ) allocate( ylm(nlmr) )

! Loop over atom n1 such a symmetry S exists with:
! S(n0) = n0 and S(n1) = n', where n' is not part of p'.
      do n1 = 1,nb_rpr(ia,ib)
      
! ibb : cluster index not symmetrized
        ibb = ia_rep(n1,ia,ib)

! The atom does not scatter towards itself
        w(1:3) = posi(1:3,ia) - pos(1:3,ibb)
        r = sqrt( sum( w(:)**2 ) )
        if( r < eps10 ) cycle

! Search of the index "ind" corresponding to the representant n1
        do ind = 1,nb_eq(ib)
          if( ibb == ia_eq(ind,ib) ) exit
        end do

        fac = real( nb_eq(ia) * nb_rep_t(n1,ia,ib), db ) / I_rep_dim

! Calculation of Hankel functions
        do isp = 1,nspinr
          if( Spinorbite ) then
            jsp = isp
          else
            jsp = ispin
          endif
          if( Ereel ) then
            rkr = konder(jsp) * r
            call cbessneur(fnorm,rkr,Lmax,Lmax,besselr,neumanr)
            hankelm(:,isp) = cmplx(besselr(:),neumanr(:),db)
          else
            kr = konde(jsp) * r
            call cbessneu(fnormc,kr,Lmax,Lmax,bessel,neuman)
            hankelm(0:Lmax,isp) = bessel(0:Lmax) + img*neuman(0:Lmax)
          endif
          do L = 1,Lmax
            hankelm(L,isp) = img**L * hankelm(L,isp)
          end do
        end do

! Calculation of Ylm
        call cylm(Lmax,w,r,Ylmc,nlmc)
        if( .not. Ylm_comp ) call Ylmcr(Lmax,nlmc,nlmr,Ylmc,Ylm)

        do Lm1 = 1,nlmsa(ia)

          iligne = iligne0 + Lm1
          L1 = Lato(Lm1,ia,igrph)
          m1 = mato(Lm1,ia,igrph)
          iso1 = iato(Lm1,ia,igrph)

          do Lm2 = 1,nlmsa(ib)

            icolone = icolone0 + Lm2

            L2 = Lato(Lm2,ib,igrph)
            m2 = mato(Lm2,ib,igrph)
            iso2 = iato(Lm2,ib,igrph)
            if( iso2 /= iso1 .and. .not. Spin_axis_free ) cycle
      if( iso1 /= iso2 ) cycle
            Lma = abs( L1 - L2 )
            Lmb = L1 + L2

! Loops over is1, ms1, is2 and ms2, associated with symmetry and atom rotation
            Gmats = (0._db, 0._db)

            do is1 = 1,nspino
            do ms1 = -L1,L1

! spin index, while solution index is not necessary the spin, inside the MST matrix, it is the case. 
              isp = min( is1, nspin ) ! valid with and without spinorbit

              if( Spin_axis_free ) then
                if( abs( DLmsLms_sym_AC(Lm1,ms1,is1,1,ia) ) < eps10 ) cycle
              elseif( nb_sym_op == 1 ) then
                if( ms1 /= m1 .or. is1 /= iso1 ) cycle
              else
                if( abs( Cmat(ia,Lm1,1,ms1,is1) ) < eps10 ) cycle
              endif

              do is2 = 1,nspino
              do ms2 = -L2,L2

! There is no spin-flip in the interstitial region
                if( is1 /= is2 ) cycle
                if( Spin_axis_free ) then
                  if( abs( DLmsLms_sym_CA(Lm2,ms2,is2,ind,ib) ) < eps10 ) cycle
                elseif( nb_sym_op == 1 ) then
                  if( ms2 /= m2 .or. is2 /= iso2 ) cycle
                else
                  if( abs( Cmat(ib,Lm2,ind,ms2,is2) ) < eps10) cycle
                endif

                Gmat = ( 0._db, 0._db )
                do L = Lma,Lmb,2

                  if( Ylm_comp ) then
                    Lm0 = ( L**2 + L ) / 2 + 1
                    Gmatc = ( 0._db, 0._db )
                    do m = -L,L
! It is Y(L2,ms2) which is complexe conjugate
                      g = Gauntcp(L2,ms2,L1,ms1,L,m)
                      if( abs( g ) < eps10 ) cycle
                      Lmc = Lm0 + abs(m)
                      if( m >= 0 ) then
                        Gmatc = Gmatc + Ylmc(Lmc) * g
                      else
                        Gmatc = Gmatc + (-1)**m * conjg( Ylmc(Lmc) ) * g
                      endif
                    end do
                    Gmat = Gmat + Gmatc * hankelm(L,isp) 
                  else
                    LL = L**2 + L + 1
                    Gmatr = 0._db
                    do m = -L,L
                      g = Gaunt_r(L1,ms1,L2,ms2,L,m)
                      if( abs( g ) < eps10 ) cycle
                      Lm = LL + m
                      Gmatr = Gmatr + Ylm(Lm) * g
                    end do
                    Gmat = Gmat + Gmatr * hankelm(L,isp)    
                  endif

                end do ! end of loop over L

                Gmat = - quatre_pi * Gmat * img**(L1-L2+1)

                if( Spin_axis_free ) then ! <Atom:Cluster>  <Cluster>   <Cluster:Atom>
                  Gmats = Gmats + DLmsLms_sym_AC(Lm1,ms1,is1,1,ia) * Gmat * DLmsLms_sym_CA(Lm2,ms2,is2,ind,ib)                  
                elseif( nb_sym_op == 1 ) then
                  Gmats = Gmats + Gmat
                else
                  Gmats = Gmats + Cmat(ia,Lm1,1,ms1,is1) * Gmat * Conjg( Cmat(ib,Lm2,ind,ms2,is2) )
                endif

              end do
              end do  ! end of loop over Lms2

            end do
            end do  ! end of loop over Lms1

            mat(iligne,icolone) = mat(iligne,icolone) - fac * Gmats
            if( .not. ( nb_sym_op > 1 .or. Normaltau .or. Ylm_comp ) ) mat(icolone,iligne) = mat(iligne,icolone)

          end do  ! end of loop over Lm2
        end do  ! end of loop over Lm1

      end do ! end of loop over equivalent atoms

      if( Ereel ) then
        deallocate( besselr )
        deallocate( neumanr )
      else
        deallocate( bessel )
        deallocate( neuman )
      endif

      deallocate( hankelm )
      if( .not. Ylm_comp ) deallocate( Ylm )
      deallocate( ylmc )
    
    end do ! Loop over ib

  end do ! Loop over ia

! diagonal
  if( normaltau .or. Brouder ) then
    do iligne = 1,ndim
      mat(iligne,iligne) = mat(iligne,iligne) + (1._db,0._db)
    end do
  else
    iligne = 0
    do ia = 1,natome
      jligne0 = iligne
      do Lm = 1,nlmsa(ia)
        iligne = iligne + 1
        if( tau_nondiag ) then
          do Lm2 = 1,nlmsa(ia)
            jligne = jligne0 + Lm2
            mat(iligne,jligne) = mat(iligne,jligne) + Taup_inv(Lm,Lm2,ia)
          end do
        else
          mat(iligne,iligne) = mat(iligne,iligne) + 1/Taup(Lm,Lm,ia)
        endif
      end do
    end do
  endif

  if( icheck > 2 ) then
    write(3,120)
    iligne = 0
    do ia = 1,natome
      do Lm = 1,nlmsa(ia)
        iligne = iligne + 1
        if( Spinorbite ) then
          write(3,130) iligne, ia, Lato(Lm,ia,igrph), mato(Lm,ia,igrph), iato(Lm,ia,igrph)
        else
          write(3,130) iligne, ia, Lato(Lm,ia,igrph), mato(Lm,ia,igrph), ispin
        endif
        write(3,140) mat(iligne,1:ndim)
      end do
    end do
  endif

  if( normaltau .or. Brouder ) then
    i0 = 0
    do ia = 1,natome
      j0 = 0
      do ib = 1,natome
        if( ia == ib ) then
          j0 = j0 + nlmsa(ib)
          cycle
        endif

        allocate( matp(nlmsa(ia),nlmsa(ib)) )
        do lma = 1,nlmsa(ia)
          i = i0 + lma
          do lmb = 1,nlmsa(ib)
            matp(lma,lmb) = mat(i,j0+lmb)
          end do
        end do

        if( Brouder ) then

          if( Spinorbite ) then
            write(6,'(//3x,A)') 'Brouder not coded with Spinorbite !'
            stop
          else
            do lma = 1,nlmsa(ia)
              i = i0 + lma
              do lmb = 1,nlmsa(ib)
                mat(i,j0+lmb) = sqrt( Taup(lma,lma,ia) ) * matp(lma,lmb) * sqrt( Taup(lmb,lmb,ib) )
              end do
            end do
          endif

        else

          do lma = 1,nlmsa(ia)
            i = i0 + lma
            do lmb = 1,nlmsa(ib)
              mat(i,j0+lmb) = sum( Taup(lma,1:nlmsa(ia),ia) * matp(1:nlmsa(ia),lmb) )
            end do
          end do

        endif

        deallocate( matp )
        j0 = j0 + nlmsa(ib)
      end do
      i0 = i0 + nlmsa(ia)
    end do
  endif

  call CPU_TIME(time)
  tp2 = real(time,db)

  Time_fill = tp2 - tp1

! Filling of matrix is done, now come the resolution

  if( natome == 1 .and. nb_rpr(1,1) == 1 ) then
! When only one atom (and no equivalent by symmetry), no inversion to perform
    ia = 1
    do Lm1 = 1,nlmsa(ia)
      L = Lato(Lm1,ia,igrph)
      m1 = mato(Lm1,ia,igrph)
      Lm01 = L**2 + L + 1 + m1
      if( Spinorbite ) then
        is1 = iato(Lm1,ia,igrph)
      else
        is1 = ispin
      endif
      do Lm2 = 1,nlmsa(ia)
        L = Lato(Lm2,ia,igrph)
        m2 = mato(Lm2,ia,igrph)
        Lm02 = L**2 + L + 1 + m2
        if( Spinorbite ) then
          is2 = iato(Lm2,ia,igrph)
        else
          is2 = ispin
        endif
        Taull(Lm01,is1,Lm02,is2,ia) = Taull(Lm01,is1,Lm02,is2,ia) + Taup(Lm1,Lm2,ia)
        if( .not. Spinorbite .and. Repres_comp ) then
          lm01c = lm01 - 2 * m1
          Lm02c = lm02 - 2 * m2
          isg = (-1)**(m1+m2)
          Taull(Lm02c,is2,lm01c,is1,ia) = Taull(Lm02c,is2,lm01c,is1,ia) + isg * Taup(Lm1,Lm2,ia)
        endif
      end do
    end do

  else

! Full multiple scattering, matrix inversion

    if( Normaltau .or. nb_sym_op > 1 .or. Ylm_comp ) then
      Matrix_info = ' Main MST matrix with - normaltau .or. nb_sym_op > 1 .or. Ylm_comp - from Mat_MST' 
      call invcomp(Matrix_info,ndim,mat,ndim,ndim,0,Stop_job)
    else
      Matrix_info = ' Main MST matrix without - normaltau .or. nb_sym_op > 1 .or. Ylm_comp - from Mat_MST' 
      call invcomp(Matrix_info,ndim,mat,ndim,ndim,1,Stop_job)
    endif

    if( Brouder ) then

      do Lm = 1,ndim
        mat(Lm,Lm) = mat(Lm,Lm) - (1._db,0._db)
      end do

      i0 = 0
      do ia = 1,natome
        j0 = 0
        do ib = 1,natome

          if( Spinorbite ) then
          else
            do lma = 1,nlmsa(ia)
              i = i0 + lma
              do lmb = 1,nlmsa(ib)
                mat(i,j0+lmb) = sqrt( taup(lma,lma,ia) ) * mat(i,j0+lmb) * sqrt( taup(lmb,lmb,ib) )
              end do
              if( ia == ib ) mat(i,i) = mat(i,i) + Taup(lma,lma,ia)
            end do
          endif

          j0 = j0 + nlmsa(ib)
        end do
        i0 = i0 + nlmsa(ia)
      end do

    endif

    Lm0 = 0

    do ia = 1,natome

      if( ia > 1 ) Lm0 = Lm0 + nlmsa(ia-1)

      if( ia /= iaabsi .and. .not. State_all_r .and. Cal_xanes ) cycle

      if( Normaltau ) then

        do Lm1 = 1,nlmsa(ia)
          lmd1 = Lm0 + Lm1
          do Lm2 = 1,nlmsa(ia)
            lmd2 = Lm0 + Lm2
            Taullp(Lm2,Lm1) =  sum( mat(lmd2,Lm0+1:Lm0+nlmsa(ia)) * Taup(1:nlmsa(ia),Lm1,ia) )
          end do
        end do

      else

        do Lm1 = 1,nlmsa(ia)
          lmd1 = Lm0 + Lm1
          taullp(Lm1,Lm1) = mat(lmd1,lmd1)
! Because of complex harmonics, matrix is not symmetric
          do Lm2 = 1,nlmsa(ia)
            lmd2 = Lm0 + Lm2
            taullp(Lm1,Lm2) = mat(lmd1,lmd2)
          end do
        end do

      endif

! Setting of Tau in Xtal basis
      if( Spin_axis_free ) then

        allocate( Mat_AC(nlmsa(ia),nlmsa(ia)) )
        allocate( Mat_CA(nlmsa(ia),nlmsa(ia)) )
        Mat_AC(:,:) = ( 0._db, 0._db )
        Mat_CA(:,:) = ( 0._db, 0._db )
        do Lm = 1,nlmsa(ia)
          L = Lato(Lm,ia,igrph)
          do Lm_p = 1,nlmsa(ia)
            if( Lato(Lm_p,ia,igrph) /= L ) cycle 
            m_p = mato(Lm_p,ia,igrph)
            iso_p = iato(Lm_p,ia,igrph)
            Mat_AC(Lm,Lm_p) = DLmsLms_sym_AC(Lm,m_p,iso_p,1,ia) 
            Mat_CA(Lm_p,Lm) = DLmsLms_sym_CA(Lm,m_p,iso_p,1,ia) 
          end do
        end do
        allocate( Mat_tau(nlmsa(ia),nlmsa(ia)) )
        
        Mat_tau(1:nlmsa(ia),1:nlmsa(ia)) = Taullp(1:nlmsa(ia),1:nlmsa(ia))

        Taullp(:,:) = (0._db,0._db)
        
        Mat_tau = matmul( Mat_CA, matmul( Mat_tau, Mat_AC ) )    
        
        Taullp(1:nlmsa(ia),1:nlmsa(ia)) = Mat_tau(1:nlmsa(ia),1:nlmsa(ia)) 

        deallocate( Mat_AC, Mat_CA, Mat_tau )
        
      elseif( nb_sym_op > 1 ) then
      
        allocate( Cmatr(nlmsamax,-Lmax_g:Lmax_g,nspino) )

        Mat_rot(:,:) = Rot_atom(:,:,ia)
        Mat_rot = Transpose( Mat_rot )
        
        allocate( Dlmm_ia(-Lmaxa(ia):Lmaxa(ia),-Lmaxa(ia):Lmaxa(ia), 0:Lmaxa(ia)) )

        call Rotation_mat(Dlmm_ia,icheck,Mat_rot,Lmaxa(ia), Ylm_comp)
        do Lm = 1,nlmsa(ia)
          L = Lato(Lm,ia,igrph)
          do isp = 1,nspino
            do m = -L,L
              Cmatr(Lm,m,isp) = sum( Cmat(ia,Lm,1,-L:L,isp) * Dlmm_ia(-L:L,m,L) )
            end do
          end do
        end do

        deallocate( Dlmm_ia )

        taullq(:,:) = taullp(:,:)
        taullp(:,:) = (0._db,0._db)

        do Lm1 = 1,nlmsa(ia)
          L1 = Lato(Lm1,ia,igrph)
          m1 = mato(Lm1,ia,igrph)
          is1 = iato(Lm1,ia,igrph)
          do Lm2 = 1,nlmsa(ia)
            L2 = Lato(Lm2,ia,igrph)
            m2 = mato(Lm2,ia,igrph)
            is2 = iato(Lm2,ia,igrph)
      
            do lm1p = 1,nlmsa(ia)
              if( Lato(lm1p,ia,igrph) /= L1 ) cycle
              if( abs( Cmatr(lm1p,m1,is1) ) < eps10 ) cycle
              do lm2p = 1,nlmsa(ia)
                if( Lato(lm2p,ia,igrph) /= L2 ) cycle
                if( abs( Cmatr(lm2p,m2,is2) ) < eps10 ) cycle
      
                taullp(Lm1,Lm2) = taullp(Lm1,Lm2) + Cmatr(lm1p,m1,is1) * taullq(lm1p,lm2p) * Conjg( Cmatr(lm2p,m2,is2) )
      
              end do
            end do
      
          end do
        end do

        deallocate( Cmatr )
        
      endif

      do Lm1 = 1,nlmsa(ia)
        L1 = Lato(Lm1,ia,igrph)
        m1 = mato(Lm1,ia,igrph)
        lm01 = L1**2 + L1 + 1 + m1
        if( Spinorbite ) then
          is1 = iato(Lm1,ia,igrph)
        else
          is1 = ispin
        endif
        do Lm2 = 1,nlmsa(ia)
          L2 = Lato(Lm2,ia,igrph)
          m2 = mato(Lm2,ia,igrph)
          lm02 = L2**2 + L2 + 1 + m2
          if( Spinorbite ) then
            is2 = iato(Lm2,ia,igrph)
          else
            is2 = ispin
          endif
          Taull(lm01,is1,lm02,is2,ia) = Taull(lm01,is1,lm02,is2,ia) + Taullp(Lm1,Lm2)
 ! Even for the Atom out of the symmetry axis the contribution must be added corresponding to the conjugate representation 
          if( .not. Spinorbite .and. Repres_comp ) then
            lm01c = lm01 - 2 * m1
            Lm02c = lm02 - 2 * m2
            isg = (-1)**(m1+m2)
            Taull(Lm02c,is2,lm01c,is1,ia) = Taull(Lm02c,is2,lm01c,is1,ia) + isg * taullp(Lm1,Lm2)
          endif
        end do
      end do

    end do ! end of loop over ia

  endif

  if( natome > 1 .and. nab_coop > 0 ) then
    irep = abs( irep_util(igrph,1) )
    Sym_3 = abs(karact(49,irep)) > eps10 .or. abs(karact(53,irep)) > eps10

    call Cal_Tau_coop(Cmat,Dist_coop,ia_coop,iaprotoi,iato,igrph,ispin,Lato,Lmaxa,Lmax_g,mat,mato, &
                      n_atom_coop,n_atom_proto,nab_coop,natome,nb_sym_op,ndim,ngrph,nlmagm,nlmch,nlmsa,nlmsam,nlmsamax,nlmsmax, &
                      nspino,nspinp,Posi,Repres_comp,Rmtg,rot_atom,Spinorbite,Sym_3,Tau_coop,Ylm_comp)
  endif
 
  if( natome > 1 .or. nb_sym_op > 1 ) then
    if( Spin_axis_free ) then
      deallocate( DLmsLms_sym_AC, DLmsLms_sym_CA )
    else
      deallocate( Cmat )
    endif
  endif

  if( icheck > 2 ) then
    write(3,150)
    iligne = 0
    do ia = 1,natome
      do Lm = 1,nlmsa(ia)
        iligne = iligne + 1
        if( Spinorbite ) then
          write(3,130) iligne, ia, Lato(Lm,ia,igrph), mato(Lm,ia,igrph), iato(Lm,ia,igrph)
        else
          write(3,130) iligne, ia, Lato(Lm,ia,igrph), mato(Lm,ia,igrph), ispin
        endif
        write(3,140) mat(iligne,1:ndim)
      end do
    end do
  endif

  deallocate( mat )
  deallocate( taullp )
  deallocate( taullq )
  deallocate( taup )
  if( tau_nondiag ) deallocate( taup_inv )

  endif            ! arrival when ndim = 0

  if( recop .and. igrph == ngrph .and. ( Spinorbite .or. ispin == nspin ) ) then
    allocate( Taull_tem(nlmagm,nspinp,nlmagm,nspinp) )
    Taull_tem(:,:,:,:) = Taull(:,:,:,:,iaabsi)
    call recop_taull(nlmagm,nspinp,Sym_cubic,Taull_tem,Ylm_comp)
    Taull(:,:,:,:,iaabsi) = Taull_tem(:,:,:,:)
    deallocate( taull_tem )
  endif

! Substraction of the atomic scattering amplitude (for complex energy with irregular solution)
  if( Solsing .and. igrph == ngrph ) then
    if( icheck > 2 ) then
      write(3,161)
      call Write_Tau(iaabsi,iaabsi,Lmaxa,natome,nlmagm,nspinp,Taull)
    end if
    call Soustract_tl(Axe_Atom_gr,Cal_xanes,Full_atom,iaabsi,iaprotoi,igroupi,ispin,Lmaxa,n_atom_0,n_atom_ind, &
                natome,ngroup_m,nlmagm,nspin,nspino,nspinp,Rot_atom,Spin_axis_free,Spinorbite,State_all_r, &
                Tau_ato,Taull)
  end if

  if( icheck > 1 .and. igrph == ngrph .and. ( Spinorbite .or. ispin == nspin ) ) then
    do ia = 1,natome
      if(  ia /= iaabsi .and. .not. State_all_r .and. Cal_xanes ) cycle
      call Write_Tau(ia,iaabsi,Lmaxa,natome,nlmagm,nspinp,Taull)
    end do
  endif
  
  if( ndim == 0 ) then
    Time_tria = 0._db
  else
    call CPU_TIME(time)
    tp3 = real(time,db)
    Time_tria = tp3 - tp2
  endif

  return
  110 format(/' ---- Mat_MST ---------',100('-'))
  115 format(/ ' igrph =',i2,', Matrix dimension =',i7)
  120 format(/' Multiple scattering matrix'/)
  130 format(' iligne =',i3,', ia =',i2,', (L,m,isp) = (',3i2,')' )
  140 format(1p,12(1x,2e11.4))
  150 format(/' Inverse of the multiple scattering matrix'/)
  161 format(/' Scattering amplitude of the central atom (Taull) before eliminating the contribution of the singular solution:')
  218 format(3i3,2x,1p,34(1x,2e11.3))
  220 format(8x,a1,2x,102i3)
  230 format(3x,2i3,2x,102a3)
  240 format(3i3,2x,102a3)
  250 format(1p,5(1x,a3,' = ',2e15.7))    

end

!***********************************************************************
 
subroutine Write_Tau(ia,iaabsi,Lmaxa,natome,nlmagm,nspinp,Taull)

  use declarations
  implicit none

  integer, parameter:: nletm = 52
  
  integer:: i, ia, iaabsi, is1, is2, ival, L, L1, L2, letu, letd, lm1, lm2, lmw, m, m1, m2, natome, nlet, nlmagm, nlmw, nspinp
  
  integer, dimension(natome):: Lmaxa

  character(len=1) let(0:nletm)
  character(len=3), dimension(nletm**2):: motval
  character(len=3), dimension(nlmagm,nspinp,nlmagm,nspinp):: mot
  
  complex(kind=db), dimension(nletm**2):: val
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,natome):: Taull

  data let/' ','a','b','c','d','e','f','g','h','i','j','k','l','m', &
   'n','o','p','q','r','s','t','u','v','w','x','y','z','A','B','C', &
   'D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z'/

  if( ia == iaabsi ) then
    write(3,160)
  else
    write(3,162) ia
  endif
  lmw = min( 3, Lmaxa(ia) )
  nlmw = min( nlmagm, (lmw+1)**2 )

  val(:) = (0._db, 0._db)
  motval(:) = ' '
  mot(:,:,:,:) = ' '
  nlet = 0
  letu = 0 
  letd = 0
  do is1 = 1,nspinp
    lm1 = 0
    do L1 = 0,lmw
      do m1 = -L1,L1
        lm1 = lm1 + 1
        do is2 = 1,nspinp
          lm2 = 0
          do L2 = 0,lmw
            boucle_m2: do m2 = -L2,L2
              lm2 = lm2 + 1
              if( abs( Taull(lm1,is1,lm2,is2,ia) ) < eps15 ) then
                mot(lm1,is1,lm2,is2) = '  .'
              else
                do ival = 1,nlet
                  if( abs( Taull(lm1,is1,lm2,is2,ia) - val(ival) ) < eps10 ) then
                    mot(lm1,is1,lm2,is2) = motval(ival)
                    cycle boucle_m2
                  elseif( abs( Taull(lm1,is1,lm2,is2,ia) + val(ival) ) < eps10 ) then
                    mot(lm1,is1,lm2,is2) = motval(ival)
                    L = len_trim( adjustl( motval(ival) ) )
                    mot(lm1,is1,lm2,is2)(3-L:3-L) = '-'
                    cycle boucle_m2
                  endif
                end do

                if( nlet == nletm**2 ) then
                  mot(lm1,is1,lm2,is2) = '  ?'
                else
                  nlet = nlet + 1
                  letu = letu + 1
                  if( letu > nletm ) then
                    letu = 1
                    letd = letd + 1
                  endif
                  val(nlet) = Taull(lm1,is1,lm2,is2,ia)
                  motval(nlet)(2:3) = let(letd) // let(letu) 
                  mot(lm1,is1,lm2,is2) =  motval(nlet) 
                endif
              endif  
            end do boucle_m2
          end do
        end do
      end do
    end do
  end do

  write(3,220) 'l', ((( L, m = -L,L), L = 0,lmw ), i = 1,nspinp )
  write(3,220) 'm', ((( m, m = -L,L), L = 0,lmw ), i = 1,nspinp )
  if( nspinp == 2 ) then
    write(3,220) 's', ((( i, m = -L,L), L = 0,lmw ), i = 1,nspinp )
    write(3,'(A)') '  l  m  s'
  else
    write(3,'(A)') '     l  m'
  endif
  do is1 = 1,nspinp
    lm1 = 0
    do L1 = 0,lmw
      do m1 = -L1,L1
        lm1 = lm1 + 1
        if( nspinp == 1 ) then
          write(3,230) L1, m1, (( mot(lm1,is1,lm2,is2), lm2 = 1,nlmw), is2 = 1,nspinp )
        else
          write(3,240) L1, m1, is1, (( mot(lm1,is1,lm2,is2), lm2 = 1,nlmw), is2 = 1,nspinp )
        endif
      end do
    end do
  end do
  
  write(3,*)
  write(3,250) ( motval(i), val(i), i = 1,nlet )

  return
  160 format(/' Scattering amplitude of the central atom (Taull):')
  162 format(/' Scattering amplitude of the atom ia =',i3,' (Taull):')
  218 format(3i3,2x,1p,34(1x,2e11.3))
  220 format(8x,a1,2x,102i3)
  230 format(3x,2i3,2x,102a3)
  240 format(3i3,2x,102a3)
  250 format(1p,5(1x,a3,' = ',2e15.7))    
end

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

! Calculation of the matrix coefficients due to the projection to lower symmetry.

! Il y a peut-etre une erreur dans cette routine --> resultat different
! pour le cuivre cfc en symetrie 4/mmm et mmm (ou 1). A cause de ca, le
! groupe ponctuel est pris de symetrie plus basse en cas de green dans
! igrpt_sg_cal appele par la routine agregat.

subroutine Cmat_cal(Cmat,iato,icheck,igrph,iopsymr,irep_util,is_eq,karact,Lato,Lmax_g,mato,natome,nb_eq,nb_sym_op,ngrph, &
                nlmsamax,nlmsa,nlmsam,nspino,rot_atom,Ylm_comp)

  use declarations
  implicit none
  
  integer:: ia, icheck, igrph, ind, irep, is, isp, js, ks, L, Lm, lmax, Lmax_g, m, mp, natome, nb_sym_op, ngrph, nlmsamax, &
            nlmsam, nspino 

  complex(kind=db), dimension(nb_sym_op,nspino):: C_Gamma
  complex(kind=db), dimension(-Lmax_g:Lmax_g):: Ctem
  complex(kind=db), dimension(natome,nlmsamax,nb_sym_op,-Lmax_g:Lmax_g,nspino):: Cmat
  complex(kind=db), dimension(nopsm,nrepm):: karact
  complex(kind=db), dimension(:,:,:), allocatable:: Dlmm_ia, Dlmm_ib, Dlmm_is
  complex(kind=db), dimension(:,:), allocatable:: Mat_ia, Mat_is

  integer, dimension(natome):: nb_eq, nlmsa
  integer, dimension(nopsm):: iopsymr
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato
  integer, dimension(nb_sym_op,natome):: is_eq
  integer, dimension(nrepm,2):: irep_util

  logical:: Rotation, Sym_rot, Ylm_comp

  real(kind=db):: C_norm
  real(kind=db), dimension(3,3):: Mat_rot
  real(kind=db), dimension(3,3,natome):: rot_atom

  do isp = 1,nspino
    js = 0
    irep = abs( irep_util(igrph,isp) )
    if( irep == 0 ) cycle
    do is = 1,nopsm
      if( iopsymr(is) == 0 ) cycle
      js = js + 1
      if( irep_util(igrph,isp) > 0 ) then
        C_Gamma(js,isp) = Conjg( karact(is,irep) )
      else
        C_Gamma(js,isp) = karact(is,irep)
      endif
    end do
  end do
  Cmat(:,:,:,:,:) = (0._db, 0._db)

  do ia = 1,natome
    if( nlmsa(ia) == 0 ) cycle
    lmax = Lato(nlmsa(ia),ia,igrph)
    allocate( Dlmm_ia(-lmax:lmax,-lmax:lmax,0:lmax) )
    allocate( Dlmm_ib(-lmax:lmax,-lmax:lmax,0:lmax) )
    allocate( Dlmm_is(-lmax:lmax,-lmax:lmax,0:lmax) )

    if( abs( rot_atom(1,1,ia) - 1._db ) > eps10 .or. abs( rot_atom(2,2,ia) - 1._db ) > eps10 .or. &
        abs( rot_atom(3,3,ia) - 1._db ) > eps10 ) then
      rotation = .true.
    else
      rotation = .false.
    endif
    Mat_rot(:,:) = rot_atom(:,:,ia)
    call Rotation_mat(Dlmm_ia,icheck,Mat_rot,lmax,Ylm_comp)

    if( icheck > 3 ) then
      write(3,100) ia
      do L = 0,lmax
        write(3,102) L
        do m = -L,L
          write(3,104) Dlmm_ia(m,-L:L,L)
        end do
      end do
    endif

! Loop over equivalent atoms
    do ind = 1,nb_eq(ia)
      is = abs( is_eq(ind,ia) )

      js = 0
      do ks = 1,nopsm
        if( iopsymr(ks) == 0 ) cycle
        js = js + 1
        if( ks == is ) exit
      end do

      if( is == 1 .or. is == 25 ) then
        sym_rot = .false.
      else
        sym_rot = .true.
      endif

      call Dlmm_cal_is(icheck,is,lmax,Dlmm_is,Ylm_comp)

! La matrice de la symetrie est mise dans la base de l'atome ia.
! On prend tous les atomes equivalents selon la meme base.
      if( rotation .or. sym_rot ) then
        do L = 0,lmax
          allocate( Mat_is(-L:L,-L:L) )
          allocate( Mat_ia(-L:L,-L:L) )
          Mat_is(-L:L,-L:L) = Dlmm_is(-L:L,-L:L,L)
          Mat_ia(-L:L,-L:L) = Dlmm_ia(-L:L,-L:L,L)
          Mat_is = Matmul( Mat_is, Conjg(Transpose(Mat_ia)) )
          Mat_is = Matmul( Mat_ia, Mat_is )
          Dlmm_is(-L:L,-L:L,L) = Mat_is(-L:L,-L:L)
          deallocate( Mat_is )
          deallocate( Mat_ia )
        end do
      endif

      if( icheck > 3 ) then
        write(3,106) is, ind
        do L = 0,lmax
          write(3,102) L
          do m = -L,L
            write(3,104) Dlmm_is(m,-L:L,L)
          end do
        end do
      endif

      do Lm = 1,nlmsa(ia)
        L = Lato(Lm,ia,igrph)
        m = mato(Lm,ia,igrph)
        isp = iato(Lm,ia,igrph)
        do mp = -L,L
          Cmat(ia,Lm,ind,mp,isp) = Cmat(ia,Lm,ind,mp,isp) + C_Gamma(js,isp) * conjg(Dlmm_is(m,mp,L))
        end do
      end do
    end do

! Normalisation
    do Lm = 1,nlmsa(ia)
      do isp = 1,nspino
        C_norm = sum( abs( Cmat(ia,Lm,:,:,isp) )**2 )
        if( C_norm < eps10 ) cycle
        Cmat(ia,Lm,:,:,isp) = Cmat(ia,Lm,:,:,isp) / sqrt( C_norm )
      end do
    end do

    if( rotation ) then
      do ind = 1,nb_eq(ia)

! Supprime car tous les atomes 'ind' sont dans la base de 'ia'
!            isym = abs( is_eq(ind,ia) )
!            if( isym /= 1 ) then
!              call Dlmm_cal_is(icheck,isym,lmax,Dlmm_is,Ylm_comp)
!              do L = 0,lmax
!                allocate( Mat_is(-L:L,-L:L) )
!                allocate( Mat_ia(-L:L,-L:L) )
!                Mat_is(-L:L,-L:L) = Dlmm_is(-L:L,-L:L,L)
!                Mat_ia(-L:L,-L:L) = Dlmm_ia(-L:L,-L:L,L)
!                Mat_ia = Matmul( Mat_is, Mat_ia )
!                Dlmm_ib(-L:L,-L:L,L) = Mat_ia(-L:L,-L:L)
!                deallocate( Mat_is )
!                deallocate( Mat_ia )
!              end do
!            else
!              Dlmm_ib(:,:,:) = Dlmm_ia(:,:,:)
!            endif

! Rotation pour se mettre dans la base du cluster (et de la diffusion).
        do Lm = 1,nlmsa(ia)
          L = Lato(Lm,ia,igrph)
          do isp = 1,nspino
            Ctem(-L:L) = Cmat(ia,Lm,ind,-L:L,isp)
            do m = -L,L
              Cmat(ia,Lm,ind,m,isp) = sum( Ctem(-L:L) * Dlmm_ia(-L:L,m,L) )
            end do
          end do
        end do

      end do

    endif

    deallocate( Dlmm_ia )
    deallocate( Dlmm_ib )
    deallocate( Dlmm_is )

  end do

  if( icheck > 2 ) then
    write(3,110) igrph
    do ia = 1,natome
      write(3,*)
      do Lm = 1,nlmsa(ia)
          L = Lato(Lm,ia,igrph)
          do ind = 1,nb_eq(ia)
            do m = -L,L
              if( sum( abs(Cmat(ia,Lm,ind,m,:))) < eps10 ) cycle
              write(3,120) ia, Lm, ind, L, m, Cmat(ia,Lm,ind,m,:)
            end do
          end do
      end do
    end do
  endif

  return
  100 format(/' Dlmm_ia (Atom =',i2,') :')
!  100 format(/' Atom =',i2,', ind =',i2,', isym =',i3,//' Dlmm_ia :')
  102 format(/' l =',i2)
  104 format(7(1x,2f8.5))
  106 format(/' Dlmm_is ( is =',i3,', ind =',i2,' ) :')
  110 format(/' ia Lm  ind  l  m   Cmat = <rep,ia,Lm=(l,s)|ind,l,m>,    rep =',i2)
  120 format(2i3,2x,3i3,4f10.5)
end

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

! Calculation of the matrix coefficients due to the projection to lower symmetry.
! Case for spinors

subroutine DLmsLms_sym_cal(DLmsLms_sym_AC,DLmsLms_sym_CA,ia,iato,icheck,igrph,iopsymr,is_eq,Lato,Lmax,Lmax_g,mato,natome,nb_eq, &
                       nb_sym_op,ngrph,nlmsa,nlmsam,nlmsamax,nspino,Rot_atom)

  use declarations
  implicit none
  
  integer:: ia, icheck, igrph, ind, is, iso, iso_p, iso_pp, js, ks, L, Lm, Lm_p, Lmax, Lmax_g, m, m_p, m_pp, natome, &
            nb_sym_op, ngrph, nlm, nlmsam, nlmsamax, nspino 

  complex(kind=db), dimension(nlmsamax,-Lmax_g:Lmax_g,nspino,nb_sym_op,natome):: DLmsLms_sym_AC, DLmsLms_sym_CA
  complex(kind=db), dimension(:,:), allocatable:: Mat_ia, Mat_ia_i, Mat_is, Mat_is_i
  complex(kind=db), dimension(:,:,:,:,:), allocatable:: Dlmsms_ia, Dlmsms_ia_i, Dlmsms_ib, Dlmsms_is, Dlmsms_is_i
  complex(kind=db), dimension(:,:,:,:,:,:), allocatable:: DLmsLms_sym, DLmsLms_sym_i

  integer, dimension(natome):: nb_eq, nlmsa
  integer, dimension(nopsm):: iopsymr
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato
  integer, dimension(nb_sym_op,natome):: is_eq

  logical:: Rotation, Sym_rot

  real(kind=db):: C_norm
  real(kind=db), dimension(3,3):: Mat_rot, Matopsym
  real(kind=db), dimension(3,3,natome):: rot_atom

  if( nlmsa(ia) == 0 ) return
  
  allocate( DLmsLms_sym(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax,nb_sym_op) ) 
  allocate( DLmsLms_sym_i(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax,nb_sym_op) ) 
  DLmsLms_sym(:,:,:,:,:,:) = ( 0._db, 0._db )
  DLmsLms_sym_i(:,:,:,:,:,:) = ( 0._db, 0._db )
  
  if( abs( rot_atom(1,1,ia) - 1._db ) > eps10 .or. abs( rot_atom(2,2,ia) - 1._db ) > eps10 .or. &
      abs( rot_atom(3,3,ia) - 1._db ) > eps10 ) then
    Rotation = .true.
  else
    Rotation = .false.
  endif
  Mat_rot(:,:) = Rot_atom(:,:,ia)

  allocate( Dlmsms_ia(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax) )
  allocate( Dlmsms_ia_i(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax) )
  allocate( Dlmsms_ib(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax) )
  allocate( Dlmsms_is(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax) )
  allocate( Dlmsms_is_i(-Lmax:Lmax,nspino,-Lmax:Lmax,nspino,0:Lmax) )

  if( icheck > 2 ) write(3,'(/a6,i3,A)') ' Atom =', ia, ', Dlmsms_ia, Xtal to Atom basis'  
  call Dlmsms_half_integer(Dlmsms_ia,icheck,.false.,Lmax,Mat_rot) ! false = not inverse rotation = cluster basis towards Atom basis
  if( icheck > 2 ) write(3,'(/a6,i3,A)') ' Atom =', ia, ', Dlmsms_ia_i, Atom to Xtal basis'  
  call Dlmsms_half_integer(Dlmsms_ia_i,icheck,.true.,Lmax,Mat_rot) ! true = inverse rotation = Atom basis towards cluster basis

! Loop over equivalent atoms
  do ind = 1,nb_eq(ia)
    is = abs( is_eq(ind,ia) )

    js = 0
    do ks = 1,nopsm
      if( iopsymr(ks) == 0 ) cycle
      js = js + 1
      if( ks == is ) exit
    end do

!    if( is == 1 .or. is == 25 ) then
    if( is == 1 ) then
      sym_rot = .false.
    else
      sym_rot = .true.
    endif

    call opsym(is,Matopsym)
    if( icheck > 2 ) write(3,'(/a11,i3,A)') ' Symmetry =', is, ', Dlmsms_is'  
    call Dlmsms_half_integer(Dlmsms_is,icheck,.false.,Lmax,Matopsym) ! no inverse rotation
    if( is == 1 ) then
      Dlmsms_is_i(:,:,:,:,:) = Dlmsms_is(:,:,:,:,:)
    else
      if( icheck > 2 ) write(3,'(/a11,i3,A)') ' Symmetry =', is, ', Dlmsms_is, inverse'  
      call Dlmsms_half_integer(Dlmsms_is_i,icheck,.true.,Lmax,Matopsym) ! inverse rotation
    endif

! Symmetry matrix is set (at this step) in the atom basis
    if( Sym_rot ) then
      do L = 0,Lmax
        nlm = nspino * ( 2 * L + 1 )
        allocate( Mat_ia(nlm,nlm) )
        allocate( Mat_ia_i(nlm,nlm) )
        allocate( Mat_is(nlm,nlm) )
        allocate( Mat_is_i(nlm,nlm) )
        Lm = 0
        do iso = 1,nspino
          do m = - L,L
            Lm = Lm + 1
            Lm_p = 0
            do iso_p = 1,nspino
              do m_p = - L, L
                Lm_p = Lm_p + 1
                Mat_ia(Lm,Lm_p) = Dlmsms_ia(m,iso,m_p,iso_p,L)
                Mat_ia_i(Lm,Lm_p) = Dlmsms_ia_i(m,iso,m_p,iso_p,L)
                Mat_is(Lm,Lm_p) = Dlmsms_is(m,iso,m_p,iso_p,L)
                Mat_is_i(Lm,Lm_p) = Dlmsms_is_i(m,iso,m_p,iso_p,L)
              end do
            end do
          end do
        end do

        Mat_is = Matmul( Mat_ia_i, Matmul( Mat_is, Mat_ia ) )
        Mat_is_i = Matmul( Mat_ia_i, Matmul( Mat_is_i, Mat_ia ) )

        Lm = 0
        do iso = 1,nspino
          do m = - L,L
            Lm = Lm + 1
            Lm_p = 0
            do iso_p = 1,nspino
              do m_p = - L,L
                Lm_p = Lm_p + 1
                Dlmsms_is(m,iso,m_p,iso_p,L) = Mat_is(Lm,Lm_p)
                Dlmsms_is_i(m,iso,m_p,iso_p,L) = Mat_is_i(Lm,Lm_p)
              end do
            end do
          end do
        end do

        deallocate( Mat_ia, Mat_ia_i, Mat_is, Mat_is_i )
        
      end do
    endif

    DLmsLms_sym(:,:,:,:,:,ind) = Dlmsms_is(:,:,:,:,:)
    DLmsLms_sym_i(:,:,:,:,:,ind) = Dlmsms_is_i(:,:,:,:,:)
    
  end do ! end of loop over equivalent atoms (ind)

! Normalisation
  do L = 0,Lmax
    do m = -L,L
      do iso = 1,nspino
        C_norm = sum( abs( DLmsLms_sym(m,iso,-L:L,:,L,:) )**2 )
        if( C_norm > eps10 ) DLmsLms_sym(m,iso,-L:L,:,L,:) = DLmsLms_sym(m,iso,-L:L,:,L,:) / sqrt( C_norm )
        C_norm = sum( abs( DLmsLms_sym_i(m,iso,-L:L,:,L,:) )**2 )
        if( C_norm > eps10 ) DLmsLms_sym_i(m,iso,-L:L,:,L,:) = DLmsLms_sym_i(m,iso,-L:L,:,L,:) / sqrt( C_norm )
      end do
    end do
  end do

  do ind = 1,nb_eq(ia)
  
    do Lm = 1,nlmsa(ia)
      L = Lato(Lm,ia,igrph)
      m = mato(Lm,ia,igrph)
      iso = iato(Lm,ia,igrph)

! Rotation to go in the cluster basis
      if( Rotation ) then
  
        do m_p = -L,L
          do iso_p = 1,nspino
  
            do m_pp = -L,L
              do iso_pp = 1,nspino
            ! Cluster basis on the right side
                DLmsLms_sym_AC(Lm,m_p,iso_p,ind,ia) = DLmsLms_sym_AC(Lm,m_p,iso_p,ind,ia) &
                                                    + DLmsLms_sym(m,iso,m_pp,iso_pp,L,ind) * Dlmsms_ia_i(m_pp,iso_pp,m_p,iso_p,L)
            ! Cluster basis on the left side
                DLmsLms_sym_CA(Lm,m_p,iso_p,ind,ia) = DLmsLms_sym_CA(Lm,m_p,iso_p,ind,ia) &
                                                    + Dlmsms_ia(m_p,iso_p,m_pp,iso_pp,L) * DLmsLms_sym_i(m_pp,iso_pp,m,iso,L,ind)
              end do
            end do
          
          end do
        end do
      
      else
 
        do m_p = -L,L       
          do iso_p = 1,nspino
            DLmsLms_sym_AC(Lm,m_p,iso_p,ind,ia) = DLmsLms_sym(m,iso,m_p,iso_p,L,ind) 
            DLmsLms_sym_CA(Lm,m_p,iso_p,ind,ia) = DLmsLms_sym_i(m,iso,m_p,iso_p,L,ind)
          end do
        end do
      
      endif
        
    end do
  end do

  if( icheck > 2 ) then
    write(3,110) igrph
    write(3,*)
    do ind = 1,nb_eq(ia)
      do Lm = 1,nlmsa(ia)
        L = Lato(Lm,ia,igrph)
        m = mato(Lm,ia,igrph)
        iso = iato(Lm,ia,igrph)
        do m_p = -L,L
          do iso_p = 1,nspino 
            if( abs(DLmsLms_sym_AC(Lm,m_p,iso_p,ind,ia)) < eps10 ) cycle
            write(3,120) ia, ind, Lm, L, m, iso, m_p, iso_p, DLmsLms_sym_AC(Lm,m_p,iso_p,ind,ia), &
                                                             DLmsLms_sym_CA(Lm,m_p,iso_p,ind,ia)
          end do
        end do
      end do
    end do
  endif

  deallocate( DLmsLms_sym, DLmsLms_sym_i )
  
  return
  110 format(/' ia eq Lm    L  m is mp ip DLmsLms = <rep,ia,Lm=(l,m,s)|ind,L,m,s>,    Representation =',i2)
  120 format(3i3,2x,5i3,1p,2(1x,2e13.5))
end

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

subroutine Dlmm_cal_is(icheck,is,lmax,Dlmm_is,Ylm_comp)

  use declarations
  implicit none

  integer:: icheck, is, js, L, lmax, m, mp

  complex(kind=db), dimension(-lmax:lmax,-lmax:lmax,0:lmax):: Dlmm_is

  logical:: Ylm_comp

  real(kind=db):: Dlmm_cal
  real(kind=db), dimension(3,3):: matopsym

  select case(is)

    Case(1,22,23,24,25,40,41,42)

       do L = 0,lmax
         do m = -L,L
           do mp = -L,L
             Dlmm_is(m,mp,L) = cmplx( Dlmm_cal(is,L,m,mp), 0._db,db)
           end do
         end do
       end do

    Case default

      select case(is)

        case(26,27,28)
          js = is - 7
        case(29,30,31)
          js = is - 13
        case(32,33,34,35,36,37,38,39)
          js = is - 30
        case(43)
          js = 15
        case(44)
          js = 13
        case(45)
          js = 11
        case(46)
          js = 14
        case(47)
          js = 12
        case(48)
          js = 10
        case(53)
          js = 52
        case(54)
          js = 51
        case(55)
          js = 50
        case(56)
          js = 49
        case(57,59,61,63)
          js = is + 1
        case default
          js = is
      end select

      call opsym(js,matopsym)

      call Rotation_mat(Dlmm_is,icheck,matopsym,lmax,Ylm_comp)

! On ajoute l'inversion (changement de signe pour l impair).
      if( is /= js ) then
       do L = 1,lmax,2
         do m = -L,L
           do mp = -L,L
             Dlmm_is(m,mp,L) = - Dlmm_is(m,mp,L)
           end do
         end do
       end do
      endif

  end select

  return
end

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

function Dlmm_cal(isym,L,m,mp)

  use declarations
  implicit none
  
  integer:: im, isym, L, m, m2, mm, mp

  logical:: pair

  real(kind=db):: Dlmm_cal

  Dlmm_cal = 0._db

  if( m /= mp ) return

  pair = .false.
  im = abs( m )

  select case(isym)

    case(1)
      pair = .true.

    case(22)
      mm = mod(L+im,2)
      pair = (mm == 0 .and. m >= 0) .or. (mm == 1 .and. m < 0)

    case(23)
      mm = mod(L+2*im,2)
      pair = (mm == 0 .and. m >= 0) .or. (mm == 1 .and. m < 0)

    case(24)
      mm = mod(im,2)
      pair = mm == 0

    case(25)
      pair = mod(L,2) == 0

    case(40)
      m2 = mod(im,2)
      pair = ( m <= 0 .or. m2 == 0 ) .and. ( m >= 0 .or. m2 == 1 )

    case(41)
      pair = m >= 0

    case(42)
      pair = mod(L+m,2) == 0

  end select

  if( pair ) then
    Dlmm_cal = 1._db
  else
    Dlmm_cal = -1._db
  endif

  return
end

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

! Calculation of Wigner Dlmm versus cartesian rotation matrices.
! Comes from Didier Sebillaud
! Modified by Yves Joly in February 2008

subroutine Rotation_mat(Dlmm,icheck,Mat_rot,lmax,Ylm_comp)

  use declarations
  implicit none

  integer:: icheck, is, L, lmax, m, mp

  real(kind=db):: alfa, beta, gamma

  complex(kind=db), dimension(-lmax:lmax,-lmax:lmax):: Mat, Trans, Transi
  complex(kind=db), dimension(-lmax:lmax,-lmax:lmax,0:lmax):: Dlmm

  Logical:: Ylm_comp

  real(kind=db), dimension(3,3):: Mat_rot
  real(kind=db), dimension(-lmax:lmax,-lmax:lmax,0:lmax):: R

  call Euler_mat(icheck,Mat_rot,alfa,beta,gamma)

  call DJMN(beta,R,lmax)

  do L = 0,lmax
    do m = -L,L
      do mp = -L,L
        Dlmm(m,mp,L) = exp( - img * ( m*alfa + mp*gamma) ) * R(m,mp,L)
      end do
    end do
  end do

  if( .not. Ylm_comp ) then

! Case of real harmonics:
  
  Trans(:,:) = (0._db, 0._db)
  Transi(:,:) = (0._db, 0._db)
  do m = -lmax,lmax
    if( m > 0 ) then
      is = (-1)**m
      Trans(m,m) = cmplx( is * sqrt_1o2, 0._db, db)
      Trans(m,-m) = cmplx( 0._db, is * sqrt_1o2, db)
    elseif( m == 0 ) then
      Trans(0,0) = (1._db, 0._db)
    else
      Trans(m,-m) = cmplx(sqrt_1o2, 0._db, db)
      Trans(m,m) = cmplx(0._db, -sqrt_1o2, db)
    endif
  end do
  
  Transi = Conjg( Transpose( Trans ) )

  do L = 0,lmax

    do m = -L,L
      do mp = -L,L
        Mat(m,mp) = sum( Dlmm(m,-L:L,L) * Trans(-L:L,mp) )
      end do
    end do

    do m = -L,L
      do mp = -L,L
        Dlmm(m,mp,L) = sum( Transi(m,-L:L) * Mat(-L:L,mp) )
      end do
    end do

  end do
  
  endif

  if( icheck > 2 ) then
  
    write(3,100) alfa*180/pi, beta*180/pi, gamma*180/pi
    do L = 0,lmax
      Write(3,110) L
        write(3,120 ) (m, m =-L,L)
      do m = -L,L
        write(3,130 ) m, Dlmm(m,-L:L,L)
      end do
    end do
  endif

  return
  100 format(/' alfa, beta, gamma =',3f10.4)
  110 format(/' Dlmm,  L =',i3)
  120 format(5x,'m',100(13x,i3,11x))  
  130 format(i6,1p,100(1x,2e13.5))  
end

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

!  This routine calculates Wigner rotation matrices R^{L}_{M1 M2} up to
!             order LMAX, following Messiah's convention.
!                   They are stored as R(M2,M1,L).

Subroutine DJMN(beta,R,LMAX)

  use declarations
  implicit none

  integer, parameter:: Nfac = 125

  integer:: i, j, L, L1, Lmax, M, M1, M2, NL_M, EPS0

  real(kind=db):: beta, C, C1, CC, CF1, CF2, CL, CM, CM2, Cmul, Coef, D0, FLL1, FLL2, FLLM, ProdL, PRODM, Small, R_1, R_2, &
    R_A, R_B, S, S1, Sqr2 
  
  real(kind=db), dimension(-lmax:lmax,-lmax:lmax,0:lmax):: R
  real(kind=db), dimension(0:2*lmax,0:2*lmax,0:2*lmax):: CF
  real(kind=db), dimension(0:2*lmax,0:2*lmax):: expr
  real(kind=db), dimension(0:Nfac):: GLG

!  data SMALL,SQR2 /0.001_db,1.414213562_db/

  Small = 0.001_db
  Sqr2 = sqrt( 2._db )
  
!  Logarithm of the Gamma function : Log(N!)
  GLG(1) = 0._db
  do I = 2,Nfac
    J = I - 1
    GLG(I) = GLG(J) + Log( real(J, db) )
  end do

  NL_M = lmax + 1

!  Coefficients appearing in Wigner's rotation matrices

  expr(0,0) = 1._db
  do L = 1,2*NL_M-2
    do M1 = 0,L
      expr(M1,L) =  exp(0.5_db*(GLG(L+L+1)-GLG(L+M1+1)-GLG(L-M1+1)))
      if( L < NL_M ) then
        do M2 = 0,L
          CF(L,M1,M2) = sqrt( real( (L*L-M1*M1)*(L*L-M2*M2), db ) ) / L
        end do
      endif
    end do
  end do

  C = COS( beta ) * 0.5_db
  S = SIN( beta ) * 0.5_db
  CC = C + C
  CMUL = -1._db
  if( ABS(S) < SMALL ) then
    if( C > 0._db ) EPS0  =  1
    if( C < 0._db ) EPS0  =  -1
    do L = 0,LMAX
      do M1 = -L,L
        do M2 = -L,L
          if( M1 /=  M2 * EPS0 ) then
            R(M2,M1,L) = 0._db
          else
            if( EPS0 == 1 ) then
              R(M2,M1,L) = 1._db
            else
              if( MOD(L+M1,2) == 0 ) then
                R(M2,M1,L) = 1._db
              else
                R(M2,M1,L) = -1._db
              endif
            endif
          endif
        end do
      end do
    end do
  else
    S1 = S * SQR2
    C1 = 0.5_db + C
    R(0,0,0) = 1._db
    if( Lmax > 0 ) then
      R(-1,-1,1) = C1
      R(0,-1,1) = S1
      R(1,-1,1) = 1._db - C1
      R(-1,0,1) = -S1
      R(0,0,1) = CC
      R(1,0,1) = S1
      R(-1,1,1) = 1._db - C1 
      R(0,1,1) = -S1
      R(1,1,1) = C1
    endif

    PRODL = -S
    COEF = -S / C1
    CL = -1._db
    do L = 2,LMAX
      CL = -CL
      L1 = L - 1
      FLL1 = CC * ( L + L1 )
      FLL2 = 1._db / ( L * L1 * CC )
      PRODL = - PRODL * S

!  Case M = 0

      R_1 = expr(0,L) * PRODL
      R(-L,0,L) = R_1

      R(L,0,L) = R_1 * CL
      R(0,-L,L) = R_1 * CL

      R(0,L,L) = R_1

      CM2 = CL
      do M2 = -L1,-1
        CM2 = CM2 * CMUL
        CF1 = CF(L1,0,-M2) / FLL1
        CF2 = FLL1 / CF(L,0,-M2)
        if( -M2 < L1 ) then
          R_A = CF2 * ( R(M2,0,L1) - R(M2,0,L-2) * CF1 )
        else
          R_A = CF2 * R(M2,0,L1)
        endif

        R(M2,0,L) = R_A

        R(-M2,0,L) = R_A * CM2
        R(0,M2,L) = R_A * CM2

        R(0,-M2,L) = R_A

      end do

      R(0,0,L) = FLL1 * R(0,0,L1) / CF(L,0,0) - R(0,0,L-2) * CF(L1,0,0) / CF(L,0,0)

!  Case M > 0

      PRODM = 1._db
      CM = CL
      FLLM = 0._db
      do M = 1,L1
        CM = -CM
        PRODM = PRODM * COEF
        FLLM = FLLM + FLL2

        R_1 = expr(M,L) * PRODL * PRODM
        R_2 = R_1 / ( PRODM * PRODM )

        R(-L,M,L) = R_1
        R(-L,-M,L) = R_2

        R(L,-M,L) = R_1 * CM
        R(M,-L,L) = R_1 * CM
        R(L,M,L) = R_2 * CM
        R(-M,-L,L) = R_2 * CM

        R(-M,L,L) = R_1
        R(M,L,L) = R_2

        CM2 = CM
        do M2 = -L1,-M
          CM2 = -CM2
          D0 = M2 * FLLM
          CF1 = CF(L1,M,-M2) / FLL1
          CF2 = FLL1 / CF(L,M,-M2)
          if( (M < L1) .AND. (-M2 < L1) ) then
            R_A = CF2 * ( (1._db-D0)*R(M2,M,L1) - R(M2,M,L-2)*CF1 )
            R_B = CF2 * ( (1._db+D0)*R(M2,-M,L1) - R(M2,-M,L-2)*CF1)
          else
            R_A = CF2 * ( 1._db - D0 ) * R(M2,M,L1)
            R_B = CF2 * ( 1._db + D0 ) * R(M2,-M,L1)
          endif

          R(M2,M,L) = R_A
          R(M2,-M,L) = R_B

          R(-M2,-M,L) = R_A * CM2
          R(M,M2,L) = R_A * CM2
          R(-M,M2,L) = R_B * CM2
          R(-M2,M,L) = R_B * CM2

          R(-M,-M2,L) = R_A
          R(M,-M2,L) = R_B

        end do
      end do

      PRODM = PRODM * COEF
      R_1 = PRODL * PRODM
      R_2 = PRODL / PRODM
      R(-L,L,L) = R_1
      R(L,-L,L) = R_1
      R(L,L,L) = R_2
      R(-L,-L,L) = R_2

    end do
  endif

  return
end

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

subroutine Soustract_tl(Axe_Atom_gr,Cal_xanes,Full_atom,iaabsi,iaprotoi,igroupi,ispinin,Lmaxa,n_atom_0,n_atom_ind, &
                natome,ngroup_m,nlmagm,nspin,nspino,nspinp,Rot_atom,Spin_axis_free,Spinorbite,State_all_r, &
                Tau_ato,Taull)

  use declarations
  implicit none

  integer:: i1, i2, ia, iaabsi, iapr, icheck, iga, igr, is1, is2, isp1, isp2, ispinin, its, L1, L2, &
     lm10, lm20, lm1, lm2, m1, m2, n_atom_0, n_atom_ind, natome, ngroup_m, nlm, nlmagm, nspin, nspino, nspinp

  integer, dimension(natome):: iaprotoi, igroupi, Lmaxa

  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,n_atom_0:n_atom_ind):: tau_ato
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,natome):: Taull

  complex(kind=db), dimension(:,:,:,:,:), allocatable:: Taull_p

  logical:: Cal_xanes, Full_atom, Spin_axis_free, Spinorbite, State_all_r

  real(kind=db):: cosang
  real(kind=db), dimension(3,ngroup_m):: Axe_atom_gr
  real(kind=db), dimension(3,3,natome):: Rot_atom

  icheck = 0

  if( Spin_axis_free ) then
    allocate( Taull_p(nlmagm,nspinp,nlmagm,nspinp,natome) )
    Taull_p(:,:,:,:,:) = (0._db, 0._db)
    do ia = 1,natome
      nlm = ( Lmaxa(ia) + 1 )**2
      if( Full_atom ) then
        iapr = ia 
      else
        iapr = iaprotoi( ia )
      endif
      Taull_p(1:nlm,:,1:nlm,:,ia) = Tau_ato(1:nlm,:,1:nlm,:,iapr) 
    end do

    call Rotate_TauLL(icheck,Lmaxa,natome,nlmagm,nspinp,Rot_atom,TauLL_p,.false.)
  endif

  do ia = 1,natome

    if( ia /= iaabsi .and. .not. State_all_r .and. Cal_xanes ) cycle

    if( nspin == 2 .and. .not. Spin_axis_free ) then
      igr = igroupi(ia)
      iga = igroupi(iaabsi)

      cosang = sum( Axe_Atom_gr(:,iga) * Axe_Atom_gr(:,igr) )

      if( abs( cosang - 1 ) < eps4 ) then
        its = 1
      elseif( abs( cosang + 1 ) < eps4 ) then
        its = -1
      else
        its = 0
      endif
    else
      its = 1
    endif

    if( Full_atom ) then
      iapr = ia
    else
      iapr = iaprotoi(ia)
    endif

    lm1 = 0
    do L1 = 0,Lmaxa(ia)
      do m1 = -L1,L1
        lm1 = lm1 + 1
        do is1 = 1,nspino
          if( Spinorbite ) then
            isp1 = is1
          else
            isp1 = ispinin
          endif
          if( its == - 1 ) then
            i1 = 3 - isp1
            lm10 = L1**2 + L1 + 1 - m1
          else
            i1 = isp1
            lm10 = L1**2 + L1 + 1 + m1
          endif
          lm2 = 0
          do L2 = 0,Lmaxa(ia)
! Even in Hubbard, different L do not mix
            do m2 = -L2,L2
              lm2 = lm2 + 1
              if( L2 /= L1 ) cycle
              do is2 = 1,nspino
!                    if( is1 /= is2 ) cycle
                if( Spinorbite ) then
                  isp2 = is2
                else
                  isp2 = ispinin
                endif
                if( its == - 1 ) then
                  i2 = 3 - isp2
                  lm20 = L2**2 + L2 + 1 - m2
                else
                  i2 = isp2
                  lm20 = L2**2 + L2 + 1 + m2
                endif
! if( Lm10 /= Lm20 ) cycle               
                if( Spin_axis_free ) then
                  Taull(lm1,isp1,lm2,isp2,ia) = Taull(lm1,isp1,lm2,isp2,ia) - Taull_p(lm1,isp1,lm2,isp2,ia)
                elseif( its == 0 ) then
                  if( i1 == i2 ) Taull(lm1,isp1,lm2,isp2,ia) = Taull(lm1,isp1,lm2,isp2,ia) - 0.5_db &
                     * ( Tau_ato(lm10,1,lm20,1,iapr) + Tau_ato(lm10,nspinp,lm20,nspinp,iapr) )
                else
                  Taull(lm1,isp1,lm2,isp2,ia) = Taull(lm1,isp1,lm2,isp2,ia) - Tau_ato(lm10,i1,lm20,i2,iapr)
                end if
              end do
            end do
          end do
        end do
      end do
    end do

  end do
  
  if( Spin_axis_free ) deallocate( Taull_p )

  return
end

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

subroutine Cal_Tau_coop(Cmat,Dist_coop,ia_coop,iaprotoi,iato,igrph,ispin,Lato,Lmaxa,Lmax_g,mat,mato, &
                      n_atom_coop,n_atom_proto,nab_coop,natome,nb_sym_op,ndim,ngrph,nlmagm,nlmch,nlmsa,nlmsam,nlmsamax,nlmsmax, &
                      nspino,nspinp,Posi,Repres_comp,Rmtg,rot_atom,Spinorbite,Sym_3,Tau_coop,Ylm_comp)

  use declarations
  implicit none
  
  integer:: i_trans, ia, ia_tr, iab, ib, ib_tr, ic, igrph, ipra, iprb, is1, is2, ispin, isg, isp, j, L, L1, L2, Lm, lm01, lm01c, &
    lm02, Lm02c, lm1, lm1p, lm2, lm2p, lma, lma0, lma0_tr, Lmax_g, lmb, lmb0, lmb0_tr, m, m1, m2, n_atom_coop, n_atom_proto, &
    natome, nb_sym_op, ndim, ngrph, nlmch, nlmagm, nlmsam, nlmsamax, nlmsmax, nspino, nspinp, nab_coop

  integer, dimension(n_atom_coop):: ia_coop
  integer, dimension(natome):: iaprotoi, Lmaxa, nlmsa
  integer, dimension(nlmsam,natome,ngrph):: Lato, mato, iato
  
  logical:: Repres_comp, Spinorbite, Sym_3, Ylm_comp
  
  real(kind=db):: Dist
  real(kind=db), dimension(2):: Dist_coop
  real(kind=db), dimension(0:n_atom_proto):: Rmtg
  real(kind=db), dimension(3,3):: Mat_rot
  real(kind=db), dimension(3,natome):: Posi
  real(kind=db), dimension(3,3,natome):: rot_atom

  complex(kind=db), dimension(ndim,nlmch):: mat
  complex(kind=db), dimension(nlmsmax,nlmsmax):: taullp, taullq
  complex(kind=db), dimension(natome,nlmsamax,nb_sym_op,-Lmax_g:Lmax_g,nspino):: Cmat
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,nab_coop,nspino):: Tau_coop
  complex(kind=db), dimension(:,:,:), allocatable:: Dlmm_ia
  complex(kind=db), dimension(:,:,:,:), allocatable:: Cmatr

  iab = 0

  do ia = 1,natome

    ipra = iaprotoi(ia)
    lma0 = sum( nlmsa(1:ia-1) )

    lmb0 = 0
    do ib = ia+1,natome

      if( n_atom_coop > 0 ) then
        do j = 1,n_atom_coop
          if( ia_coop(j) == 0 ) cycle
          if( ia == ia_coop(j) .or. ib == ia_coop(j) ) exit
        end do
        if( j > n_atom_coop ) cycle
      endif

      iprb = iaprotoi(ib)
      Dist = sqrt( sum( ( Posi(:,ia) - Posi(:,ib) )**2 ) )
      if( Dist < Dist_coop(1) - eps10 ) cycle
      if( Dist_coop(2) > eps10 ) then
        if( Dist > Dist_coop(2) + eps10 ) cycle
      else
        if( Dist > Rmtg( ipra ) + Rmtg( iprb ) ) cycle
      endif

      iab = iab + 1

      lmb0 = sum( nlmsa(1:ib-1) )

      if( nb_sym_op > 1 ) then

        allocate( Cmatr(nlmsamax,-Lmax_g:Lmax_g,nspino,ia:ib) )
    
        do ic = ia,ib,ib-ia
          Mat_rot(:,:) = rot_atom(:,:,ic)
          Mat_rot = Transpose( Mat_rot )
          allocate( Dlmm_ia(-Lmaxa(ic):Lmaxa(ic),-Lmaxa(ic):Lmaxa(ic),0:Lmaxa(ic)) )
          call Rotation_mat(Dlmm_ia,0,Mat_rot,Lmaxa(ic),Ylm_comp)
          do Lm = 1,nlmsa(ic)
            L = Lato(Lm,ic,igrph)
            do isp = 1,nspino
              do m = -L,L
                Cmatr(Lm,m,isp,ic) = sum( Cmat(ic,Lm,1,-L:L,isp) * Dlmm_ia(-L:L,m,L) )
              end do
            end do
          end do
          deallocate( Dlmm_ia )
        end do
    
      endif

      do i_trans = 1,nspino ! loop over transpose (usefull when spinorbit)

        if( i_trans == 1 ) then
          ia_tr = ia
          ib_tr = ib
          lma0_tr = lma0
          lmb0_tr = lmb0
        else
          ia_tr = ib
          ib_tr = ia
          lma0_tr = lmb0
          lmb0_tr = lma0
        endif
              
        do lma = 1,nlmsa(ia_tr)
          do lmb = 1,nlmsa(ib_tr)
            taullp(lma,lmb) = mat(lma0_tr+lma,lmb0_tr+lmb)
          end do
        end do
    
        if( nb_sym_op > 1 ) then
    
          taullq(:,:) = taullp(:,:)
          taullp(:,:) = (0._db,0._db)
    
          do lm1 = 1,nlmsa(ia_tr)
            L1 = Lato(lm1,ia_tr,igrph)
            m1 = mato(lm1,ia_tr,igrph)
            is1 = iato(lm1,ia_tr,igrph)
            do lm2 = 1,nlmsa(ib_tr)
              L2 = Lato(lm2,ib_tr,igrph)
              m2 = mato(lm2,ib_tr,igrph)
              is2 = iato(lm2,ib_tr,igrph)
    
              do lm1p = 1,nlmsa(ia_tr)
                if( Lato(lm1p,ia_tr,igrph) /= L1 ) cycle
                if( abs( Cmatr(lm1p,m1,is1,ia_tr) ) < eps10 ) cycle
                do lm2p = 1,nlmsa(ib_tr)
                  if( Lato(lm2p,ib_tr,igrph) /= L2 ) cycle
                  if( abs( Cmatr(lm2p,m2,is2,ib_tr) ) < eps10 ) cycle
    
                  Taullp(lm1,lm2) = Taullp(lm1,lm2) + Cmatr(lm1p,m1,is1,ia_tr) * Taullq(lm1p,lm2p) &
                                                                               * Conjg( Cmatr(lm2p,m2,is2,ib_tr) )
    
                end do
              end do
    
            end do
          end do
    
          if( i_trans == nspino ) deallocate( Cmatr )
    
        endif
    
        do lm1 = 1,nlmsa(ia_tr)
          L = Lato(lm1,ia_tr,igrph)
          m1 = mato(lm1,ia_tr,igrph)
          lm01 = L**2 + L + 1 + m1
          if( Spinorbite ) then
            is1 = iato(lm1,ia_tr,igrph)
          else
            is1 = ispin
          endif
          do lm2 = 1,nlmsa(ib_tr)
            L = Lato(lm2,ib_tr,igrph)
            m2 = mato(lm2,ib_tr,igrph)
            lm02 = L**2 + L + 1 + m2
            if( Spinorbite ) then
              is2 = iato(lm2,ib_tr,igrph)
            else
              is2 = ispin
            endif
            Tau_coop(lm01,is1,lm02,is2,iab,i_trans) = Tau_coop(lm01,is1,lm02,is2,iab,i_trans) + taullp(lm1,lm2)
 
            if( .not. Spinorbite .and. Repres_comp ) then
              lm01c = lm01 - 2 * m1
              Lm02c = lm02 - 2 * m2
              isg = (-1)**(m1+m2)
!              Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) = Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) + isg * taullp(lm1,lm2)
              if( i_trans == 1 ) then
                if( Sym_3 ) then
! Empirical finding of the power of img in order Sym 3 correct
                  Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) = Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) &
                                                            + isg * taullp(lm1,lm2)
                else  
! Empirical finding of the power of img in order sym 2/m with Ycomp and 4/m give same result
                  Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) = Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) &
                                                            + isg * img**(-m2+1) * taullp(lm1,lm2)
                endif

              else
           ! Not used because only with spinorbite n_trans can be = 2
                Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) = Tau_coop(lm01c,is1,Lm02c,is2,iab,i_trans) &
                                                          + isg * img**(m1-1) * taullp(lm1,lm2)
              endif
            endif
          end do
        end do

      end do ! end of loop over i_trans
      
    end do
  end do ! end of loop over ia

  return
end

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

! Writing or lecture of the scattering amplitudes when calculation is done using One run.

subroutine Data_one_run(First_run,i_layer,ia_eq,iabsm,icheck,igr_layer,igreq,index_e,igroup,igroupi,ipr1,lmax_probe, &
              Lmaxa,multi_run,n_abs_bulk,n_atom_proto,n_layer_abs,n_multi_run,n_proto_layer,n_proto_layer_max,natome,natomp,nb_eq, &
              nb_sym_op,neqm,nge,ngreq,nlm_probe,nlmagm,nspinp,One_run_s,Rot_atom,Rot_cluster,Spinorbite,Taull,Taull_stk,Ylm_comp)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: i, i_layer, ia, ib, icheck, igr, index_e, ipr, ipr1, is1, is2, j, L, L1, L2, lm1, lm2, lmax_probe, lms1, &
    lms2, lmx, lp1, lp2, m, m1, m2, mu, mv, multi_run, n_abs_bulk, n_atom_proto, &
    n_layer_abs, n_multi_run, n_proto_layer_max, natome, natomp, nb_sym_op, neqm, nge, nlm_probe, nlmagm, nlmw, nspinp
            
  integer, dimension(natome):: igroupi, Lmaxa, nb_eq
  integer, dimension(natomp):: igroup
  integer, dimension(n_multi_run):: iabsm
  integer, dimension(0:n_atom_proto):: ngreq
  integer, dimension(0:n_atom_proto,neqm):: igreq
  integer, dimension(n_layer_abs):: n_proto_layer
  integer, dimension(n_proto_layer_max,n_layer_abs):: igr_layer
  integer, dimension(nb_sym_op,natome):: ia_eq

  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp,n_multi_run-n_abs_bulk-1,nge):: Taull_stk

  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,natome):: Taull
  complex(kind=db), dimension(:,:), allocatable:: Mat, Mat_1, Mat_2
  complex(kind=db), dimension(:,:,:), allocatable:: Dlmm

  logical:: First_run, Found, One_run_s, Spinorbite, Ylm_comp

  real(kind=db), dimension(3,3):: Mat_rot, Rot_cluster
  real(kind=db), dimension(3,3,natome):: Rot_atom

!  if( index_e == 1 .and. mpinodes == 1 ) then
!    if( multi_run == 1 ) then
!      Open(50, status = 'SCRATCH')
!    else
!      Rewind(50)
!    endif
!  endif

  mv = 0
  do mu = 1,n_multi_run
    if( n_abs_bulk == 0 .and. mu == 1 ) cycle
    if( n_abs_bulk /= 0 .and. mu >= n_multi_run - n_abs_bulk ) exit
    mv = mv + 1
    if( One_run_s ) then
      do i = 1,n_proto_layer(i_layer)
        if( iabsm(mu) == igr_layer(i,i_layer) ) exit
      end do
      if( i > n_proto_layer(i_layer) ) cycle  
    endif

    Found = .false.

    do ipr = ipr1,n_atom_proto
      if( igreq(ipr,1) == iabsm(mu) ) exit
    end do

    do igr = 1,ngreq(ipr)
      Loop_ia: do ia = 1,natome
        if( igreq(ipr,igr) == igroupi(ia) ) then
          Found = .true.
          exit Loop_ia
        endif
        do i = 1,nb_eq(ia)
          ib = ia_eq(i,ia)
          if( igreq(ipr,igr) == igroup(ib) ) then
            Found = .true.
            exit Loop_ia
          endif
        end do
      end do Loop_ia
      if( .not. Found ) cycle

      lmx = min(Lmaxa(ia), lmax_probe)

      if( abs( rot_atom(1,1,ia) - 1._db ) > eps10 .or. abs( rot_atom(2,2,ia) - 1._db ) > eps10 .or. &
        abs( rot_atom(3,3,ia) - 1._db ) > eps10 ) then
        Mat_rot(:,:) = rot_atom(:,:,ia)
        Mat_rot = Matmul( Rot_cluster, Mat_rot )
      else
        Mat_rot = Rot_cluster
      endif

! When writing, one applies the inverse rotation --> Transpose
      if( First_run ) Mat_rot = Transpose( Mat_rot )

      Allocate( Dlmm(-lmx:lmx,-lmx:lmx,0:lmx) )
      call Rotation_mat(Dlmm,icheck,Mat_rot,lmx,Ylm_comp)

      do is1 = 1,nspinp
        do is2 = 1,nspinp
          if( .not. spinorbite .and. is1 /= is2 ) cycle

          do L1 = 0,lmx
            lp1 = L1**2 + L1 + 1
            Allocate( Mat_1(-L1:L1,-L1:L1) )
            Mat_1(-L1:L1,-L1:L1) = Dlmm(-L1:L1,-L1:L1,L1)
            Mat_1 = Transpose( Mat_1 )

            do L2 = 0,lmx
              lp2 = L2**2 + L2 + 1
              Allocate( Mat_2(-L2:L2,-L2:L2) )
              Mat_2(-L2:L2,-L2:L2) = Dlmm(-L2:L2,-L2:L2,L2)
              Allocate( Mat(-L1:L1,-L2:L2) )

              lm1 = lp1
              lms1 = lp1
              lm2 = lp2
              lms2 = lp2

              if( First_run ) then
            
                do m1 = -L1,L1
                  do m2 = -L2,L2
                    Mat(m1,m2) = Taull(lm1+m1,is1,lm2+m2,is2,ia)
                  end do
                end do
                Mat = Matmul( Mat_1, Matmul( Mat, Mat_2) )
!                if( mpinodes == 1 ) then
!                  write(50,*) Mat(:,:)
!                else
                  do m1 = -L1,L1
                    do m2 = -L2,L2
                      Taull_stk(lms1+m1,is1,lms2+m2,is2,mv,index_e) = Mat(m1,m2)
                    end do
                  end do
!                endif
            
              else
            
!                if( mpinodes == 1 ) then
!                  Read(50,*) Mat(:,:)
!                else
                  do m1 = -L1,L1
                    do m2 = -L2,L2
                      Mat(m1,m2) = Taull_stk(lms1+m1,is1,lms2+m2,is2,mv,index_e)
                    end do
                  end do
!                endif
                Mat = Matmul( Mat_1, Matmul( Mat, Mat_2) )
                do m1 = -L1,L1
                  do m2 = -L2,L2
                    Taull(lm1+m1,is1,lm2+m2,is2,ia) = Mat(m1,m2)
                  end do
                end do
            
              endif
              
              Deallocate( Mat, Mat_2 )
            end do

            Deallocate( Mat_1 )
          end do

        end do
      end do

      deallocate( Dlmm )

      exit
    end do

    if( .not. Found ) then
      call write_error
      do ipr = 6,9,3
        write(ipr,100) mu
      end do
      stop
    endif

  end do

  if( icheck > 1 .and. .not. First_run ) then

    do ia = 1,natome
      do i = 1,nb_eq(ia)
        ib = ia_eq(i,ia)
        if( igroup(ib) /= iabsm(multi_run) ) cycle

        nlmw = min(nlmagm,9)
        write(3,110)
        write(3,120)
        if( nspinp == 2 ) then
          write(3,210) ((( L, m, j, j = 1,2 ), m = -L,L), L = 0,2 )
        else
          write(3,220) (( L, m, m = -L,L), L = 0,2 )
        endif

        lm1 = 0
        do L1 = 0,lmx
          do m1 = -L1,L1
            lm1 = lm1 + 1
            do is1 = 1,nspinp
              write(3,230) L1, m1, is1, (( Taull(lm1,is1,lm2,is2,ia), is2 = 1,nspinp ), lm2 = 1,nlmw)
            end do
          end do
        end do
        exit

      end do
    end do
    
  endif

  return
  100 format(//' Prototypical atom number ',i3,' not in the sphere of calculation !'//&
               ' Increase the cluster radius !'/)
  110 format(/' ---- Data_One_Run ',100('-'))
  120 format(/' Taull: Multiple scattering amplitude:')
  210 format('( l, m, s)',50(7x,3i3,6x))
  220 format('( l, m, s)',50(10x,2i3,7x))
  230 format(3i3,2x,1p,50(1x,2e11.3))
end

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

! Reading of Taull when Extract_tau

subroutine E_reading(icheck,Eimag,Energ,ie,multi_run,nom_fich_extract)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: i, icheck, ie, istat, multi_run
  
  character(len=32):: mot32
  character(len=Length_name):: nom_fich_extract
 
  real(kind=db):: Energ, Eimag

  if( ie == 1 ) then
  
    open(1, file = nom_fich_extract, status='old', iostat=istat)
    if( istat /= 0 ) call write_open_error(nom_fich_extract,istat,1)

    i = 0
    do
     read(1,'(A)' ) mot32
      if( mot32(2:15) /= 'Absorbing atom' ) cycle
      i = i + 1
      if( i == multi_run ) exit
    end do

  endif

  do
    read(1,'(A)' ) mot32
    if( mot32 == '        Energ              Eimag' ) exit
  end do
  
  read(1,*) Energ, Eimag
  
  if( icheck > 1 ) then
    write(3,110)
    write(3,260)
    write(3,270) Energ, Eimag
  endif

  Eimag = Eimag / Rydb
  Energ = Energ / Rydb
  
  return
  110 format(/' ---- Tau_reading ---------',100('-'))
  260 format(/'        Energy          Eimag')
  270 format(1p,2e19.11)
end

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

! Writing for extract

subroutine Tau_writing(Eimag,Energ,iprabs,ie_computer,itypepr,lmaxat,lmax_probe,mpinodes,n_atom_proto,nlmmax, &
                       nspinp,ntype,Numat,Taull_dft)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: ie_computer, ipr, iprabs, isp, jsp, L, Lm, lmax, lmax_probe, m, mpinodes, n_atom_proto, nlm, nlmmax, nspinp, ntype
  integer, dimension(0:ntype):: Numat
  integer, dimension(0:n_atom_proto):: itypepr
  integer, dimension(0:n_atom_proto,0:mpinodes-1):: lmaxat
 
  complex(kind=db), dimension(nlmmax,nspinp,nlmmax,nspinp,0:mpinodes-1):: Taull_dft

  real(kind=db):: Energ, Eimag
   
  write(3,110)
  write(3,'(/A)') '        Energ              Eimag'   
  write(3,'(1p,2e19.11)') Energ*rydb, Eimag*rydb
  write(3,'(/A)') ' Z, lmax  (1,... n_atom_proto)'   
  write(3,'(10000(i3,i2))') ( Numat(itypepr(ipr)), lmaxat(ipr,ie_computer), ipr = 1,n_atom_proto )

! To save space in the bav file one limits the printing down to 3.
! Tau_dft is just used to calculate cross section after extraction...
  lmax = max( lmax_probe, 3 )
  lmax = min( lmaxat(iprabs,ie_computer), lmax )
  nlm = ( lmax + 1 )**2

  write(3,260)
  if( nspinp == 1 ) then
    write(3,270) lmax, nspinp
    write(3,280) (( L, m, m = -L,L ), L = 0,lmax )
  else
    write(3,290) lmax, nspinp
    write(3,300) ((( L, m, isp, m = -L,L ), L = 0,lmax ), isp = 1,nspinp )
  endif
  do isp = 1,nspinp
    do Lm = 1,nlm
      write(3,310) ( Taull_dft(Lm,isp,1:nlm,jsp,ie_computer), jsp = 1,nspinp )
    end do
  end do

  return
  110 format(/' ---- Tau ---------',100('-'))
  260 format(/' Multiple scattering amplitude')
  270 format(2i4,' = lmax, nspinp / ( l, m)')
  280 format(2420(17x,2i3,16x))
  290 format(2i4,' = lmax, nspinp / ( l, m, s)')
  300 format(2420(16x,3i3,14x))
  310 format(1p,2420(1x,2e19.11))

end

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

subroutine Tau_reading(icheck,ie,iaabsi,Lmaxa,natome,nenerg,nlmagm,nspinp,RIXS,Taull)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: i, iaabsi, icheck, ie, isp, jsp, L, Lm, lmax, Lmaxa, Lmp, m, natome, nenerg, nlm, nlma, nlmagm, nspinp

  character(len=20):: mot20
    
  complex(kind=db), dimension(nlmagm,nspinp,nlmagm,nspinp,natome):: Taull

  logical:: RIXS
  
  real(kind=db):: p
  real(kind=db), dimension(nlmagm,nspinp):: Taulli, Taullr

  do
    read(1,'(A)' ) mot20
    if( mot20 == ' Multiple scattering' ) exit
  end do

  read(1,*) lmax, nspinp

! Lmaxa can be < lmax because lmax is defined by the max of the lmax of the atoms with same atomic numer than the absorbing one 
  nlm = ( lmax + 1 )**2
  if( Lmaxa < lmax ) then
    nlma = ( Lmaxa + 1 )**2
  else
    nlma = nlm
  endif

  Taull(:,:,:,:,:) = ( 0._db, 0._db)

  read(1,*)
  
  do isp = 1,nspinp
    do Lm = 1,nlma
      read(1,*) ( ( Taullr(Lmp,jsp),  Taulli(Lmp,jsp), Lmp = 1,nlma ), ( p, p, Lmp = nlma+1,nlm ), jsp = 1,nspinp )
      Taull(Lm,isp,1:nlma,:,iaabsi) = cmplx( Taullr(1:nlma,:),  Taulli(1:nlma,:), db )
    end do
    do Lm = nlma+1,nlm
      read(1,*)
    end do
  end do

  if( icheck > 1 ) then
    write(3,*)
    if( nspinp == 1 ) then
      write(3,270) Lmaxa, nspinp
      write(3,280) (( L, m, m = -L,L ), L = 0,Lmaxa )
    else
      write(3,290) Lmaxa, nspinp
      write(3,300) ((( L, m, i, m = -L,L ), L = 0,Lmaxa ), i = 1,nspinp)
    endif
    do isp = 1,nspinp
      do Lm = 1,nlma
        write(3,310) ( Taull(Lm,isp,1:nlma,jsp,iaabsi), jsp = 1,nspinp )
      end do
    end do
  endif
  
  if( ie == nenerg .and. .not. RIXS ) Close(1)
  
  return
  270 format(2i4,' = lmax, nspinp / ( l, m)')
  280 format(2420(13x,2i3,12x))
  290 format(2i4,' = lmax, nspinp / ( l, m, s)')
  300 format(2420(12x,3i3,10x))
  310 format(1p,2420(1x,2e15.7))
end

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

! Writing for extract

subroutine Ampl_writing(Ampl_dft,ie_computer,lmaxso,mpinodes,nlm_rixs,nlmomax,nspinp)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: ie_computer, isp, jsp, L, Lm, lmax, lmaxso, m, mpinodes, nlm_rixs, nlmo, &
            nlmomax, nspinp
 
  complex(kind=db), dimension(nlm_rixs,nspinp,nlmomax,nspinp,0:mpinodes-1):: Ampl_dft
  
  write(3,110)

! To save space in the bav file one limits the printing down to 3.
  lmax = nint( sqrt( real( nlm_rixs ) ) ) - 1 
  nlmo = (lmaxso + 1)**2

  write(3,260)
  if( nspinp == 1 ) then
    write(3,270) lmax, lmaxso, nspinp, nlmomax
    write(3,280) (( L, m, m = -L,L ), L = 0,lmaxso )
  else
    write(3,290) lmax, lmaxso, nspinp, nlmomax
    write(3,300) ((( L, m, isp, m = -L,L ), L = 0,lmaxso ), isp = 1,nspinp )
  endif
  do isp = 1,nspinp
    do Lm = 1,nlm_rixs
      write(3,310) ( Ampl_dft(Lm,isp,1:nlmo,jsp,ie_computer), jsp = 1,nspinp )
    end do
  end do

  return
  110 format(/' ---- Ampl --------',100('-'))
  260 format(/' Atomic amplitude')
  270 format(3i4,i6,' = lmax, lmxso, nspinp, nlmomax / ( l, m)')
  280 format(24200(17x,2i3,16x))
  290 format(3i4,i6,' = lmax, lmaxso, nspinp, nlmomax / ( l, m, s)')
  300 format(24200(16x,3i3,14x))
  310 format(1p,24200(1x,2e19.11))

end

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

function extract_nlmomax(multi_run,nom_fich_extract)

  use declarations
  implicit none

  integer:: i, istat, multi_run, extract_nlmomax, nlmomax

  character(len=17):: mot17
  character(len=132):: mot
  character(len=Length_name):: nom_fich_extract

  open(1, file = nom_fich_extract, status='old', iostat=istat)
  if( istat /= 0 ) call write_open_error(nom_fich_extract,istat,1)

  i = 0
  do
    read(1,'(A)' ) mot
    if( mot(2:15) /= 'Absorbing atom' ) cycle
    i = i + 1
    if( i == multi_run ) exit
  end do

  do
    read(1,'(A)' ) mot17
    if( mot17 == ' Atomic amplitude' ) exit
  end do

  read(1,*) i, i, i, nlmomax

  extract_nlmomax = nlmomax
     
  Close(1)

  return
end

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

subroutine Ampl_reading(Ampl_dft,icheck,ie,lmaxso,mpinodee,mpirank0,nenerg,nlm_probe,nlm_rixs,nlmomax,nspinp)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: icheck, ie, ipr, isp, jsp, L, Lm, lmax, lmax_m, lmax_probe, lmaxso, lmaxso_m, Lmp, m, mpinodee, mpirank0, nenerg, &
            nlm, nlm_probe, nlm_rixs, nlma, nlmo, nlmso_reg, nlmomax, nspinp

  character(len=17):: mot17
    
  complex(kind=db), dimension(nlm_rixs,nspinp,nlmomax,nspinp,0:mpinodee-1):: Ampl_dft

  real(kind=db):: p
  real(kind=db), dimension(nlmomax,nspinp):: Ampli, Amplr

  do
    read(1,'(A)' ) mot17
    if( mot17 == ' Atomic amplitude' ) exit
  end do

  read(1,*) lmax, lmaxso, nspinp
  
! Test on lmax value
  lmax_probe = nint( sqrt( real( nlm_probe ) ) ) - 1
  if( lmax_probe > lmax ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,110) lmax_probe, lmax
    end do
    stop
  endif

! Lmaxa can be < lmax because lmax is defined by the max of the lmax of the atoms with same atomic numer than the absorbing one 
  nlma = ( lmax + 1 )**2
  nlm = min( nlma, nlm_rixs )
  lmax_m = nint( sqrt( real( nlm ) ) ) - 1
  
  nlmso_reg = ( lmaxso + 1 )**2
  nlmo = min( nlmso_reg, nlmomax )
  lmaxso_m = nint( sqrt( real( nlmo ) ) ) - 1

  Ampl_dft(:,:,:,:,mpirank0) = ( 0._db, 0._db)

  read(1,*)
  
  do isp = 1,nspinp
    do Lm = 1,nlm
      read(1,*) ( ( Amplr(Lmp,jsp),  Ampli(Lmp,jsp), Lmp = 1,nlmo ), ( p, p, Lmp = nlmo+1,nlmso_reg ), jsp = 1,nspinp )
      Ampl_dft(Lm,isp,1:nlmo,:,mpirank0) = cmplx( Amplr(1:nlmo,:),  Ampli(1:nlmo,:), db )
    end do
    do Lm = nlm+1,nlma
      read(1,*)
    end do
  end do

  if( icheck > 1 ) then
    write(3,*)
    if( nspinp == 1 ) then
      write(3,270) lmax_m, lmaxso_m, nspinp
      write(3,280) (( L, m, m = -L,L ), L = 0,lmaxso_m )
    else
      write(3,290) lmax_m, lmaxso_m, nspinp
      write(3,300) ((( L, m, isp, m = -L,L ), L = 0,lmaxso_m ), isp = 1,nspinp)
    endif
    do isp = 1,nspinp
      do Lm = 1,nlm
        write(3,310) ( Ampl_dft(Lm,isp,1:nlmo,jsp,mpirank0), jsp = 1,nspinp )
      end do
    end do
  endif
  
  if( ie == nenerg ) Close(1)
  
  return
  110 format(//' the Lmax used in the extract file for the atomic amplitude in RIXS = ',i2, &
               ' is smaller than the one needed by the demanded multipolar component =',i2,//)
  270 format(3i4,' = lmax, lmaxso, nspinp / ( l, m)')
  280 format(24200(13x,2i3,12x))
  290 format(3i4,' = lmax, lmaxso, nspinp / ( l, m, s)')
  300 format(24200(12x,3i3,10x))
  310 format(1p,24200(1x,2e15.7))
end

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

! Writing of Tau

subroutine Tau_wout(iabsorig,Energ,ie,ie_computer,lmax,mpinodes,n_multi_run,nlm,nomfich,nspinp,Taull)

  use declarations
  implicit none
  include 'mpif.h'

  integer:: iabsorig, ie, ie_computer, L, lmax, mpinodes, n_multi_run, nlm, nspinp

  character(len=Length_name):: nomfich, nomficht
 
  complex(kind=db), dimension(nlm,nspinp,nlm,nspinp,0:mpinodes-1):: Taull

  real(kind=db):: Energ

  nomficht = nomfich
  if( n_multi_run > 1 ) then
    L = len_trim(nomficht)
    nomficht(L+1:L+1) = '_'
    call ad_number(iabsorig,nomficht,Length_name)
  endif
  L = len_trim(nomficht)
  nomficht(L+1:L+8) = '_tau.txt'

  if( ie == 1 ) then
    open(10,file = nomficht)
    write(10,110) (L, L, L = 0,lmax)
  else
    open(10,file = nomficht, position='append')
  endif   

  write(10,120) Energ*Rydb, ( Taull(1+L+L**2,1,1+L+L**2,1,ie_computer), L = 0,lmax )

  Close(10)
  
  return
  110 format('     Energy ',100('    Tau_r(',i1,')  ','   Tau_i(',i1,')  '))
  120 format(f12.5,1p,100(1x,2e13.5))

end

