! FDMNES subroutines

! Calculation of the cartesian tensors

! < g_1 | o*( L_s m_s, k_s, j_s, L_s, index_rank_s ) | f_1 > < f_2 | o*( L_i, m_i, k_i, j_i, L_i, index_rank_i  ) | g_2 >  

! n_Ec = ninit if TDDFT and XANES and Core_resolved
!      = nbseuil if TDDFT and XANES without core_resolved
!      = 1 in DFT XANES
! n_V = ninit if TDDFT and Core_resolved
!     = nbseuil if TDDFT without core_recolved
!     = 1 in DFT

subroutine Tenseur_car( Chi, Coef_g, Core_resolved, Ecinetic, Eimag, Energ, &
                  Enervide, Eseuil, Final_tddft, Full_potential, Green, Green_int, Hubb_a, Hubb_d, &
                  icheck, ie, ip_max, ip0, is_g, Lmax, Lmax_pot, Ldip, Lmoins1, Loct, Lplus1, Lqua, Lseuil, m_g, &
                  m_hubb, mpinodes, mpirank, mpirank0, msymdd, msymddi, msymdq, msymdqi, msymdo, msymdoi, msymoo, msymooi, &
                  msymqq, msymqqi, Multipole, n_Ec, n_oo, n_V, nb_test, nbseuil, ns_dipmag, nenerg_tddft, ninit1, &
                  ninit, ninit_tddft, ninitr, ninitv, nlm_pot, nlm_probe, nlm_p_fp, nlmamax, nr, nrm, nspin, nspino, &
                  nspinp, Numat, Psii, r, Relativiste, Renorm, Rmtg, Rmtsd, rof0, Rot_atom_abs, Rot_cluster, &
                  T_E1E1, T_E1E1_m, T_E1E2, T_E1E2_m, T_E1E3, T_E1E3_m, T_E1M1, T_E1M1_m, T_E1M2, T_E1M2_m, T_E1SP, T_E1SP_m, &
                  T_E2E2, T_E2E2_m, T_E2M2, T_E2M2_m, T_E3E3, T_E3E3_m, T_M1M1, T_M1M1_m, T_SPSP, T_SPSP_m, &
                  Solsing, Solsing_only, Spinorbite, Taull, Tddft, Test_case, V_hubb, V_intmax, V0bd, Vrato, Ylm_comp )

  use declarations
  implicit none

  integer, parameter:: Lomax = 3
  integer, parameter:: Lmomax = ( Lomax + 1 )**2

  real(kind=db), parameter:: Four_mc2 = 4 * 2 / alfa_sf**2  ! in Rydberg  

  integer:: h_s, hh_s, hh_i, h_i, i, icheck, ie, index_Multipole, initr, initlt, ip, ip_max, ip0, index_rank_i, index_rank_s, &
    irang1, is, iseuil, isi, isol, isp, j, j_i, j_s, jh_i, jh_s, jj_i, jj_s, jjh_i, jjh_s, k_i, k_i_SpinPos, k_s, &
    k_s_SpinPos, kk_i, kk_s, L_i, L_s, lm, Lm_i, Lm_s, Lmax, Lmax_pot, Lseuil, m_hubb, m_i, m_k_i, m_k_s, m_s, mpinodes, &
    mpirank, mpirank0, n_Ec, n_oo, n_V, nb_test, nbseuil, nenerg_tddft, ninit1, ninit, ninit_tddft, &
    ninitr, ninitv, nh_i, nh_s, nj_i, nj_s, nlm_p_fp, nlm_pot, nlm_probe, &
    nlmamax, nr, nrang, nrm, ns_dipmag, nspin, nspino, nspinp, Numat

  integer, dimension(ninit,2):: m_g
  integer, dimension(ninit):: is_g
  integer, dimension(3):: Ldip
  integer, dimension(3,3):: Lqua, msymdd, msymddi
  integer, dimension(3,3,3):: Loct, msymdq, msymdqi
  integer, dimension(3,3,3,3):: msymdo, msymdoi, msymqq, msymqqi
  integer, dimension(3,n_oo,3,n_oo):: msymoo, msymooi

  Character(len=4), dimension(n_Multipole_max+1):: Multipole_name
  
  complex(kind=db), dimension(ninitr):: Te, Te_m, Ten, Ten_m, Tens, Tens_m
  complex(kind=db), dimension(3,3) :: Mat2
  complex(kind=db), dimension(3,3,3) :: Mat3
  complex(kind=db), dimension(3,3,3,3) :: Mat4
  complex(kind=db), dimension(3,3,3,3,3,3) :: Mat6
  complex(kind=db), dimension(3,3,ninitr,0:mpinodes-1):: T_E1E1, T_E1E1_m, T_E1SP, T_E1SP_m, T_E1M1, T_E1M1_m, T_M1M1, T_M1M1_m, &
                                                         T_SPSP, T_SPSP_m
  complex(kind=db), dimension(3,3,3,ninitr,0:mpinodes-1):: T_E1E2, T_E1E2_m, T_E1M2, T_E1M2_m
  complex(kind=db), dimension(3,3,3,3,ninitr,0:mpinodes-1):: T_E1E3, T_E1E3_m, T_E2E2, T_E2E2_m, T_E2M2, T_E2M2_m
  complex(kind=db), dimension(3,n_oo,3,n_oo,ninitr,0:mpinodes-1):: T_E3E3, T_E3E3_m
  complex(kind=db), dimension(lmomax,lmomax,ninitr):: Tens_lm, Tens_lm_m
  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp):: Taull
  complex(kind=db), dimension(nenerg_tddft,nlmamax,nspinp,nspino,nbseuil):: rof0
  complex(kind=db), dimension(-m_hubb:m_hubb,nspinp,-m_hubb:m_hubb,nspinp):: V_hubb
  complex(kind=db), dimension(nlm_probe*nspino,nlm_probe*nspino,ninit_tddft,ninit_tddft,2,2,ns_dipmag):: Chi

  complex(kind=db), dimension(:,:,:,:,:,:,:), allocatable:: Singul
  complex(kind=db), dimension(:,:,:,:,:,:), allocatable:: rof

  logical:: Core_resolved, Dip_rel, E1E1, E1E2, E1E3, E1M1, E1M2, E1SP, E2E2, E2M2, E3E3, Final_tddft, Full_potential, &
    Green, Green_int, Hubb_a, Hubb_d, Lmoins1, Lplus1, M1M1, No_recopy, NRIXS, Relativiste, &
    Renorm, Solsing, Solsing_only, Spinorbite, SpinPos_i, SpinPos_s, SPSP, Tddft, Ylm_comp

  logical, dimension(n_multipole_max):: Multipole
  logical, dimension(nb_test):: Test_case

  real(kind=db):: clm_s, clm_i, Eimag, Energ, Enervide_t, Ephoton, Rmtg, Rmtsd, V_intmax

  real(kind=db), dimension(nspin):: Ecinetic_e, V0bd_e
  real(kind=db), dimension(nbseuil):: Eseuil, Omega
  real(kind=db), dimension(ninitr):: Tensi, Tensr
  real(kind=db), dimension(-1:lomax,3,3,3,lmomax):: Clm
  real(kind=db), dimension(3,3):: Rot_atom_abs, Rot_cluster, Rot_tem
  real(kind=db), dimension(ninit,2):: Coef_g
  real(kind=db), dimension(nspin,n_Ec):: Ecinetic
  real(kind=db), dimension(n_Ec):: Enervide
  real(kind=db), dimension(nspin,n_V):: V0bd
  real(kind=db), dimension(nr):: r
  real(kind=db), dimension(nr,nlm_pot,nspin):: Vrato_e
  real(kind=db), dimension(nrm,nbseuil):: psii
  real(kind=db), dimension(nr,nlm_pot,nspin,n_V):: Vrato
  real(kind=db), dimension(nr,ip0:ip_max):: r_or_bess

  if( icheck > 1 ) write(3,100)

  data Multipole_name/ 'E1E1', 'E1E2', 'E1E3', 'E1M1', 'E1M2', 'E2E2', 'E3E3', 'M1M1', 'M1M2', 'M2M2', 'E1SP', 'SPSP', 'E2M2', &
                       'SPE1' /

  E1E1 = Multipole(1); E1E2 = Multipole(2); E1E3 = Multipole(3)
  E1M1 = Multipole(4); E1M2 = Multipole(5)
  E2E2 = Multipole(6);
  E3E3 = Multipole(7); M1M1 = Multipole(8)
  E1SP = Multipole(11); SPSP = Multipole(12)
  E2M2 = Multipole(13)
  
  NRIXS = .false.
  Dip_rel = E1SP .or. SPSP

! Calculation of radial integrals for regular and irregular solutions

  do ip = ip0,ip_max
    select case(ip)
      case(0)
        r_or_bess(1:nr,ip) = 1._db
      case(1)
        r_or_bess(1:nr,ip) = r(1:nr)
      case default
        r_or_bess(1:nr,ip) = r(1:nr)**ip
    end select
  end do    
  
  allocate( Singul(nlm_probe,nspinp,nlm_probe,nspinp,ip0:ip_max,ip0:ip_max,ninitv) )
  allocate( rof(nlm_probe,nlm_p_fp,nspinp,nspino,ip0:ip_max,ninitv) )

  rof(:,:,:,:,:,:) = (0._db, 0._db)
  Singul(:,:,:,:,:,:,:) = (0._db, 0._db)

  do initlt = 1,n_Ec
    if( Final_tddft ) then
      Ecinetic_e(:) = Ecinetic(:,initlt)
    else
      Ecinetic_e(:) = Ecinetic(:,1)
    endif
    if( Final_tddft ) then
      do isp = 1,nspin
        V0bd_e(isp) = V0bd(isp,initlt)
        Vrato_e(1:nr,:,isp) = Vrato(1:nr,:,isp,initlt)
      end do
    else
      do isp = 1,nspin
        V0bd_e(isp) = V0bd(isp,1)
        Vrato_e(1:nr,:,isp) = Vrato(1:nr,:,isp,1)
      end do
    endif
    Enervide_t = Enervide(initlt)
    
    call Radial_integral(Ecinetic_e,Eimag,Energ,Enervide_t,Eseuil,Final_tddft,Full_potential,Hubb_a,Hubb_d,icheck, &
       initlt,ip_max,ip0,Lmax,Lmax_pot,m_hubb,nb_test,nbseuil,ninit1,ninitv,nlm_pot,nlm_probe,nlm_p_fp,nr,NRIXS,nrm, &
       nspin,nspino,nspinp,Numat,psii,r,r_or_bess,Relativiste,Renorm,Rmtg,Rmtsd,rof,Singul,Solsing,Spinorbite,Test_case, &
       V_hubb,V_intmax,V0bd_e,Vrato_e,Ylm_comp)

  end do

  if( Tddft .and. .not. Final_tddft ) then
    do isol = 1,nspino
      do isp = 1,nspinp
        do lm = 1,nlm_probe
          do i = 1,nbseuil
            if( nlm_p_fp == 1 ) then
              rof0(ie,lm,isp,isol,i) = rof(lm,1,isp,isol,0,i)
            else
! ce qui sert de base en tddft, ce sont les (L_i,m_i,s) et pas les (l,m) vrais.
              rof0(ie,lm,isp,isol,:) = sum( rof(1:nlm_probe,lm,isp,isol,0,i) )
            endif
          end do
        end do
      end do
    end do
  endif

  if( icheck > 1 ) write(3,110)

  if( E1E3 .or. E3E3 ) then
    nrang = 3
  elseif( E1E2 .or. E2E2 .or. E2M2 ) then
    nrang = 2
  else
    nrang = 1
  endif
  if( E1M2 .or. E2M2 ) then
    irang1 = -1
  elseif( E1M1 .or. M1M1 ) then
    irang1 = 0
  else
    irang1 = 1
  endif

! Calculation of the componants of the transformation cartesian basis - spherical basis
  call Tr_Cartesian_Spherical_basis( Clm, irang1, Lmomax, Lomax, mpirank0, nrang )

  T_E1E1(:,:,:,mpirank) = (0._db,0._db)
  T_E1E2(:,:,:,:,mpirank) = (0._db,0._db)
  T_E1E3(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_E1M1(:,:,:,mpirank) = (0._db,0._db)
  T_E1M2(:,:,:,:,mpirank) = (0._db,0._db)
  T_E1SP(:,:,:,mpirank) = (0._db,0._db)
  T_E2E2(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_E2M2(:,:,:,:,:,mpirank) = (0._db,0._db)
  if( E3E3 ) T_E3E3(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_M1M1(:,:,:,mpirank) = (0._db,0._db)
  T_SPSP(:,:,:,mpirank) = (0._db,0._db)

! At this stage, only useful if magnetic calculation and Greenint (that is tddft calculation) 
  T_E1E1_m(:,:,:,mpirank) = (0._db,0._db)
  T_E1E2_m(:,:,:,:,mpirank) = (0._db,0._db)
  T_E1E3_m(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_E1M1_m(:,:,:,mpirank) = (0._db,0._db)
  T_E1M2_m(:,:,:,:,mpirank) = (0._db,0._db)
  T_E1SP_m(:,:,:,mpirank) = (0._db,0._db)
  T_E2E2_m(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_E2M2_m(:,:,:,:,:,mpirank) = (0._db,0._db)
  if( E3E3 ) T_E3E3_m(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_M1M1_m(:,:,:,mpirank) = (0._db,0._db)
  T_SPSP_m(:,:,:,mpirank) = (0._db,0._db)
  
  if( E1SP .or. SPSP ) then
    do iseuil = 1,nbseuil
      Ephoton = Energ + Eseuil(iseuil)
! For very low energy edges (to avoid negative photon energy) :
      Ephoton = max( 0.001_db/Rydb, Ephoton)
      Omega(iseuil) = Ephoton / Four_mc2
    end do 
  endif

  if( icheck > 1 ) write(3,120) Lseuil

! Loop over multipoles (the +1 is to calculate SPE1 which is add to E1SP)
  do index_Multipole = 1,n_Multipole_max+1

    if( index_Multipole == 14 ) then
      if( .not. Multipole( 11 ) ) cycle
    else
      if( .not. Multipole( index_Multipole ) ) cycle
    endif

    SpinPos_i = .false.
    SpinPos_s = .false.

    select case( index_Multipole )
      case(1) ! E1E1
        index_rank_i = 1
        index_rank_s = 1
      case(2) ! E1E2
        index_rank_i = 2
        index_rank_s = 1
      case(3) ! E1E3
        index_rank_i = 3
        index_rank_s = 1
      case(4) ! E1M1
        index_rank_i = 0
        index_rank_s = 1
      case(5) ! E1M2
        index_rank_i = -1
        index_rank_s = 1
      case(6) ! E2E2
        index_rank_i = 2
        index_rank_s = 2
      case(7) ! E3E3
        index_rank_i = 3
        index_rank_s = 3
      case(8) ! M1M1
        index_rank_i = 0
        index_rank_s = 0
      case(9) ! M1M2  ! not coded
        index_rank_i = -1
        index_rank_s = 0
      case(10) ! M2M2  not coded
        index_rank_i = -1
        index_rank_s = -1
      case(11) ! E1SP
        index_rank_i = 1
        index_rank_s = 1
        SpinPos_i = .true.
      case(12) ! SPSP
        index_rank_i = 1
        index_rank_s = 1
        SpinPos_i = .true.
        SpinPos_s = .true.
      case(13) ! E2M2
        index_rank_i = -1
        index_rank_s = 2
      case(14) ! SPE1 which will be add to E1SP
        index_rank_i = 1
        index_rank_s = 1
        SpinPos_s = .true.
    end select

    Tens_lm(:,:,:) = (0._db,0._db)
    Tens_lm_m(:,:,:) = (0._db,0._db)

    if( icheck > 1 ) write(3,130) Multipole_name(index_Multipole), index_rank_s, index_rank_i

    No_recopy = .not. ( index_Multipole <= 3 .or. index_Multipole == 6 .or. index_Multipole == 7 )
    
! note that in tensors outgoing (s) are first index and incoming (i) second index  ( line and column when rank 2 )
! Thus "E1E2" is "E1_s-E2_i", same all the others, but E1SP which is E1SP + SPE1 

! Loops over tensor indices
    do k_s = 1,3

      if( index_Multipole == 1 .and. Ldip(k_s) /= 1 ) cycle

! Dipole position index
      do k_s_SpinPos = 1,3
        if( ( .not. SpinPos_s .and. k_s_SpinPos /= k_s ) .or. ( SpinPos_s .and. k_s_SpinPos == k_s ) ) cycle

        if( index_rank_s > 1 .or. index_rank_s < 0 ) then
          nj_s = 3
        else
          nj_s = 1
        endif

        do j_s = 1,nj_s

          if( index_rank_s == -1 ) then
 ! value of "m" for the pseudo Gaunt coefficient in M2, corresponding the the wave vector orientation
            if( j_s == 1 ) then
              m_k_s = 1
            elseif( j_s == 2 ) then
              m_k_s = -1
            else
              m_k_s = 0
            endif
          else
            m_k_s = 0 !not used in that case
          endif
             
          if( index_Multipole == 6 .and. Lqua(k_s,j_s) /= 1 ) cycle

          if( index_rank_s == 3 ) then
            nh_s = 3
          else
            nh_s = 1
          endif

          do h_s = 1,nh_s

            if( index_rank_s == 3 .and. Loct(k_s,j_s,h_s) /= 1 ) cycle
            jh_s = 3 * ( j_s - 1 ) + h_s

       ! incoming
            do k_i = 1,3

              if( index_Multipole <= 3 .and. Ldip(k_i) /= 1 ) cycle

! Dipole position index
              do k_i_SpinPos = 1,3
                if( ( .not. SpinPos_i .and. k_i_SpinPos /= k_i ) .or. ( SpinPos_i .and. k_i_SpinPos == k_i ) ) cycle

                if( index_rank_i > 1 .or. index_rank_i < 0 ) then
                  nj_i = 3
                else
                  nj_i = 1
                endif

                do j_i = 1,nj_i
            
                  if( index_rank_i == -1 ) then
 ! value of "m" for the pseudo Gaunt coefficient in M2, corresponding the the wave vector orientation
                    if( j_i == 1 ) then
                      m_k_i = 1
                    elseif( j_i == 2 ) then
                      m_k_i = -1
                    else
                      m_k_i = 0
                    endif
                  else
                    m_k_i = 0 ! not used in that case
                  endif

                  if( index_multipole == 6 .and. Lqua(k_i,j_i) /= 1  ) cycle
                
                  if( index_rank_i == 3 ) then
                    nh_i = 3
                  else
                    nh_i = 1
                  endif
               
                  do h_i = 1,nh_i
                
                    if( index_rank_i == 3 .and. Loct(k_i,j_i,h_i) /= 1 ) cycle
                    jh_i = 3 * ( j_i - 1 ) + h_i
                
                    if( ( index_multipole == 1 .or. index_Multipole == 8 ) .and. k_i < k_s ) cycle
 
                    if( index_multipole == 1 ) then
                      if( msymdd(k_s,k_i) == 0 .and. msymddi(k_s,k_i) == 0 ) cycle
                      if( sum( abs( T_E1E1(k_s,k_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                    elseif( index_rank_s == 1 .and. index_rank_i == 2 ) then
                      if( msymdq(k_s,k_i,j_i) == 0 .and. msymdqi(k_s,k_i,j_i) == 0 ) cycle
                      if( sum( abs( T_E1E2(k_s,k_i,j_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                    elseif( index_rank_s == 2 .and. index_rank_i == 2 ) then
                      if( msymqq(k_s,j_s,k_i,j_i) == 0 .and. msymqqi(k_s,j_s,k_i,j_i) == 0 ) cycle
                      if( sum( abs( T_E2E2(k_s,j_s,k_i,j_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                    elseif( index_rank_s == 1 .and. index_rank_i == 3 ) then
                      if( msymdo(k_s,k_i,j_i,h_i) == 0 .and. msymdoi(k_s,k_i,j_i,h_i) == 0 ) cycle
                    elseif( index_rank_s == 3 .and. index_rank_i == 3 ) then
                      if( msymoo(k_s,jh_s,k_i,jh_i) == 0 .and. msymooi(k_s,jh_s,k_i,jh_i) == 0 ) cycle
                      if( sum( abs( T_E3E3(k_s,jh_s,k_i,jh_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                    endif
                
                    Tens(:) = (0._db,0._db)
                    Tens_m(:) = (0._db,0._db)
               
                    if( icheck > 1 ) write(3,140) k_s, j_s, h_s, k_i, j_i, h_i
               
! Loops over spherical components of the tensors
                    Lm_s = 0
                    do L_s = 0,max(index_rank_s,1)
                      do m_s = -L_s,L_s
                        Lm_s = Lm_s + 1    ! for M2, correponds only to wave vector k_s

                        if( index_rank_s == - 1 ) then
                          if( L_s == 0 ) cycle
                          if( ( k_s == 1 .and. m_s /= 1 ) .or. ( k_s == 2 .and. m_s /= -1 ) .or. ( k_s == 3 .and. m_s /= 0 ) ) cycle 
                          clm_s = clm( index_rank_s, k_s, j_s, h_s, L_s**2 + L_s + 1 + m_k_s )
                        elseif( SpinPos_s ) then
                          clm_s = clm( index_rank_s, k_s_SpinPos, j_s, h_s, Lm_s )
                        else
                          clm_s = clm( index_rank_s, k_s, j_s, h_s, Lm_s )
                        endif
                        if( abs(clm_s) < eps10 ) cycle
               
                        Lm_i = 0
                        do L_i = 0,max(index_rank_i,1)
                          do m_i = -L_i,L_i
                            Lm_i = Lm_i + 1
               
                            if( index_rank_i == - 1 ) then
                              if( L_i == 0 ) cycle
                          ! here  m_i is the one of index k_i that is of the direction of (2/3L + 2S) in M2
                          ! m_k_i is the one of j_i, that is the one of the wave vector in M2
                              if( ( k_i == 1 .and. m_i /= 1 ) .or. ( k_i == 2 .and. m_i /= -1 ) &
                                                               .or. ( k_i == 3 .and. m_i /= 0 ) ) cycle 
                              clm_i = clm( index_rank_i, k_i, j_i, h_i, L_i**2 + L_i + 1 + m_k_i )
                            elseif( SpinPos_i ) then
                              clm_i = clm( index_rank_i, k_i_SpinPos, j_i, h_i, Lm_i )
                            else
                              clm_i = clm( index_rank_i, k_i, j_i, h_i, Lm_i )
                            endif
                            if( abs(clm_i) < eps10 ) cycle

! Calculation of the tensor component
                            if( sum( abs( Tens_lm(Lm_s,Lm_i,:) ) ) < 1.e-15_db .or. No_recopy ) then

                              if( icheck > 1 ) write(3,150) L_s, m_s, L_i, m_i, clm_s, clm_i
                              call Tens_ab( Chi, Coef_g, Core_resolved, Final_tddft, Full_potential, &
                                Green, Green_int, icheck, index_rank_i, index_rank_s, ip_max, ip0, is_g, k_i, k_s, L_s, &
                                m_k_i, m_k_s, m_s, L_i, m_g, m_i, Lmax, Lmoins1, Lplus1, Lseuil, nbseuil, ns_dipmag, ninit1, &
                                ninit, ninit_tddft, ninitv, ninitr, nlm_probe, nlm_p_fp, NRIXS, nspinp, nspino, Omega, rof, &
                                Singul, Solsing, Solsing_only, Spinorbite, SpinPos_i, SpinPos_s, Taull, Ten, Ten_m, Ylm_comp )
                          
                              if( Numat == 1) then
                         ! For hydrogen there is only 1 core state
                                Ten(:) = Ten(:) / 2
                                Ten_m(:) = Ten_m(:) / 2
                              endif

                              Tens_lm(Lm_s,Lm_i,:) = Ten(:)
                              Tens_lm_m(Lm_s,Lm_i,:) = Ten_m(:)

                            endif

                            Tens(:) = Tens(:) + clm_s * clm_i * Tens_lm(Lm_s,Lm_i,:)
                            Tens_m(:) = Tens_m(:) + clm_s * clm_i * Tens_lm_m(Lm_s,Lm_i,:)
                          end do
                        end do ! end of loop over L_i
                      end do
                    end do ! end of loop over L_s

! Filling of the calculted value in all the tensor components equivalent by symmetry 

                    Tensr(:) = real( Tens(:), db )
                    Tensi(:) = aimag( Tens(:) )

! M1-M1 (dipole magnetique - dipole magnetique)
                    if( index_rank_s == 0 .and. index_rank_i == 0 ) then

                      T_M1M1(k_s,k_i,:,mpirank) = Tens(:)
                      if( Green_int ) then
                        T_M1M1(k_i,k_s,:,mpirank) = Tens(:)
                        T_M1M1_m(k_s,k_i,:,mpirank) = Tens_m(:)
                        T_M1M1_m(k_i,k_s,:,mpirank) = - Tens_m(:)
                      else
                        T_M1M1(k_i,k_s,:,mpirank) = conjg( Tens(:) )
                      endif
                    endif

! E1-M1 (dipole electric - dipole magnetic )
                    if( index_rank_s == 1 .and. index_rank_i == 0 ) then
                      T_E1M1(k_s,k_i,:,mpirank) = Tens(:)
                      if( Green_int ) T_E1M1_m(k_s,k_i,:,mpirank) = Tens_m(:)
                    endif

! E1-M2 ( dipole electric - quadrupole magnetic)
                    if( index_rank_s == 1 .and. index_rank_i == -1 ) then
                      T_E1M2(k_s,k_i,j_i,:,mpirank) = Tens(:)
                      if( Green_int ) T_E1M2_m(k_s,k_i,j_i,:,mpirank) = Tens_m(:)
                    endif

! E2-M2 ( quadrupole electric - quadrupole magnetic)
                    if( index_rank_s == 2 .and. index_rank_i == -1 ) then
                      T_E2M2(k_s,j_s,k_i,j_i,:,mpirank) = Tens(:)
                      if( Green_int ) T_E2M2_m(k_s,j_s,k_i,j_i,:,mpirank) = Tens_m(:)
                    endif

! E1E1
                    if( index_Multipole == 1 ) then

                      do kk_s = 1,3
                        do kk_i = kk_s,3
                          if( abs(msymdd(kk_s,kk_i)) /= abs(msymdd(k_s,k_i)) .or. abs(msymddi(kk_s,kk_i)) &
                             /= abs(msymddi(k_s,k_i)) ) cycle
                          if( msymdd(kk_s,kk_i) /= 0 ) then
                            is = msymdd(kk_s,kk_i) / msymdd(k_s,k_i)
                          else
                            is = 0
                          endif
                          if( msymddi(kk_s,kk_i) /= 0 ) then
                            isi = msymddi(kk_s,kk_i) / msymddi(k_s,k_i)
                          else
                            isi = 0
                          endif
                          if( Green_int ) then
                            T_E1E1(kk_s,kk_i,:,mpirank) = is * Tens(:)
                            T_E1E1(kk_i,kk_s,:,mpirank) = is * Tens(:)
                            T_E1E1_m(kk_s,kk_i,:,mpirank) = isi * Tens_m(:)
                            T_E1E1_m(kk_i,kk_s,:,mpirank) = -isi * Tens_m(:)
                          else
! tenseur hermitique
                            Te(:) = cmplx( is*Tensr(:), isi*Tensi(:), db)
                            T_E1E1(kk_s,kk_i,:,mpirank) = Te(:)
                            T_E1E1(kk_i,kk_s,:,mpirank) = conjg( Te(:) )
                          endif
                        end do
                      end do

! E1-SP
                    elseif( index_rank_s == 1 .and. index_rank_i == 1 .and. SpinPos_i .and. .not. SpinPos_s ) then

                      T_E1SP(k_s,k_i,:,mpirank) = T_E1SP(k_s,k_i,:,mpirank) + Tens(:)
                      if( Green_int ) T_E1SP_m(k_s,k_i,:,mpirank) = T_E1SP_m(k_s,k_i,:,mpirank) + Tens_m(:)

                    elseif( index_rank_s == 1 .and. index_rank_i == 1 .and. .not. SpinPos_i .and. SpinPos_s ) then

                      T_E1SP(k_s,k_i,:,mpirank) = T_E1SP(k_s,k_i,:,mpirank) + Tens(:)
                      if( Green_int ) T_E1SP_m(k_s,k_i,:,mpirank) = T_E1SP_m(k_s,k_i,:,mpirank) + Tens_m(:)

! SP-SP
                    elseif( index_Multipole == 12 ) then

                      T_SPSP(k_s,k_i,:,mpirank) = T_SPSP(k_s,k_i,:,mpirank) + Tens(:)
                      if( Green_int ) T_SPSP_m(k_s,k_i,:,mpirank) = T_SPSP_m(k_s,k_i,:,mpirank) + Tens_m(:)

! Dipole-quadrupole
                    elseif( index_Multipole == 2 ) then
                      do kk_s = 1,3
                        do kk_i = 1,3
                          do jj_i = kk_i,3
                            if( abs(msymdq(kk_s,kk_i,jj_i)) /= abs(msymdq(k_s,k_i,j_i)) .or. abs(msymdqi(kk_s,kk_i,jj_i)) &
                               /= abs(msymdqi(k_s,k_i,j_i)) ) cycle
                            if( msymdq(kk_s,kk_i,jj_i) /= 0 ) then
                              is = msymdq(kk_s,kk_i,jj_i) / msymdq(k_s,k_i,j_i)
                            else
                              is = 0
                            endif
                            if( msymdqi(kk_s,kk_i,jj_i) /= 0 ) then
                              isi = msymdqi(kk_s,kk_i,jj_i) / msymdqi(k_s,k_i,j_i)
                            else
                              isi = 0
                            endif
                            if( Green_int ) then
                              T_E1E2(kk_s,kk_i,jj_i,:,mpirank) = is * Tens(:)
                              T_E1E2(kk_s,jj_i,kk_i,:,mpirank) = is * Tens(:)
                              T_E1E2_m(kk_s,kk_i,jj_i,:,mpirank) = isi * Tens_m(:)
                              T_E1E2_m(kk_s,jj_i,kk_i,:,mpirank) = isi * Tens_m(:)
                            else
                              Te(:) = cmplx( is*Tensr(:), isi*Tensi(:), db )
                              T_E1E2(kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                              T_E1E2(kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                            endif
                          end do
                        end do
                      end do

! Dipole-octupole
                    elseif( index_multipole == 3 ) then
                      do kk_s = 1,3
                        do kk_i = 1,3
                          do jj_i = 1,3
                            do hh_i = 1,3
                              if( abs( msymdo(kk_s,kk_i,jj_i,hh_i) ) /= abs( msymdo(k_s,k_i,j_i,h_i) )  &
                                .or. abs( msymdoi(kk_s,kk_i,jj_i,hh_i) ) /= abs( msymdoi(k_s,k_i,j_i,h_i) ) ) cycle
                              if( msymdo(kk_s,kk_i,jj_i,hh_i) /= 0 ) then
                                is = msymdo(kk_s,kk_i,jj_i,hh_i) / msymdo(k_s,k_i,j_i,h_i)
                              else
                                is = 0
                              endif
                              if( msymdoi(kk_s,kk_i,jj_i,hh_i)/=0) then
                                isi = msymdoi(kk_s,kk_i,jj_i,hh_i) / msymdoi(k_s,k_i,j_i,h_i)
                              else
                                isi = 0
                              endif
                              if( Green_int ) then
                                T_E1E3(kk_s,kk_i,jj_i,hh_i,:,mpirank) = is * Tens(:)
                                T_E1E3(kk_s,kk_i,hh_i,jj_i,:,mpirank) = is * Tens(:)
                                T_E1E3(kk_s,jj_i,kk_i,hh_i,:,mpirank) = is * Tens(:)
                                T_E1E3(kk_s,hh_i,kk_i,jj_i,:,mpirank) = is * Tens(:)
                                T_E1E3(kk_s,jj_i,hh_i,kk_i,:,mpirank) = is * Tens(:)
                                T_E1E3(kk_s,hh_i,jj_i,kk_i,:,mpirank) = is * Tens(:)
                                T_E1E3_m(kk_s,kk_i,jj_i,hh_i,:,mpirank) = isi * Tens_m(:)
                                T_E1E3_m(kk_s,kk_i,hh_i,jj_i,:,mpirank) = isi * Tens_m(:)
                                T_E1E3_m(kk_s,jj_i,kk_i,hh_i,:,mpirank) = isi * Tens_m(:)
                                T_E1E3_m(kk_s,hh_i,kk_i,jj_i,:,mpirank) = isi * Tens_m(:)
                                T_E1E3_m(kk_s,jj_i,hh_i,kk_i,:,mpirank) = isi * Tens_m(:)
                                T_E1E3_m(kk_s,hh_i,jj_i,kk_i,:,mpirank) = isi * Tens_m(:)
                              else
                                Te(:) = cmplx(is*Tensr(:),isi*Tensi(:),db)
                                T_E1E3(kk_s,kk_i,jj_i,hh_i,:,mpirank) = Te(:)
                                T_E1E3(kk_s,kk_i,hh_i,jj_i,:,mpirank) = Te(:)
                                T_E1E3(kk_s,jj_i,kk_i,hh_i,:,mpirank) = Te(:)
                                T_E1E3(kk_s,hh_i,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E1E3(kk_s,jj_i,hh_i,kk_i,:,mpirank) = Te(:)
                                T_E1E3(kk_s,hh_i,jj_i,kk_i,:,mpirank) = Te(:)
                              endif
                            end do
                          end do
                        end do
                      end do

! Quadrupole-quadrupole
                    elseif( index_rank_s == 2 .and. index_rank_i == 2 ) then
                      do kk_s = 1,3
                        do jj_s = 1,3
                          do kk_i = 1,3
                            do jj_i = 1,3
                              if( sum( abs( T_E2E2(kk_s,jj_s,kk_i,jj_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                              if( abs( msymqq(kk_s,jj_s,kk_i,jj_i) ) /= abs( msymqq(k_s,j_s,k_i,j_i) ) &
                                  .or. abs( msymqqi(kk_s,jj_s,kk_i,jj_i) ) /= abs( msymqqi(k_s,j_s,k_i,j_i) ) ) cycle
                              if( msymqq(kk_s,jj_s,kk_i,jj_i) /= 0 ) then
                                is = msymqq(kk_s,jj_s,kk_i,jj_i) / msymqq(k_s,j_s,k_i,j_i)
                              else
                                is = 0
                              endif
                              if( msymqqi(kk_s,jj_s,kk_i,jj_i) /= 0 ) then
                                isi = msymqqi(kk_s,jj_s,kk_i,jj_i) / msymqqi(k_s,j_s,k_i,j_i)
                              else
                                isi = 0
                              endif
                              if( Green_int ) then
                                Te(:) = is * Tens(:)
                                T_E2E2(kk_s,jj_s,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E2E2(jj_s,kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E2E2(kk_s,jj_s,jj_i,kk_i,:,mpirank) = Te(:)
                                T_E2E2(jj_s,kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                                T_E2E2(kk_i,jj_i,kk_s,jj_s,:,mpirank) = Te(:)
                                T_E2E2(jj_i,kk_i,kk_s,jj_s,:,mpirank) = Te(:)
                                T_E2E2(kk_i,jj_i,jj_s,kk_s,:,mpirank) = Te(:)
                                T_E2E2(jj_i,kk_i,jj_s,kk_s,:,mpirank) = Te(:)
                                Te(:) = isi * Tens_m(:)
                                T_E2E2_m(kk_s,jj_s,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E2E2_m(jj_s,kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E2E2_m(kk_s,jj_s,jj_i,kk_i,:,mpirank) = Te(:)
                                T_E2E2_m(jj_s,kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                                T_E2E2_m(kk_i,jj_i,kk_s,jj_s,:,mpirank) = - Te(:)
                                T_E2E2_m(jj_i,kk_i,kk_s,jj_s,:,mpirank) = - Te(:)
                                T_E2E2_m(kk_i,jj_i,jj_s,kk_s,:,mpirank) = - Te(:)
                                T_E2E2_m(jj_i,kk_i,jj_s,kk_s,:,mpirank) = - Te(:)
                              else
                                Te(:)= cmplx( is*Tensr(:), isi*Tensi(:),db)
                                T_E2E2(kk_s,jj_s,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E2E2(jj_s,kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                                T_E2E2(kk_s,jj_s,jj_i,kk_i,:,mpirank) = Te(:)
                                T_E2E2(jj_s,kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                                Te(:) = Conjg( Te(:) )
                                T_E2E2(kk_i,jj_i,kk_s,jj_s,:,mpirank) = Te(:)
                                T_E2E2(jj_i,kk_i,kk_s,jj_s,:,mpirank) = Te(:)
                                T_E2E2(kk_i,jj_i,jj_s,kk_s,:,mpirank) = Te(:)
                                T_E2E2(jj_i,kk_i,jj_s,kk_s,:,mpirank) = Te(:)
                              endif
                            end do
                          end do
                        end do
                      end do
! Octupole-octupole
                    elseif( index_rank_s == 3 .and. index_rank_i == 3 ) then
                
                      do kk_s = 1,3
                        do jj_s = 1,3
                          do hh_s = 1,3
                            jjh_s = 3 * ( jj_s - 1 ) + hh_s
                            do kk_i = 1,3
                              do jj_i = 1,3
                                do hh_i = 1,3
                                  jjh_i = 3 * ( jj_i - 1 ) + hh_i
                                  i = msymoo(kk_s,jjh_s,kk_i,jjh_i)
                                  j = msymoo(k_s,jh_s,k_i,jh_i)
                                  if( abs(i) /= abs(j) ) cycle
                                  if( j /= 0 ) then
                                    is = i / j
                                  else
                                    is = 0
                                  endif
                                  i = msymooi(kk_s,jjh_s,kk_i,jjh_i)
                                  j = msymooi(k_s,jh_s,k_i,jh_i)
                                  if( abs(i) /= abs(j) ) cycle
                                  if( j /= 0 ) then
                                    isi = i / j
                                  else
                                    isi = 0
                                  endif
                                  if( Green_int ) then
                                    Te(:) = is * Tens(:)
                                    Te_m(:) = isi * Tens_m(:)
                                  else
                                    Te(:)= cmplx( is * Tensr(:), isi * Tensi(:), db )
                                    Te_m(:) = Conjg( Te(:) )
                                  endif
                
                                  T_E3E3(kk_s,jjh_s,kk_i,jjh_i,:,mpirank) = Te(:)
                                  if( Green_int ) T_E3E3_m(kk_s,jjh_s,kk_i,jjh_i,:,mpirank) = Te_m(:)
                
                                end do
                              end do
                            end do
                          end do
                        end do
                      end do
                
                    endif

                  end do
                end do
              end do
            end do
          end do
        end do
      end do
    end do ! end of loop over k_s

    if( icheck < 2 ) cycle

    write(3,160)
    Lm_s = 0
    do L_s = 0,max(index_rank_s,1)
      do m_s = -L_s,L_s
        Lm_s = Lm_s + 1
        Lm_i = 0
        do L_i = 0,max(index_rank_i,1)
          do m_i = -L_i,L_i
            Lm_i = Lm_i + 1
            if( sum( abs(Tens_lm(Lm_s,Lm_i,:)) ) > eps15 ) write(3,170) index_rank_s, index_rank_i, L_s, m_s, L_i, m_i, &
                                                      ( Tens_lm(Lm_s,Lm_i,initr), initr = 1,ninitr )
            if( .not. Green_int ) cycle
            if( sum( abs(Tens_lm_m(Lm_s,Lm_i,:)) ) > eps15 ) write(3,180) index_rank_s, index_rank_i, L_s, m_s, L_i, m_i, &
                 ( Tens_lm_m(Lm_s,Lm_i,initr), initr = 1,ninitr )
          end do
        end do
      end do
    end do

  end do

! Rotation to get the tensors in R1 basis

!  Rot_tem = matmul( Rot_cluster, transpose(Rot_atom_abs) )
  Rot_tem = matmul( Rot_cluster, Rot_atom_abs )

  if( abs( Rot_tem(1,1) - 1._db ) > eps10 .or. abs( Rot_tem(2,2) - 1._db ) > eps10 .or. abs( Rot_tem(3,3) - 1._db ) > eps10 ) then
  
    do initr = 1,ninitr
      if( E1E1 ) then
        mat2(:,:) = T_E1E1(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_E1E1(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E1SP ) then
        mat2(:,:) = T_E1SP(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_E1SP(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( SPSP ) then
        mat2(:,:) = T_SPSP(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_SPSP(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E1E2 ) then
        mat3(:,:,:) = T_E1E2(:,:,:,initr,mpirank)
        call Rot_tensor_3( mat3, rot_tem )
        T_E1E2(:,:,:,initr,mpirank) = mat3(:,:,:)
      endif
  
      if( E1M2 ) then
        mat3(:,:,:) = T_E1M2(:,:,:,initr,mpirank)
        call Rot_tensor_3( mat3, rot_tem )
        T_E1M2(:,:,:,initr,mpirank) = mat3(:,:,:)
      endif
  
      if( E2E2 ) then
        mat4(:,:,:,:) = T_E2E2(:,:,:,:,initr,mpirank)
        call Rot_tensor_4( mat4, rot_tem )
        T_E2E2(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E2M2 ) then
        mat4(:,:,:,:) = T_E2M2(:,:,:,:,initr,mpirank)
        call Rot_tensor_4( mat4, rot_tem )
        T_E2M2(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1E3 ) then
        mat4(:,:,:,:) = T_E1E3(:,:,:,:,initr,mpirank)
        call Rot_tensor_4( mat4, rot_tem )
        T_E1E3(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1M1 ) then
        mat2(:,:) = T_E1M1(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_E1M1(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( M1M1 ) then
        mat2(:,:) = T_M1M1(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_M1M1(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E3E3 ) then
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                Mat6(:,j_s,h_s,:,j_i,h_i) = T_E3E3(:,jh_s,:,jh_i,initr,mpirank)
              end do
            end do
          end do
        end do
        call Rot_tensor_6( Mat6, Rot_tem )
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                T_E3E3(:,jh_s,:,jh_i,initr,mpirank) = Mat6(:,j_s,h_s,:,j_i,h_i)
              end do
            end do
          end do
        end do
      endif
  
      if( .not. Green_int ) Cycle
  
      if( E1E1 ) then
        mat2(:,:) = T_E1E1_m(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_E1E1_m(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E1E2 ) then
        mat3(:,:,:) = T_E1E2_m(:,:,:,initr,mpirank)
        call Rot_tensor_3( mat3, rot_tem )
        T_E1E2_m(:,:,:,initr,mpirank) = mat3(:,:,:)
      endif
  
      if( E1M2 ) then
        mat3(:,:,:) = T_E1M2_m(:,:,:,initr,mpirank)
        call rot_tensor_3( mat3, rot_tem )
        T_E1M2_m(:,:,:,initr,mpirank) = mat3(:,:,:)
      endif
  
      if( E2E2 ) then
        mat4(:,:,:,:) = T_E2E2_m(:,:,:,:,initr,mpirank)
        call Rot_tensor_4( mat4, rot_tem )
        T_E2E2_m(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E2M2 ) then
        mat4(:,:,:,:) = T_E2M2_m(:,:,:,:,initr,mpirank)
        call Rot_tensor_4( mat4, rot_tem )
        T_E2M2_m(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1E3 ) then
        mat4(:,:,:,:) = T_E1E3_m(:,:,:,:,initr,mpirank)
        call Rot_tensor_4( mat4, rot_tem )
        T_E1E3_m(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1M1 ) then
        mat2(:,:) = T_E1M1_m(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_E1M1_m(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( M1M1 ) then
        mat2(:,:) = T_M1M1_m(:,:,initr,mpirank)
        call Rot_tensor_2( mat2, rot_tem )
        T_M1M1_m(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E3E3 ) then
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                Mat6(:,j_s,h_s,:,j_i,h_i) = T_E3E3_m(:,jh_s,:,jh_i,initr,mpirank)
              end do
            end do
          end do
        end do
        call Rot_tensor_6( Mat6, Rot_tem )
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                T_E3E3_m(:,jh_s,:,jh_i,initr,mpirank) = Mat6(:,j_s,h_s,:,j_i,h_i)
              end do
            end do
          end do
        end do
      endif
  
    end do

  endif
    
  deallocate( rof, Singul )
  
  return
  100 format(/' ---- Tenseur_car Radial ---------',100('-'))
  110 format(/' ---- Tenseur_car Tens_ab --------',100('-'))
  120 format(/' Lseuil =',i2)
  130 format(/' -- ',a4,' -- index_rank_s =',i2,', index_rank_i =',i2,' --')
  140 format(/' k_s, j_s, h_s =',3i3,',  k_i, j_i, h_i =',3i3)
  150 format(/' L_s, m_s =',2i3,',  L_i, m_i =',2i3,', clm_s, clm_i =',1p, 2e13.5)
  160 format(/' Tensor by harmonics (basis R4) :'/, ' ir jr  L_s m_s  L_i m_i      Tens(Ylm,i=1,ninitr)')
  170 format(2i3,i4,i3,2x,i4,i3,1p,20e15.7)
  180 format('i',i2,i3,i4,i3,2x,i4,i3,1p,20e15.7)

end

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

! Calculation of the componants of the transformation cartesian basis - spherical basis

subroutine Tr_Cartesian_Spherical_basis( Clm, irang1, Lmomax, Lomax, mpirank0, nrang )

  use declarations
  implicit none

  integer:: i, ipr, index_rank_s, irang1, j, k, Lmomax, Lomax, mpirank0, nrang
  
  real(kind=db):: c, c0, c1, c3, c5, c8, c12, c120
  
  real(kind=db), dimension(-1:Lomax,3,3,3,Lmomax):: Clm

  if( nrang > Lomax .and. mpirank0 == 0 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,'(/A)') ' nrang > Lomax in Tenseur_car !'
    end do
    stop
  endif
 
  do index_rank_s = irang1,nrang

    clm(index_rank_s,:,:,:,:) = 0._db

    select case(index_rank_s)

      case(-1)  
! M2
        c = sqrt( 4 * pi / 3 )
        clm(index_rank_s,:,1,:,4) = c
        clm(index_rank_s,:,2,:,2) = c
        clm(index_rank_s,:,3,:,3) = c

      case(0)
! Magnetic dipole M1
! This factor is for the integral on the solid angles. In all the other transitions including M2, this one is done in the Gaunt coefficients
        c = sqrt( 4 * pi ) 
! In that case corresponds to (L,m) = (0,0), but used to define x, y or z
        clm(index_rank_s,1,:,:,4) = c
        clm(index_rank_s,2,:,:,2) = c
        clm(index_rank_s,3,:,:,3) = c

      case(1)  
! Dipole
        c = sqrt( 4 * pi / 3 )
        clm(index_rank_s,1,:,:,4) = c
        clm(index_rank_s,2,:,:,2) = c
        clm(index_rank_s,3,:,:,3) = c

      case(2)
! Quadrupole
        c0 = sqrt( 4 * pi ) / 3
        c = sqrt( 4 * pi / 15 )
        c3 = c / sqrt( 3._db )

        clm(index_rank_s,1,1,:,1) = c0
        clm(index_rank_s,1,1,:,7) = - c3;   clm(index_rank_s,1,1,:,9) = c
        clm(index_rank_s,1,2,:,5) = c
        clm(index_rank_s,1,3,:,8) = c
        clm(index_rank_s,2,2,:,1) = c0
        clm(index_rank_s,2,2,:,7) = - c3;   clm(index_rank_s,2,2,:,9) = - c
        clm(index_rank_s,2,3,:,6) = c
        clm(index_rank_s,3,3,:,1) = c0
        clm(index_rank_s,3,3,:,7) = 2 * c3

        do i = 1,3
          do j = i+1,3
            clm(index_rank_s,j,i,:,:) = clm(index_rank_s,i,j,:,:)
          end do
        end do

      case(3)
! Octupole
        c1 = sqrt( 4 * pi / 75 )

        c = sqrt( 4 * pi / 35 )
        c5 = c / sqrt( 5._db )
        c8 = c / sqrt( 2._db )
        c12 = c / sqrt( 3._db )
        c120 = c / sqrt( 30._db )

        clm(index_rank_s,1,1,1,4) = 3 * c1
        clm(index_rank_s,1,1,1,14) = - 3 * c120; clm(index_rank_s,1,1,1,16) = c8
        clm(index_rank_s,2,2,2,2) = 3 * c1
        clm(index_rank_s,2,2,2,12) = - 3 * c120; clm(index_rank_s,2,2,2,10) = - c8
        clm(index_rank_s,3,3,3,3) = 3 * c1
        clm(index_rank_s,3,3,3,13) = 2 * c5
        clm(index_rank_s,1,1,2,2) = c1
        clm(index_rank_s,1,1,2,10) = c8;        clm(index_rank_s,1,1,2,12) = - c120
        clm(index_rank_s,1,2,2,4) = c1
        clm(index_rank_s,1,2,2,16) = - c8;      clm(index_rank_s,1,2,2,14) = - c120
        clm(index_rank_s,1,1,3,3) = c1
        clm(index_rank_s,1,1,3,13) = - c5;      clm(index_rank_s,1,1,3,15) = c12
        clm(index_rank_s,2,2,3,3) = c1
        clm(index_rank_s,2,2,3,13) = - c5;      clm(index_rank_s,2,2,3,15) = - c12
        clm(index_rank_s,2,3,3,2) = c1
        clm(index_rank_s,2,3,3,12) = 4 * c120
        clm(index_rank_s,1,3,3,4) = c1
        clm(index_rank_s,1,3,3,14) = 4 * c120
        clm(index_rank_s,1,2,3,11) = c12

        do i = 1,3
          do j = i,3
            do k = j,3
              if( i == j .and. i == k ) cycle
              clm(index_rank_s,i,k,j,:) = clm(index_rank_s,i,j,k,:)
              clm(index_rank_s,j,i,k,:) = clm(index_rank_s,i,j,k,:)
              clm(index_rank_s,k,i,j,:) = clm(index_rank_s,i,j,k,:)
              clm(index_rank_s,j,k,i,:) = clm(index_rank_s,i,j,k,:)
              clm(index_rank_s,k,j,i,:) = clm(index_rank_s,i,j,k,:)
            end do
          end do
        end do

    end select

  end do

  return
end

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

subroutine Tens_ab( Chi, Coef_g, Core_resolved, Final_tddft, Full_potential, &
                                  Green, Green_int, icheck, index_rank_i, index_rank_s, ip_max, ip0, is_g, k_i, k_s, L_s, &
                                  m_k_i, m_k_s, m_s, L_i, m_g, m_i, Lmax, Lmoins1, Lplus1, Lseuil, nbseuil, ns_dipmag, ninit1, &
                                  ninit, ninit_tddft, ninitv, ninitr, nlm_probe, nlm_p_fp, NRIXS, nspinp, nspino, Omega, rof, &
                                  Singul, Solsing, Solsing_only, Spinorbite, SpinPos_i, SpinPos_s, Taull, Ten, Ten_m, Ylm_comp )

  use declarations
  implicit none

  integer, intent(in):: icheck, index_rank_i, index_rank_s, ip_max, ip0, L_i, L_s, m_i, m_s, Lseuil, ninit1, &
      ninit, ninit_tddft, ninitv, ninitr, nlm_probe, nlm_p_fp, ns_dipmag, nspinp, nspino
  integer, dimension(ninit,2), intent(in):: m_g
  integer, dimension(ninit), intent(in):: is_g

  integer:: i_g_1, i_g_2, initl1, initl2, initr, is_dipmag, is_r1, is_r2, iseuil1, iseuil2, iso1, iso2, &
    ispf1, ispf2, ispinf1, ispinf2, isping1, isping2, ispp_f1, ispp_f2, k_i, k_s, &
    L_f1, L_f2, L_g, L_o_i, L_o_s, Lm01, Lm02, Lm_f1, Lm_f2, Lmax, lmp01, lmp02, lmp_f1, &
    lmp_f2, Lmps_f1, Lmps_f2, Lmv1, Lmv2, lp_f1, lp_f2, m_f1, m_f2, m_k_i, m_k_s, mi1, mi2, mp_f1, mp_f2, mv1, mv2, &
    nbseuil

  complex(kind=db):: cfe, cfs, Cg, Fac_SpinPos, Fac_SpinPos_i, Fac_SpinPos_s, Gaunt_i, Gaunt_nrixs, Gaunt_s,&
                      Gaunt_xrc, PseudoGaunt_M1, PseudoGaunt_M2, rof_1, rof_2, Tau_rad, Tau_rad_i

  complex(kind=db), dimension(ninitr):: Ten, Ten_m
  complex(kind=db), dimension(nlm_probe,nlm_p_fp,nspinp,nspino,ip0:ip_max,ninitv):: rof
  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp):: Taull
  complex(kind=db), dimension(nlm_probe*nspino,nlm_probe*nspino,ninit_tddft,ninit_tddft,2,2,ns_dipmag):: Chi
  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp,ip0:ip_max,ip0:ip_max,ninitv):: Singul

  logical:: Core_resolved, Final_tddft, Full_potential, Green, Green_int, Lmoins1, Lplus1, NRIXS, &
            Solsing, Solsing_only, Spinorbite, SpinPos_i, SpinPos_s, Titre, Ylm_comp

  real(kind=db):: Ci_1, Ci_2, Ci2, J_initl1, J_initl2, Jz1, Jz2

  real(kind=db), dimension(nbseuil):: Omega
  real(kind=db), dimension(ninit,2):: Coef_g
  
  L_g = Lseuil
  L_o_i = abs( index_rank_i )
  L_o_s = abs( index_rank_s )

  Ten(:) = (0._db,0._db)
  Ten_m(:) = (0._db,0._db)

! Index "1" and "s" are for outgoing, that is left side matrix
!       "2" and "i"         incoming,         right

! Loop over core states < g_1 ...  > < ... >
  do i_g_1 = 1,ninit

    if( Final_tddft ) then
      initl1 = i_g_1
    else
      initl1 = 1
    endif

    J_initl1 = L_g + 0.5_db * is_g(i_g_1)
    if( i_g_1 <= ninit1 ) then
      Jz1 = - J_initl1 + i_g_1 - 1
    else
      Jz1 = - J_initl1 + i_g_1 - ninit1 - 1
    endif

    if( i_g_1 <= ninit1 ) then
      iseuil1 = 1
    else
      iseuil1 = 2
    endif

    if( Final_tddft ) then
      initr = 1
    elseif( Core_resolved ) then
      initr = i_g_1
    else
      initr = iseuil1
    endif

    if( ninitv == ninit .and. Final_tddft ) then
      is_r1 = i_g_1
    else
      is_r1 = iseuil1
    endif

! Loop over core state spin < g_1 ...  > < ... >
    do isping1 = 1,2 
    
      mi1 = m_g(i_g_1,isping1)
      Ci_1 = Coef_g(i_g_1,isping1)

      if( abs( Ci_1 ) < eps6 ) cycle

      do i_g_2 = 1,ninit

        if( .not. Final_tddft .and. i_g_1 /= i_g_2 ) cycle

        if( Final_tddft ) then
          initl2 = i_g_2
        else
          initl2 = 1
        endif

        J_initl2 = L_g + 0.5_db * is_g(i_g_2)
        if( i_g_2 <= ninit1 ) then
          Jz2 = - J_initl2 + i_g_2 - 1
        else
          Jz2 = - J_initl2 + i_g_2 - ninit1 - 1
        endif

        if( i_g_2 <= ninit1 ) then
          iseuil2 = 1
        else
          iseuil2 = 2
        endif

        if( ninitv == ninit .and. Final_tddft ) then
          is_r2 = i_g_2
        else
          is_r2 = iseuil2
        endif

! Loop over core state spin < ...  > < ... g_2 >
        do isping2 = 1,2  
        
!          if( .not. Final_tddft .and. isping1 /= isping2 ) cycle
          if( .not. ( Final_tddft .or. L_g /= 0 ) .and. isping1 /= isping2 ) cycle

          mi2 = m_g(i_g_2,isping2)
          Ci_2 = Coef_g(i_g_2,isping2)

          if( abs( Ci_2 ) < eps6 ) cycle

          Ci2 = Ci_1 * Ci_2

         if( icheck > 1 ) then
           write(3,'(/A)') ' ini1 isg1 Jz1 mi1  Ci1  ini2 isg2 Jz2 mi2  Ci2'
           write(3,120) i_g_1, isping1, Jz1, mi1, Ci_1, i_g_2, isping2, Jz2, mi2, Ci_2
         endif
         Titre = .true.

! Loop over spherical harmonics of final states < ... f_1 > < ... >
          do L_f1 = 0,Lmax

            if( Lplus1 .and. L_f1 < L_g + 1 ) cycle
            if( Lmoins1 .and. L_f1 > L_g ) cycle
            if( L_f1 > L_g + L_o_s .or. L_f1 < L_g - L_o_s .or. mod(L_f1,2) /= mod( L_g + L_o_s, 2 ) ) cycle

            Lm01 = L_f1**2 + L_f1 + 1
            do m_f1 = -L_f1,L_f1
              Lm_f1 = Lm01 + m_f1
              if( Lm_f1 > nlm_probe ) cycle

! Loop over spin of final states < ... f_1 > < ... >
              do ispinf1 = 1,2

                if( SpinPos_s ) then
                  Fac_SpinPos_s = Fac_SpinPos( ispinf1, isping1, k_s, m_s, 1 )
                  if( abs( Fac_SpinPos_s ) < eps10 ) cycle
                  Fac_SpinPos_s = img * Omega(iseuil1) * Fac_SpinPos_s
                else
! Spin-flip only for magnetic dipole or spin position
                  if( ispinf1 /= isping1 .and. ( NRIXS .or. index_rank_s > 0 ) ) cycle
                endif
                  
                ispf1 = min( ispinf1, nspinp )

                if( NRIXS ) then
                  Gaunt_s = Gaunt_nrixs(L_f1,m_f1,L_s,m_s,L_g,mi1,Ylm_comp)
                elseif( index_rank_s == 0 ) then
                  Gaunt_s = PseudoGaunt_M1(L_f1,m_f1,ispinf1,m_s,L_g,mi1,isping1,Ylm_comp,.true.)
                elseif( index_rank_s == -1 ) then
                  Gaunt_s = PseudoGaunt_M2(L_f1,m_f1,ispinf1,m_s,m_k_s,L_g,mi1,isping1,Ylm_comp,.true.)
                else
                  Gaunt_s = Gaunt_xrc(L_f1,m_f1,L_s,m_s,L_g,mi1,Ylm_comp)
                endif
                if( SpinPos_s ) Gaunt_s = Fac_SpinPos_s * Gaunt_s
                  
                if( abs(Gaunt_s) < eps10 ) cycle

! Loop over spin of final states < ... > < f_2 ... >
                do L_f2 = 0,Lmax
                  if( Lplus1 .and. L_f2 < L_g + 1 ) cycle
                  if( Lmoins1 .and. L_f2 > L_g ) cycle
                  if( L_f2 > L_g + L_o_i .or. L_f2 < L_g - L_o_i .or. mod(L_f2,2) /= mod( L_g + L_o_i, 2 ) ) cycle

                  Lm02 = L_f2**2 + L_f2 + 1
                  do m_f2 = -L_f2,L_f2
                    Lm_f2 = Lm02 + m_f2
                    if( Lm_f2 > nlm_probe ) cycle

                    do ispinf2 = 1,2  ! spin final state incoming

                      if( SpinPos_i ) then
                        Fac_SpinPos_i = Fac_SpinPos( ispinf2, isping2, k_i, m_i, 2 )
                        if( abs( Fac_SpinPos_i ) < eps10 ) cycle
                        Fac_SpinPos_i = img * Omega(iseuil2) * Fac_SpinPos_i 
                      else
! Spin-flip only for M1, M2 or spin position
                        if( ispinf2 /= isping2 .and. ( NRIXS .or. index_rank_i > 0 ) ) cycle
                      endif

                      ispf2 = min( ispinf2, nspinp )

                      if( .not. ( Final_tddft .or. Spinorbite ) .and. ispinf2 /= ispinf1 ) cycle

                      if( NRIXS ) then
                        Gaunt_i = Gaunt_nrixs(L_f2,m_f2,L_i,m_i,L_g,mi2,Ylm_comp)
                      elseif( index_rank_i == 0 ) then
                        Gaunt_i = PseudoGaunt_M1(L_f2,m_f2,ispinf2,m_i,L_g,mi2,isping2,Ylm_comp,.true.)
                      elseif( index_rank_i == -1 ) then
                        Gaunt_i = PseudoGaunt_M2(L_f2,m_f2,ispinf2,m_i,m_k_i,L_g,mi2,isping2,Ylm_comp,.true.)
                      else
                        Gaunt_i = Gaunt_xrc(L_f2,m_f2,L_i,m_i,L_g,mi2,Ylm_comp)
                      endif

                      if( abs(Gaunt_i) < eps10 ) cycle
                      if( SpinPos_i ) Gaunt_i = Fac_SpinPos_i * Gaunt_i

                      Cg = Ci2 * conjg( Gaunt_s ) * Gaunt_i

                      if( Final_tddft .and. ispinf1 /= isping1 ) then
                        is_dipmag = 2
                      else
                        is_dipmag = 1
                      endif

                      if( .not. Solsing_only ) then
                          
! Loop over harmonics of non spherical final states < ... f_1 > < ... >
          do lp_f1 = 0,Lmax
            if( .not. Full_potential .and. lp_f1 /= L_f1 ) cycle
            lmp01 = lp_f1**2 + lp_f1 + 1

            do mp_f1 = -lp_f1,lp_f1
              if( nlm_p_fp == 1 .and. mp_f1 /= m_f1 ) cycle
              lmp_f1 = min( lmp01 + mp_f1, nlm_p_fp ) 

! Loop over spin-orbit solution index < ... f_1 > < ... >
              do ispp_f1 = 1,nspinp
                if( .not. Spinorbite .and. ispp_f1 /= ispf1 ) cycle
                iso1 = min( ispp_f1, nspino )

                if( Spinorbite .and. nlm_p_fp == 1 ) then
                  mv1 = mp_f1 - ispinf1 + iso1
                  if( mv1 > L_f1 .or. mv1 < -L_f1 ) cycle
                else
                  mv1 = mp_f1
                endif
                Lmv1 = Lmp01 + mv1 
                Lmps_f1 = ( lmp01 + mv1 - 1 ) * nspino + iso1

! Loop over harmonics of non spherical final states < ... > < f_2 ... >
                do lp_f2 = 0,Lmax
                  if( .not. Full_potential .and. lp_f2 /= L_f2 ) cycle
                  lmp02 = lp_f2**2 + lp_f2 + 1

                  do mp_f2 = -lp_f2,lp_f2
                    if( nlm_p_fp == 1 .and. mp_f2 /= m_f2 ) cycle
                    lmp_f2 = min( lmp02 + mp_f2, nlm_p_fp ) 

! Loop over spin-orbit solution index < ... > < f_2 ... >
                    do ispp_f2 = 1,nspinp
                      if( .not. Spinorbite .and. ispp_f2 /= ispf2 ) cycle
                      iso2 = min( ispp_f2, nspino )

                      if( Spinorbite .and. nlm_p_fp == 1 ) then
                        mv2 = mp_f2 - ispinf2 + iso2
                        if( mv2 > L_f2 .or. mv2 < -L_f2 ) cycle
                      else
                        mv2 = mp_f2
                      endif
                      Lmv2 = Lmp02 + mv2 
                      Lmps_f2 = ( lmp02 + mv2 - 1 ) * nspino + iso2

                      cfe = ( 0._db, 0._db ) 
                      cfs = ( 0._db, 0._db ) 
                                              
 ! Lmps_f index the states, Lm_f index the harmonics 
                      rof_1 = rof(Lm_f1,lmp_f1,ispf1,iso1,L_o_s,is_r1)
                      rof_2 = rof(Lm_f2,lmp_f2,ispf2,iso2,L_o_i,is_r2)

                      if( Final_TDDFT ) then
                        if( Green ) then
                          cfe = cfe + rof_1 * rof_2 * Chi(Lmps_f1,Lmps_f2,initl1,initl2,ispinf1,ispinf2,is_dipmag)
                          cfs = cfs + rof_2 * rof_1 * Chi(Lmps_f2,Lmps_f1,initl2,initl1,ispinf2,ispinf1,is_dipmag)
                        else
                          cfe = cfe + conjg( rof_1 ) * rof_2 * Chi(Lmps_f1,Lmps_f2,initl1,initl2,ispinf1,ispinf2,is_dipmag)
                          cfs = cfs + conjg( rof_2 ) * rof_1 * Chi(Lmps_f2,Lmps_f1,initl2,initl1,ispinf2,ispinf1,is_dipmag)
                        endif
                        
                      else
                      
                        if( Green ) then
                          cfe = cfe + rof_1 * rof_2 * Taull(Lmv1,ispp_f1,Lmv2,ispp_f2)
                          cfs = cfs + rof_2 * rof_1 * Taull(Lmv2,ispp_f2,Lmv1,ispp_f1)
                        else
                          cfe = cfe + conjg( rof_1 ) * rof_2 * Taull(Lmv1,ispp_f1,Lmv2,ispp_f2)
                          cfs = cfs + conjg( rof_2 ) * rof_1 * Taull(Lmv2,ispp_f2,Lmv1,ispp_f1)
                        endif

                      endif

! Here one does not divide by pi, so convenient results seems multiplied by pi, when calculating atomic form factor.
! So, the normalization by pi must not be done in coabs to get the atomic form factor in DAFS.
                      if( Green_int ) then
! In this case all tensorial terms are complex, thus one separates the magnetic part explicitely.
! by this way furthur (in coabs) when time reversal apply, the corresponding term change of sign.  
                        Tau_rad = 0.5_db * ( cfe + cfs )
                        Tau_rad_i = 0.5_db * (cfe - cfs)
                        Ten(initr) = Ten(initr) + Real(Cg,db) * Tau_rad + img * aimag(Cg)* Tau_rad_i
                        Ten_m(initr) = Ten_m(initr) + Real(Cg,db) * Tau_rad_i + img * aimag(Cg) * Tau_rad
                      else
                        Tau_rad = 0.5_db * ( cfe - conjg( cfs ) )
                        Ten(initr) = Ten(initr) + Cg * Tau_rad
                      endif

                     if( icheck > 1 .and. abs( Tau_rad ) > eps15 ) then
                        if( Titre ) then
                          Titre = .false.
                          if( nlm_p_fp == 1 ) then
                            write(3,130)
                          else
                            write(3,135)
                          endif
                        endif
                        if( nlm_p_fp == 1 ) then
                          write(3,140) L_f1, m_f1, ispinf1, iso1, L_f2, m_f2, ispinf2, iso2, &
                            Ten(initr), Cg*Tau_rad, cfe, cfs, Gaunt_i, Gaunt_s, &
                            rof_1, rof_2, Taull(Lmv2,ispp_f2,Lmv1,ispp_f1), Taull(Lmv1,ispp_f1,Lmv2,ispp_f2)
                        else
                          write(3,145) L_f1, m_f1, lp_f1, mp_f1, ispinf1, iso1, L_f2, m_f2, lp_f2, mp_f2, ispinf2, iso2, &
                            Ten(initr), Cg*Tau_rad, cfe, cfs, Gaunt_i, Gaunt_s, &
                            rof_1, rof_2, &
                            Chi(Lmps_f2,Lmps_f1,initl2,initl1,ispp_f2,ispp_f1,is_dipmag), &
                            Chi(Lmps_f1,Lmps_f2,initl1,initl2,ispp_f1,ispp_f2,is_dipmag)
                        endif 
                      endif

                    end do ! end of loop over spin-orbit solution f_i ( incoming )
                  end do ! end of loop over mp_f2 non-spherical
                end do ! end of loop over lp_f2 non-spherical

              end do ! end of loop over spin-orbit solution f_s ( outgoing )
            end do ! end of loop over mp_f1 non-spherical
          end do ! end of loop over lp_f1 non-spherical

                      endif
                      
                      if( Solsing .and. i_g_1 == i_g_2 ) then
                        Ten(initr) = Ten(initr) + Cg * Singul(Lm_f1,ispf1,Lm_f2,ispf2,L_o_s,L_o_i,is_r1)

                        if( icheck > 1 .and. abs( Singul(Lm_f1,ispf1,Lm_f2,ispf2,index_rank_s,L_o_i,is_r1) ) > eps15 ) &
                          write(3,140) L_f1, m_f1, ispinf1, 0, L_f2, m_f2, ispinf2, 0, &
                            Ten(initr), Cg*Singul(Lm_f1,ispf1,Lm_f2,ispf2,index_rank_s,L_o_i,is_r1)

                      endif

                    end do ! end of loop over ispinf2, spin of final state ( incoming )
                  end do !  end of loop over  m_f2 final state
                end do !  end of loop over  L_f2 final state

              end do  !  end of loop over spin of final state ( outgoing )
            end do  !  end of loop over  m_f1 final state
          end do !  end of loop over  L_f1 final state

        end do    !  end of loop over spin of core states ( incoming )
      end do    !  end of loop over core states ( incoming )
    end do    !  end of loop over spin of core states ( outgoing )
  end do   !  end of loop over core states  ( outgoing )

! One multiplies by "i" in order the real part of the tensor is the absorption.
  Ten(:) = img * Ten(:)
  Ten_m(:) = img * Ten_m(:)

  return
  120 format(2(2i4,f6.1,i3,f7.3))
  130 format('  l1  m1 is1 io1  l2  m2 is2 io2',15x,'Ten',27x,'dTen',28x,'cfe',27x,'cfs',27x, &
                'Gaunt_i',24x,'Gaunt_s',25x,'rof_s',26x,'rof_i',26x,'Tau_s',26x,'Tau_i')
  135 format('  l1  m1 lp1 mp1 is1 io1  l2  m2 lp2 mp2 is2 io2',15x,'Ten',27x,'dTen',28x,'cfe',27x,'cfs',27x, &
                'Gaunt_i',24x,'Gaunt_s',25x,'rof_s',26x,'rof_i',26x,'Tau_s',26x,'Tau_i')
  140 format(8i4,1p,16(1x,2e15.7))
  145 format(12i4,1p,16(1x,2e15.7))
end

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

function Fac_SpinPos( ispinf, isping, k, m, is )

  use declarations
  implicit none

  integer:: is, ispinf, isping, k, m

  complex(kind=db):: Fac_SpinPos

  if( k == 1 ) then   ! along epsilon_x
    
    if( m == -1 ) then   ! y
      if( ispinf == isping ) then
        if( ispinf == 1 ) then
          Fac_SpinPos = ( 1._db, 0._db )
        else
          Fac_SpinPos = ( - 1._db, 0._db )
        endif
      else
        Fac_SpinPos = ( 0._db, 0._db )
      endif

    elseif( m == 0 ) then   ! z

      if( ispinf == isping ) then
        Fac_SpinPos = ( 0._db, 0._db )
      elseif( ispinf == 1 ) then
        Fac_SpinPos = ( 0._db, 1._db )
      else
        Fac_SpinPos = ( 0._db, - 1._db )
      endif
    
    else
      
      write(6,'(/A)') ' Error in Fac_SpinPos !'
      if( is == 1 ) then
        write(6,'(a6,i3)') ' k_s =', k
      else
        write(6,'(a6,i3)') ' k_i =', k
      endif
      write(6,'(a17,2i3)') ' ispinf, isping =', ispinf, isping   
      if( is == 1 ) then
        write(6,'(a6,i3)') ' m_s =', m
      else
        write(6,'(a6,i3)') ' m_i =', m
      endif   
      stop ! one cannot be here
   
    endif
   
  elseif( k == 2 ) then  ! along epsilon_y
    
    if( m == 1 ) then   ! x
    
      if( ispinf == isping ) then
        if( ispinf == 1 ) then
          Fac_SpinPos = ( - 1._db, 0._db )
        else
          Fac_SpinPos = ( 1._db, 0._db )
        endif
      else
        Fac_SpinPos = ( 0._db, 0._db )
      endif
    
    elseif( m == 0 ) then   ! z

      if( ispinf /= isping ) then
        Fac_SpinPos = ( 1._db, 0._db )
      else
        Fac_SpinPos = ( 0._db, 0._db )
      endif
    
    else
    
      write(6,'(/A)') ' Error in Fac_SpinPos !'
      if( is == 1 ) then
        write(6,'(a6,i3)') ' k_s =', k
      else
        write(6,'(a6,i3)') ' k_i =', k
      endif
      write(6,'(a17,2i3)') ' ispinf, isping =', ispinf, isping   
      if( is == 1 ) then
        write(6,'(a6,i3)') ' m_s =', m
      else
        write(6,'(a6,i3)') ' m_i =', m
      endif   
      stop ! one cannot be here
       
    endif
   
  else ! k == 3, along z
    
    if( ispinf == isping ) then
    
      Fac_SpinPos = ( 0._db, 0._db )
    
    elseif( m == 1 ) then   ! x
    
      if( ispinf == 1 ) then
        Fac_SpinPos = ( 0._db, - 1._db ) 
      else
        Fac_SpinPos = ( 0._db, 1._db )
      endif
      
    elseif( m == - 1 ) then   ! y
    
      Fac_SpinPos = ( - 1._db, 0._db )
      
    else
  
      write(6,'(/A)') ' Error in Fac_SpinPos !'
      if( is == 1 ) then
        write(6,'(a6,i3)') ' k_s =', k
      else
        write(6,'(a6,i3)') ' k_i =', k
      endif
      write(6,'(a17,2i3)') ' ispinf, isping =', ispinf, isping   
      if( is == 1 ) then
        write(6,'(a6,i3)') ' m_s =', m
      else
        write(6,'(a6,i3)') ' m_i =', m
      endif   
      stop ! one cannot be here
     
    endif
  
  endif

  return
end

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

subroutine Tenseur_car_optic( Chi, Core_resolved, Ecinetic, Eimag, Energ, Enervide, Final_tddft, &
                  Full_potential, Green_int, Hubb_a, Hubb_d, icheck, ip_max, ip0, Lmax, Lmax_pot, &
                  Ldip, Lmoins1, Loct, Lplus1, Lqua, m_hubb, mpinodes, mpirank, mpirank0, &
                  msymdd, msymddi, msymdq, msymdqi, msymdo, msymdoi, msymoo, msymooi, msymqq, msymqqi, Multipole, &
                  n_Ec, n_oo, n_V, ns_dipmag, ninitr, nlm_pot, nlm_probe, &
                  nlm_p_fp, nr, nspin, nspino, nspinp, Numat, r, Relativiste, Renorm, Rmtg, Rmtsd, Rot_atom_abs, &
                  Rot_cluster, T_E1E1, T_E1E1_m, T_E1E3, T_E1E3_m, T_E1E2, T_E1E2_m, T_E1SP, T_E1SP_m, &
                  T_E1M1, T_E1M1_m, T_M1M1, T_M1M1_m, T_E3E3, T_E3E3_m, T_E2E2, T_E2E2_m, T_SPSP, T_SPSP_m, &
                  Solsing, Spinorbite, Taull, V_hubb, V_intmax, V0bd, Vrato, Ylm_comp )

  use declarations
  implicit none

  integer, parameter:: Lomax = 3
  integer, parameter:: Lmomax = ( Lomax + 1 )**2

  integer:: h_s, hh_s, hh_i, h_i, i, icheck, ief, index_rank_i, index_rank_s, initr, ip_max, ip0, irang1, is, isi, &
    j, j_i, j_s, jh_i, jh_s, jj_i, jj_s, jjh_i, jjh_s, k_i, k_s, kk_i, kk_s, &
    L_i, L_s, Lm_i, Lm_s, Lmax, Lmax_pot, m_hubb, m_i, m_s, mpinodes, mpirank, mpirank0, &
    n_Ec, n_oo, n_V, ninitr, nh_i, nh_s, nj_i, nj_s, nlm_p_fp, nlm_pot, nlm_probe, nlm1g, &
    nlm2g, nr, nr_zet, nrang, ns_dipmag, nspin, nspino, nspinp, Numat

  integer, dimension(3):: Ldip
  integer, dimension(3,3):: Lqua, msymdd, msymddi
  integer, dimension(3,3,3):: Loct, msymdq, msymdqi
  integer, dimension(3,3,3,3):: msymdo, msymdoi, msymqq, msymqqi
  integer, dimension(3,n_oo,3,n_oo):: msymoo, msymooi

  complex(kind=db), dimension(ninitr):: Te, Te_m, Ten, Ten_m, Tens, Tens_m
  complex(kind=db), dimension(3,3) :: Mat2
  complex(kind=db), dimension(3,3,3) :: Mat3
  complex(kind=db), dimension(3,3,3,3) :: Mat4
  complex(kind=db), dimension(3,3,3,3,3,3) :: Mat6
  complex(kind=db), dimension(3,3,ninitr,0:mpinodes-1):: T_E1E1, T_E1E1_m, T_E1SP, T_E1SP_m, T_E1M1, T_E1M1_m, T_M1M1, T_M1M1_m, &
                                                         T_SPSP, T_SPSP_m 
  complex(kind=db), dimension(3,3,3,ninitr,0:mpinodes-1):: T_E1E2, T_E1E2_m
  complex(kind=db), dimension(3,3,3,3,ninitr,0:mpinodes-1):: T_E1E3, T_E1E3_m, T_E2E2, T_E2E2_m
  complex(kind=db), dimension(3,n_oo,3,n_oo,ninitr,0:mpinodes-1):: T_E3E3, T_E3E3_m
  complex(kind=db), dimension(lmomax,lmomax,ninitr):: Tens_lm, Tens_lm_m
  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp,2):: Taull
  complex(kind=db), dimension(-m_hubb:m_hubb,nspinp,-m_hubb:m_hubb,nspinp):: V_hubb
  complex(kind=db), dimension(nlm_probe*nspino,nlm_probe*nspino,nlm_probe*nspino,nlm_probe*nspino,2,2,ns_dipmag):: Chi

  logical:: Core_resolved, E1E1, E1E2, E1E3, E1M1, E1SP, E2E2, E3E3, &
    Final_tddft, Full_potential, Green_int, Hubb_a, Hubb_d, Lmoins1, Lplus1, M_depend, M1M1, No_diag, NRIXS, Relativiste, &
    Renorm, Solsing, Spinorbite, SPSP, Ylm_comp

  logical, dimension(13):: Multipole

  real(kind=db):: clm_s, clm_i, Eimag, Energ, Rmtg, Rmtsd, V_intmax
  real(kind=db), dimension(1):: Energ_t
  real(kind=db), dimension(ninitr):: Tensi, Tensr
  real(kind=db), dimension(-1:lomax,3,3,3,lmomax):: Clm
  real(kind=db), dimension(3,3):: Rot_atom_abs, Rot_cluster, Rot_tem
  real(kind=db), dimension(nspin,n_Ec):: Ecinetic
  real(kind=db), dimension(n_Ec):: Eimag_t, Enervide
  real(kind=db), dimension(nspin,n_V):: V0bd
  real(kind=db), dimension(nr):: r
  real(kind=db), dimension(nr,nlm_pot,nspin,n_V):: Vrato

  real(kind=db), dimension(:,:,:,:,:,:), allocatable:: zet
  real(kind=db), dimension(:,:,:,:,:,:,:), allocatable:: roff_rr, roff0

  if( icheck > 1 ) write(3,100)

  E1E1 = Multipole(1); E1E2 = Multipole(2); E1E3 = Multipole(3);
  E1M1 = Multipole(4); E2E2 = Multipole(6);
  E3E3 = Multipole(7); M1M1 = Multipole(8)
  E1SP = Multipole(11); SPSP = Multipole(12)

  NRIXS = .false.

! Calculation of radial integrals for regular and irregular solutions

  No_diag = Full_potential .or. ( Hubb_a .and. .not. Hubb_d )
  M_depend = Spinorbite .or. Hubb_a .or. No_diag

  if( No_diag ) then
    nlm1g = ( Lmax + 1 )**2
    nlm2g = nlm1g
  elseif( M_depend ) then
    nlm1g = ( Lmax + 1 )**2
    nlm2g = 1
  else
    nlm1g = Lmax + 1
    nlm2g = 1
  endif
  allocate( roff_rr(nlm1g,nlm1g,nlm2g,nlm2g,nspinp**2,nspino**2,ip0:ip_max) )
  nr_zet = 0
  allocate( zet(nr_zet,nlm1g,nlm2g,nspinp,nspino,n_Ec) )
  ief = 2
  allocate( roff0(ief-1,ief:n_Ec,nlm1g,nlm1g,nspinp,nspinp,nspino**2) )
  Eimag_t(:) = Eimag
  Energ_t(1) = Energ

  call Radial_optic(Ecinetic,Eimag_t,Energ_t,Enervide,Full_potential,Hubb_a,Hubb_d,icheck,ief,ip_max,ip0, &
       Lmax,Lmax_pot,m_depend,m_hubb,1,n_Ec,n_V,nlm_pot,nlm1g,nlm2g,No_diag,nr,nr_zet,nspin,nspino,nspinp,Numat,r,Relativiste, &
       Renorm,Rmtg,Rmtsd,roff_rr,roff0,Solsing,Spinorbite,V_hubb,V_intmax,V0bd,Vrato,Ylm_comp,zet)
  
  deallocate( roff0, zet )

  if( icheck > 1 ) write(3,110)

  if( E1E3 .or. E3E3 ) then
    nrang = 3
  elseif( E1E2 .or. E2E2 ) then
    nrang = 2
  else
    nrang = 1
  endif
  if( E1M1 .or. M1M1 ) then
    irang1 = 0
  else
    irang1 = 1
  endif

! Calculation of the componants of the transformation cartesian basis - spherical basis
  call Tr_Cartesian_Spherical_basis( Clm, irang1, Lmomax, Lomax, mpirank0, nrang )

  T_M1M1(:,:,:,mpirank) = (0._db,0._db)
  T_E1M1(:,:,:,mpirank) = (0._db,0._db)
  T_E1E2(:,:,:,:,mpirank) = (0._db,0._db)
  T_E2E2(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_E1E3(:,:,:,:,:,mpirank) = (0._db,0._db)
  T_E1E1(:,:,:,mpirank) = (0._db,0._db)
  T_E1SP(:,:,:,mpirank) = (0._db,0._db)
  T_SPSP(:,:,:,mpirank) = (0._db,0._db)

  if( E3E3 ) T_E3E3(:,:,:,:,:,mpirank) = (0._db,0._db)

  if( Green_int ) then
    T_M1M1_m(:,:,:,mpirank) = (0._db,0._db)
    T_E1M1_m(:,:,:,mpirank) = (0._db,0._db)
    T_E1E2_m(:,:,:,:,mpirank) = (0._db,0._db)
    T_E2E2_m(:,:,:,:,:,mpirank) = (0._db,0._db)
    T_E1E3_m(:,:,:,:,:,mpirank) = (0._db,0._db)
    T_E1E1_m(:,:,:,mpirank) = (0._db,0._db)
    T_E1SP_m(:,:,:,mpirank) = (0._db,0._db)
    T_SPSP_m(:,:,:,mpirank) = (0._db,0._db)
    if( E3E3 ) T_E3E3_m(:,:,:,:,:,mpirank) = (0._db,0._db)
  endif

! Loops over tensor rank
  do index_rank_s = irang1,nrang
    do index_rank_i = index_rank_s,nrang
      Tens_lm(:,:,:) = (0._db,0._db)
      Tens_lm_m(:,:,:) = (0._db,0._db)
      if( .not. E1E1 .and. index_rank_s == 1 .and. index_rank_i == 1 ) cycle
      if( .not. E1E2 .and. index_rank_s == 1 .and. index_rank_i == 2 ) cycle
      if( .not. E2E2 .and. index_rank_s == 2 .and. index_rank_i == 2 ) cycle
      if( .not. E1E3 .and. index_rank_s == 1 .and. index_rank_i == 3 ) cycle
      if( .not. M1M1 .and. index_rank_s == 0 .and. index_rank_i == 0 ) cycle
      if( .not. E1M1 .and. index_rank_s == 0 .and. index_rank_i == 1 ) cycle
      if( .not. E3E3 .and. index_rank_s == 3 .and. index_rank_i == 3 ) cycle
      if( index_rank_s == 0 .and. index_rank_i > 1 ) cycle
      if( index_rank_i == 3 .and. index_rank_s == 2 ) cycle

      if( icheck > 1 ) write(3,130) index_rank_s, index_rank_i

! Loops over tensor indices
      do k_s = 1,3

        if( index_rank_s == 1 .and. Ldip(k_s) /= 1 ) cycle

        if( index_rank_s > 1 ) then
          nj_s = 3
        else
          nj_s = 1
        endif

        do j_s = 1,nj_s

          if( index_rank_s == 2 .and. Lqua(k_s,j_s) /= 1 ) cycle

          if( index_rank_s == 3 ) then
            nh_s = 3
          else
            nh_s = 1
          endif

          do h_s = 1,nh_s

            if( index_rank_s == 3 .and. Loct(k_s,j_s,h_s) /= 1 ) cycle
            jh_s = 3 * ( j_s - 1 ) + h_s

            do k_i = 1,3

              if( index_rank_i == 1 .and. Ldip(k_i) /= 1 ) cycle

              if( index_rank_i > 1 ) then
                nj_i = 3
              else
                nj_i = 1
              endif

              do j_i = 1,nj_i

                if( index_rank_i == 2 .and. Lqua(k_i,j_i) /= 1 ) cycle

                if( index_rank_i == 3 ) then
                  nh_i = 3
                else
                  nh_i = 1
                endif

                do h_i = 1,nh_i

                  if( index_rank_i == 3 .and. Loct(k_i,j_i,h_i) /= 1 ) cycle
                  jh_i = 3 * ( j_i - 1 ) + h_i

                  if( index_rank_s == 0 .and. index_rank_i == 0 ) then
                    if( k_i < k_s ) cycle
                  elseif( index_rank_s == 1 .and. index_rank_i == 1 ) then
                    if( k_i < k_s ) cycle
                    if( msymdd(k_s,k_i) == 0 .and. msymddi(k_s,k_i) == 0 ) cycle
                    if( sum( abs( T_E1E1(k_s,k_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                  elseif( index_rank_s == 1 .and. index_rank_i == 2 ) then
                    if( msymdq(k_s,k_i,j_i) == 0 .and. msymdqi(k_s,k_i,j_i) == 0 ) cycle
                    if( sum( abs( T_E1E2(k_s,k_i,j_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                  elseif( index_rank_s == 2 .and. index_rank_i == 2 ) then
                    if( msymqq(k_s,j_s,k_i,j_i) == 0 .and. msymqqi(k_s,j_s,k_i,j_i) == 0 ) cycle
                    if( sum( abs( T_E2E2(k_s,j_s,k_i,j_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                  elseif( index_rank_s == 1 .and. index_rank_i == 3 ) then
                    if( msymdo(k_s,k_i,j_i,h_i) == 0 .and. msymdoi(k_s,k_i,j_i,h_i) == 0 ) cycle
                  elseif( index_rank_s == 3 .and. index_rank_i == 3 ) then
                    if( msymoo(k_s,jh_s,k_i,jh_i) == 0 .and. msymooi(k_s,jh_s,k_i,jh_i) == 0 ) cycle
                    if( sum( abs( T_E3E3(k_s,jh_s,k_i,jh_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                  endif

                    Tens(:) = (0._db,0._db)
                    Tens_m(:) = (0._db,0._db)

                    if( icheck > 1 ) write(3,140) k_s, j_s, h_s, k_i, j_i, h_i

! Loops over spherical components of the tensors
                  Lm_s = 0
                  do L_s = 0,max(index_rank_s,1)
                    do m_s = -L_s,L_s
                      Lm_s = Lm_s + 1

                      clm_s = clm(index_rank_s,k_s,j_s,h_s,Lm_s)
                      if( abs(clm_s) < eps10 ) cycle

                      Lm_i = 0
                      do L_i = 0,max(index_rank_i,1)
                        do m_i = -L_i,L_i
                          Lm_i = Lm_i + 1

                          clm_i = clm(index_rank_i,k_i,j_i,h_i,Lm_i)
                          if( abs(clm_i) < eps10 ) cycle

! Calcul de la composante du tenseur
                          if( sum( abs( Tens_lm(Lm_s,Lm_i,:) ) ) < 1.e-15_db ) then

                            if( icheck > 1 ) write(3,150) L_s, m_s, L_i, m_i, clm_s, clm_i
                            
                            call Tens_op( Chi, Core_resolved, Final_tddft, icheck, index_rank_i, index_rank_s, ip_max, ip0, &
                                L_s, m_s, L_i, m_i, Lmax, Lmoins1, Lplus1, M_depend, ns_dipmag, ninitr, nlm1g, nlm2g, nlm_probe, &
                                nlm_p_fp, nspinp, nspino, roff_rr, Spinorbite, Taull, Ten, Ten_m, Ylm_comp )

                            Tens_lm(Lm_s,Lm_i,:) = Ten(:)
                            Tens_lm_m(Lm_s,Lm_i,:) = Ten_m(:)

                          endif

                          Tens(:) = Tens(:) + clm_s * clm_i * Tens_lm(Lm_s,Lm_i,:)
                          Tens_m(:) = Tens_m(:) + clm_s * clm_i * Tens_lm_m(Lm_s,Lm_i,:)
                        end do
                      end do
                    end do
                  end do ! fin boucle L_s

! Remplissage de la valeur calculee dans tous les elements du tenseur
! equivalent par symetrie.

                  Tensr(:) = real( Tens(:), db )
                  Tensi(:) = aimag( Tens(:) )

! M1-M1 (dipole magnetique - dipole magnetique)
                  if( index_rank_s == 0 .and. index_rank_i == 0 ) then
                    T_M1M1(k_s,k_i,:,mpirank) = Tens(:)
                    if( Green_int ) then
                      T_M1M1(k_i,k_s,:,mpirank) = Tens(:)
                      T_M1M1_m(k_s,k_i,:,mpirank) = Tens_m(:)
                      T_M1M1_m(k_i,k_s,:,mpirank) = - Tens_m(:)
                    else
                      T_M1M1(k_i,k_s,:,mpirank) = conjg( Tens(:) )
                    endif
                  endif

! M1-E1 (dipole magnetique - dipole electrique)
                  if( index_rank_s == 0 .and. index_rank_i == 1 ) then
                    T_E1M1(k_s,k_i,:,mpirank) = Tens(:)
                    if( Green_int ) T_E1M1_m(k_s,k_i,:,mpirank) = Tens_m(:)
                  endif

! E1-E1 (Dipole-dipole)
                  if( index_rank_s == 1 .and. index_rank_i == 1 ) then

                    do kk_s = 1,3
                      do kk_i = kk_s,3
                        if( abs(msymdd(kk_s,kk_i)) /= abs(msymdd(k_s,k_i)) .or. abs(msymddi(kk_s,kk_i)) &
                           /= abs(msymddi(k_s,k_i)) ) cycle
                        if( msymdd(kk_s,kk_i) /= 0 ) then
                          is = msymdd(kk_s,kk_i) / msymdd(k_s,k_i)
                        else
                          is = 0
                        endif
                        if( msymddi(kk_s,kk_i) /= 0 ) then
                          isi = msymddi(kk_s,kk_i) / msymddi(k_s,k_i)
                        else
                          isi = 0
                        endif
                        if( Green_int ) then
                          T_E1E1(kk_s,kk_i,:,mpirank) = is * Tens(:)
                          T_E1E1(kk_i,kk_s,:,mpirank) = is * Tens(:)
                          T_E1E1_m(kk_s,kk_i,:,mpirank) = isi * Tens_m(:)
                          T_E1E1_m(kk_i,kk_s,:,mpirank) = -isi * Tens_m(:)
                        else
! tenseur hermitique
                          Te(:) = cmplx( is*Tensr(:), isi*Tensi(:), db)
                          T_E1E1(kk_s,kk_i,:,mpirank) = Te(:)
                          T_E1E1(kk_i,kk_s,:,mpirank) = conjg( Te(:) )
                        endif
                      end do
                    end do

! Dipole-quadrupole
                  elseif( index_rank_s == 1 .and. index_rank_i == 2 ) then
                    do kk_s = 1,3
                      do kk_i = 1,3
                        do jj_i = kk_i,3
                          if( abs(msymdq(kk_s,kk_i,jj_i)) /= abs(msymdq(k_s,k_i,j_i)) .or. abs(msymdqi(kk_s,kk_i,jj_i)) &
                             /= abs(msymdqi(k_s,k_i,j_i)) ) cycle
                          if( msymdq(kk_s,kk_i,jj_i) /= 0 ) then
                            is = msymdq(kk_s,kk_i,jj_i) / msymdq(k_s,k_i,j_i)
                          else
                            is = 0
                          endif
                          if( msymdqi(kk_s,kk_i,jj_i) /= 0 ) then
                            isi = msymdqi(kk_s,kk_i,jj_i) / msymdqi(k_s,k_i,j_i)
                          else
                            isi = 0
                          endif
                          if( Green_int ) then
                            T_E1E2(kk_s,kk_i,jj_i,:,mpirank) =is*Tens(:)
                            T_E1E2(kk_s,jj_i,kk_i,:,mpirank) =is*Tens(:)
                            T_E1E2_m(kk_s,kk_i,jj_i,:,mpirank) = isi*Tens_m(:)
                            T_E1E2_m(kk_s,jj_i,kk_i,:,mpirank) = isi*Tens_m(:)
                          else
                            Te(:)=cmplx(is*Tensr(:),isi*Tensi(:),db)
                            T_E1E2(kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                            T_E1E2(kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                          endif
                        end do
                      end do
                    end do

! Dipole-octupole
                  elseif( index_rank_s == 1 .and. index_rank_i == 3 ) then
                    do kk_s = 1,3
                      do kk_i = 1,3
                        do jj_i = 1,3
                          do hh_i = 1,3
                            if( abs( msymdo(kk_s,kk_i,jj_i,hh_i) ) /= abs( msymdo(k_s,k_i,j_i,h_i) )  &
                              .or. abs( msymdoi(kk_s,kk_i,jj_i,hh_i) ) /= abs( msymdoi(k_s,k_i,j_i,h_i) ) ) cycle
                            if( msymdo(kk_s,kk_i,jj_i,hh_i) /= 0 ) then
                              is = msymdo(kk_s,kk_i,jj_i,hh_i) / msymdo(k_s,k_i,j_i,h_i)
                            else
                              is = 0
                            endif
                            if( msymdoi(kk_s,kk_i,jj_i,hh_i)/=0) then
                              isi = msymdoi(kk_s,kk_i,jj_i,hh_i) / msymdoi(k_s,k_i,j_i,h_i)
                            else
                              isi = 0
                            endif
                            if( Green_int ) then
                              T_E1E3(kk_s,kk_i,jj_i,hh_i,:,mpirank) = is * Tens(:)
                              T_E1E3(kk_s,kk_i,hh_i,jj_i,:,mpirank) = is * Tens(:)
                              T_E1E3(kk_s,jj_i,kk_i,hh_i,:,mpirank) = is * Tens(:)
                              T_E1E3(kk_s,hh_i,kk_i,jj_i,:,mpirank) = is * Tens(:)
                              T_E1E3(kk_s,jj_i,hh_i,kk_i,:,mpirank) = is * Tens(:)
                              T_E1E3(kk_s,hh_i,jj_i,kk_i,:,mpirank) = is * Tens(:)
                              T_E1E3_m(kk_s,kk_i,jj_i,hh_i,:,mpirank) = isi * Tens_m(:)
                              T_E1E3_m(kk_s,kk_i,hh_i,jj_i,:,mpirank) = isi * Tens_m(:)
                              T_E1E3_m(kk_s,jj_i,kk_i,hh_i,:,mpirank) = isi * Tens_m(:)
                              T_E1E3_m(kk_s,hh_i,kk_i,jj_i,:,mpirank) = isi * Tens_m(:)
                              T_E1E3_m(kk_s,jj_i,hh_i,kk_i,:,mpirank) = isi * Tens_m(:)
                              T_E1E3_m(kk_s,hh_i,jj_i,kk_i,:,mpirank) = isi * Tens_m(:)
                            else
                              Te(:) = cmplx(is*Tensr(:),isi*Tensi(:),db)
                              T_E1E3(kk_s,kk_i,jj_i,hh_i,:,mpirank) = Te(:)
                              T_E1E3(kk_s,kk_i,hh_i,jj_i,:,mpirank) = Te(:)
                              T_E1E3(kk_s,jj_i,kk_i,hh_i,:,mpirank) = Te(:)
                              T_E1E3(kk_s,hh_i,kk_i,jj_i,:,mpirank) = Te(:)
                              T_E1E3(kk_s,jj_i,hh_i,kk_i,:,mpirank) = Te(:)
                              T_E1E3(kk_s,hh_i,jj_i,kk_i,:,mpirank) = Te(:)
                            endif
                          end do
                        end do
                      end do
                    end do

! Quadrupole-quadrupole
                  elseif( index_rank_s == 2 .and. index_rank_i == 2 ) then
                    do kk_s = 1,3
                      do jj_s = 1,3
                        do kk_i = 1,3
                          do jj_i = 1,3
                            if( sum( abs( T_E2E2(kk_s,jj_s,kk_i,jj_i,:,mpirank) ) ) > 1.e-15_db ) cycle
                            if( abs( msymqq(kk_s,jj_s,kk_i,jj_i) ) /= abs( msymqq(k_s,j_s,k_i,j_i) ) &
                                .or. abs( msymqqi(kk_s,jj_s,kk_i,jj_i) ) /= abs( msymqqi(k_s,j_s,k_i,j_i) ) ) cycle
                            if( msymqq(kk_s,jj_s,kk_i,jj_i) /= 0 ) then
                              is = msymqq(kk_s,jj_s,kk_i,jj_i) / msymqq(k_s,j_s,k_i,j_i)
                            else
                              is = 0
                            endif
                            if( msymqqi(kk_s,jj_s,kk_i,jj_i) /= 0 ) then
                              isi = msymqqi(kk_s,jj_s,kk_i,jj_i) / msymqqi(k_s,j_s,k_i,j_i)
                            else
                              isi = 0
                            endif
                            if( Green_int ) then
                              Te(:) = is * Tens(:)
                              T_E2E2(kk_s,jj_s,kk_i,jj_i,:,mpirank)=Te(:)
                              T_E2E2(jj_s,kk_s,kk_i,jj_i,:,mpirank)=Te(:)
                              T_E2E2(kk_s,jj_s,jj_i,kk_i,:,mpirank)=Te(:)
                              T_E2E2(jj_s,kk_s,jj_i,kk_i,:,mpirank)=Te(:)
                              T_E2E2(kk_i,jj_i,kk_s,jj_s,:,mpirank)=Te(:)
                              T_E2E2(jj_i,kk_i,kk_s,jj_s,:,mpirank)=Te(:)
                              T_E2E2(kk_i,jj_i,jj_s,kk_s,:,mpirank)=Te(:)
                              T_E2E2(jj_i,kk_i,jj_s,kk_s,:,mpirank)=Te(:)
                              Te(:) = isi * Tens_m(:)
                              T_E2E2_m(kk_s,jj_s,kk_i,jj_i,:,mpirank) = Te(:)
                              T_E2E2_m(jj_s,kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                              T_E2E2_m(kk_s,jj_s,jj_i,kk_i,:,mpirank) = Te(:)
                              T_E2E2_m(jj_s,kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                              T_E2E2_m(kk_i,jj_i,kk_s,jj_s,:,mpirank) = - Te(:)
                              T_E2E2_m(jj_i,kk_i,kk_s,jj_s,:,mpirank) = - Te(:)
                              T_E2E2_m(kk_i,jj_i,jj_s,kk_s,:,mpirank) = - Te(:)
                              T_E2E2_m(jj_i,kk_i,jj_s,kk_s,:,mpirank) = - Te(:)
                            else
                              Te(:)= cmplx( is*Tensr(:), isi*Tensi(:),db)
                              T_E2E2(kk_s,jj_s,kk_i,jj_i,:,mpirank) = Te(:)
                              T_E2E2(jj_s,kk_s,kk_i,jj_i,:,mpirank) = Te(:)
                              T_E2E2(kk_s,jj_s,jj_i,kk_i,:,mpirank) = Te(:)
                              T_E2E2(jj_s,kk_s,jj_i,kk_i,:,mpirank) = Te(:)
                              Te(:) = Conjg( Te(:) )
                              T_E2E2(kk_i,jj_i,kk_s,jj_s,:,mpirank) = Te(:)
                              T_E2E2(jj_i,kk_i,kk_s,jj_s,:,mpirank) = Te(:)
                              T_E2E2(kk_i,jj_i,jj_s,kk_s,:,mpirank) = Te(:)
                              T_E2E2(jj_i,kk_i,jj_s,kk_s,:,mpirank) = Te(:)
                            endif
                          end do
                        end do
                      end do
                    end do
! Octupole-octupole
                  elseif( index_rank_s == 3 .and. index_rank_i == 3 ) then

                    do kk_s = 1,3
                      do jj_s = 1,3
                        do hh_s = 1,3
                          jjh_s = 3 * ( jj_s - 1 ) + hh_s
                          do kk_i = 1,3
                            do jj_i = 1,3
                              do hh_i = 1,3
                                jjh_i = 3 * ( jj_i - 1 ) + hh_i
                                i = msymoo(kk_s,jjh_s,kk_i,jjh_i)
                                j = msymoo(k_s,jh_s,k_i,jh_i)
                                if( abs(i) /= abs(j) ) cycle
                                if( j /= 0 ) then
                                  is = i / j
                                else
                                  is = 0
                                endif
                                i = msymooi(kk_s,jjh_s,kk_i,jjh_i)
                                j = msymooi(k_s,jh_s,k_i,jh_i)
                                if( abs(i) /= abs(j) ) cycle
                                if( j /= 0 ) then
                                  isi = i / j
                                else
                                  isi = 0
                                endif
                                if( Green_int ) then
                                  Te(:) = is * Tens(:)
                                  Te_m(:) = isi * Tens_m(:)
                                else
                                  Te(:)= cmplx( is * Tensr(:), isi * Tensi(:), db )
                                  Te_m(:) = Conjg( Te(:) )
                                endif

                                T_E3E3(kk_s,jjh_s,kk_i,jjh_i,:,mpirank) = Te(:)
                                if( Green_int ) T_E3E3_m(kk_s,jjh_s,kk_i,jjh_i,:,mpirank) = Te_m(:)

                              end do
                            end do
                          end do
                        end do
                      end do
                    end do

                  endif

                end do
              end do
            end do
          end do
        end do
      end do

      if( icheck < 2 ) cycle

      if( .not. E1E2 .and. index_rank_s == 1 .and. index_rank_i == 2 ) cycle
      if( index_rank_i == 3 .and. index_rank_s == 2 ) cycle
      if( index_rank_s == 0 .and. index_rank_i > 1 ) cycle
      write(3,160)
      Lm_s = 0
      do L_s = 0,max(index_rank_s,1)
        do m_s = -L_s,L_s
         Lm_s = Lm_s + 1
         Lm_i = 0
         do L_i = 0,max(index_rank_i,1)
            do m_i = -L_i,L_i
              Lm_i = Lm_i + 1
              if( sum( abs(Tens_lm(Lm_s,Lm_i,:)) ) > eps15 ) write(3,170) index_rank_s, index_rank_i, L_s, m_s, L_i, m_i, &
                    ( Tens_lm(Lm_s,Lm_i,initr), initr = 1,ninitr )
              if( .not. Green_int ) cycle
              if( sum( abs(Tens_lm_m(Lm_s,Lm_i,:)) ) > eps15 ) write(3,180) index_rank_s, index_rank_i, L_s, m_s, L_i, m_i, &
                   ( Tens_lm_m(Lm_s,Lm_i,initr), initr = 1,ninitr )
            end do
          end do
        end do
      end do

    end do
  end do

! Rotation to get the tensors in R1 basis

!  Rot_tem = matmul( Rot_cluster, transpose(Rot_atom_abs) )
  Rot_tem = matmul( Rot_cluster, Rot_atom_abs )
    
  if( abs( Rot_tem(1,1) - 1._db ) > eps10 .or. abs( Rot_tem(2,2) - 1._db ) > eps10 .or. abs( Rot_tem(3,3) - 1._db ) > eps10 ) then
  
    do initr = 1,ninitr
      if( E1E1 ) then
        mat2(:,:) = T_E1E1(:,:,initr,mpirank)
        call rot_tensor_2( mat2, rot_tem )
        T_E1E1(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E1E2 ) then
        mat3(:,:,:) = T_E1E2(:,:,:,initr,mpirank)
        call rot_tensor_3( mat3, rot_tem )
        T_E1E2(:,:,:,initr,mpirank) = mat3(:,:,:)
      endif
  
      if( E2E2 ) then
        mat4(:,:,:,:) = T_E2E2(:,:,:,:,initr,mpirank)
        call rot_tensor_4( mat4, rot_tem )
        T_E2E2(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1E3 ) then
        mat4(:,:,:,:) = T_E1E3(:,:,:,:,initr,mpirank)
        call rot_tensor_4( mat4, rot_tem )
        T_E1E3(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1M1 ) then
        mat2(:,:) = T_E1M1(:,:,initr,mpirank)
        call rot_tensor_2( mat2, rot_tem )
        T_E1M1(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( M1M1 ) then
        mat2(:,:) = T_M1M1(:,:,initr,mpirank)
        call rot_tensor_2( mat2, rot_tem )
        T_M1M1(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E3E3 ) then
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                Mat6(:,j_s,h_s,:,j_i,h_i) = T_E3E3(:,jh_s,:,jh_i,initr,mpirank)
              end do
            end do
          end do
        end do
        call rot_tensor_6( Mat6, Rot_tem )
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                T_E3E3(:,jh_s,:,jh_i,initr,mpirank) = Mat6(:,j_s,h_s,:,j_i,h_i)
              end do
            end do
          end do
        end do
      endif
  
      if( .not. Green_int ) Cycle
  
      if( E1E1 ) then
        mat2(:,:) = T_E1E1_m(:,:,initr,mpirank)
        call rot_tensor_2( mat2, rot_tem )
        T_E1E1_m(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E1E2 ) then
        mat3(:,:,:) = T_E1E2_m(:,:,:,initr,mpirank)
        call rot_tensor_3( mat3, rot_tem )
        T_E1E2_m(:,:,:,initr,mpirank) = mat3(:,:,:)
      endif
  
      if( E2E2 ) then
        mat4(:,:,:,:) = T_E2E2_m(:,:,:,:,initr,mpirank)
        call rot_tensor_4( mat4, rot_tem )
        T_E2E2_m(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1E3 ) then
        mat4(:,:,:,:) = T_E1E3_m(:,:,:,:,initr,mpirank)
        call rot_tensor_4( mat4, rot_tem )
        T_E1E3_m(:,:,:,:,initr,mpirank) = mat4(:,:,:,:)
      endif
  
      if( E1M1 ) then
        mat2(:,:) = T_E1M1_m(:,:,initr,mpirank)
        call rot_tensor_2( mat2, rot_tem )
        T_E1M1_m(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( M1M1 ) then
        mat2(:,:) = T_M1M1_m(:,:,initr,mpirank)
        call rot_tensor_2( mat2, rot_tem )
        T_M1M1_m(:,:,initr,mpirank) = mat2(:,:)
      endif
  
      if( E3E3 ) then
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                Mat6(:,j_s,h_s,:,j_i,h_i) = T_E3E3_m(:,jh_s,:,jh_i,initr,mpirank)
              end do
            end do
          end do
        end do
        call rot_tensor_6( Mat6, Rot_tem )
        jh_s = 0
        do j_s = 1,3
          do h_s = 1,3
            jh_s = jh_s + 1
            jh_i = 0
            do j_i = 1,3
              do h_i = 1,3
                jh_i = jh_i + 1
                T_E3E3_m(:,jh_s,:,jh_i,initr,mpirank) = Mat6(:,j_s,h_s,:,j_i,h_i)
              end do
            end do
          end do
        end do
      endif
  
    end do

  endif
    
  deallocate( roff_rr )
  
  return
  100 format(/' ---- Tenseur_car Radial ---------',100('-'))
  110 format(/' ---- Tenseur_car Tens_op --------',100('-'))
  130 format(/' -- index_rank_s =',i2,', index_rank_i =',i2,' --')
  140 format(/' k_s, j_s, h_s =',3i3,',  k_i, j_i, h_i =',3i3)
  150 format(/' L_s, m_s =',2i3,',  L_i, m_i =',2i3,', clm_s, clm_i =',1p, 2e13.5)
  160 format(/' Tensor by harmonics (basis R4) :'/, ' ir jr  L_s m_s  L_i m_i      Tens(Ylm,i=1,ninitr)')
  170 format(2i3,i4,i3,2x,i4,i3,1p,20e15.7)
  180 format('i',i2,i3,i4,i3,2x,i4,i3,1p,20e15.7)

end

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

subroutine Tens_op(Chi, Core_resolved, Final_tddft, icheck, index_rank_i, index_rank_s, ip_max, ip0, &
                                L_s, m_s, L_i, m_i, Lmax, Lmoins1, Lplus1, M_depend, ns_dipmag, ninitr, nlm1g, nlm2g, nlm_probe, &
                                nlm_p_fp, nspinp, nspino, roff_rr, Spinorbite, Taull, Ten, Ten_m, Ylm_comp )

  use declarations
  implicit none

  integer, intent(in):: icheck, index_rank_i, index_rank_s, ip_max, ip0, L_s, m_s, &
    L_i, m_i, nlm1g, nlm2g, nlm_probe, nlm_p_fp, ns_dipmag, nspinp, nspino

  integer:: i, io1, io2, is1, is12, is2, is21, iso, iso_f1, iso_f2, iso_g1, iso_g2, isp, ispm, ispm_f1, ispm_f2, &
    ispm_g1, ispm_g2, isp_f1, isp_f2, isp_g1, isp_g2, l, L_f1, L_f2, L_g1, L_g2, lm, Lm_f1, Lm_f2, lm_g1, lm_g2, &
    Lmax, lmp_f1, lmp_f2, lmp_g1, lmp_g2, Lm_i, lms_f1, lms_f2, lms_g1, lmp, lmp0, lms_g2, lmss, lmss_f1, lmss_f2, lmss_g1, &
    lmss_g2, Lmv, Lmv_f1, Lmv_f2, Lmv_g1, Lmv_g2, lp, m, m_f1, m_f2, m_g1, m_g2, mp, mv, ninitr, nlmss

  integer, dimension(:), allocatable:: iso_val, isp_val, ispm_val, l_val, lm_val, lmp_val, lms_val, Lmv_val, m_val
 
  complex(kind=db):: cf, Gaunt_s, Gaunt_xrx, Gaunt_i, PseudoGaunt_M1, Tau_rad
  
  complex(kind=db), dimension(ninitr):: Ten, Ten_m
  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp,2):: Taull
  complex(kind=db), dimension(nlm_probe*nspino,nlm_probe*nspino,nlm_probe*nspino,nlm_probe*nspino,2,2,ns_dipmag):: Chi

  logical:: Core_resolved, Final_tddft, Lmoins1, Lplus1, M_depend, Spinorbite, Ylm_comp

  real(kind=db):: Ro
  real(kind=db), dimension(nlm1g,nlm1g,nlm2g,nlm2g,nspinp**2,nspino**2,ip0:ip_max):: roff_rr

  Ten(:) = (0._db,0._db)
  Ten_m(:) = (0._db,0._db)

! Dimension calculation -----------------------------------

  nlmss = 0
  do isp = 1,2
    do l = 0,Lmax
      do m = -l,l
! Loop over harmonics, non-spherical potential
        do lp = 0,Lmax
          if( nlm_p_fp == 1 .and. lp /= l ) cycle
          do mp = -lp,lp
            if( nlm_p_fp == 1 .and. mp /= m ) cycle
! Loop over solutions
            do iso = 1,nspino
              if( Spinorbite ) then
                mv = mp - isp + iso
                if( mv > lp .or. mv < -lp ) cycle
              endif
              nlmss = nlmss + 1
            end do
          end do
        end do
      end do
    end do
  end do

  allocate( iso_val(nlmss) )
  allocate( isp_val(nlmss) )
  allocate( ispm_val(nlmss) )
  allocate( l_val(nlmss) )
  allocate( lm_val(nlmss) )
  allocate( lmp_val(nlmss) )
  allocate( lms_val(nlmss) )
  allocate( lmv_val(nlmss) )
  allocate( m_val(nlmss) )

! Building of quantum number table -------------------------------

  lmss = 0
! Loop over harmonics
  do l = 0,Lmax

    do m = -l,l
      if( M_depend ) then
        lm = l**2 + l + 1 + m
      else
        lm = l + 1
      endif

! Loop over harmonics, non-spherical potential
      lmp = 0
      do lp = 0,Lmax
        if( nlm_p_fp == 1 .and. lp /= l ) cycle
        lmp0 = nspino * (lp**2 + lp + 1)

        do mp = -lp,lp
          if( nlm_p_fp == 1 .and. mp /= m ) cycle
          lmp = lmp + 1

          do isp = 1,2
          ispm = min( isp, nspinp )

! Loop over solutions
            do iso = 1,nspino

              if( Spinorbite ) then
                mv = mp - isp + iso
                if( mv > lp .or. mv < -lp ) cycle
              else
                mv = mp
              endif
              Lmv = Lp**2 + Lp + 1 + mv   
              
              Lm_i = lmp0 + ( mv - 1 ) * nspino + iso

              lmss = lmss + 1

              iso_val(lmss) = iso
              isp_val(lmss) = isp 
              ispm_val(lmss) = ispm
              l_val(lmss) = l
              lm_val(lmss) = lm
              lmp_val(lmss) = lmp 
              lms_val(lmss) = Lm_i 
              lmv_val(lmss) = lmv 
              m_val(lmss) = m
              
            end do
          end do
        end do
      end do
    end do
  end do

  if( icheck > 1 ) then
    if( nlm_p_fp == 1 .and. Spinorbite .and. Final_tddft ) then
      write(3,110)
    elseif( nlm_p_fp == 1 .and. Final_tddft ) then
      write(3,112)
    elseif( nlm_p_fp == 1 .and. Spinorbite ) then
      write(3,115)
    elseif( nlm_p_fp == 1 ) then
      write(3,120)
    else
      write(3,130) 
    endif
  endif

! Calculation of the tensor components -------------------------------------------

! Loop over occupied states (outgoing)
  do lmss_g1 = 1,nlmss

    iso_g1 = iso_val(lmss_g1) 
    isp_g1 = isp_val(lmss_g1) 
    ispm_g1 = ispm_val(lmss_g1)
    L_g1 = l_val(lmss_g1)
    lm_g1 = lm_val(lmss_g1)
    lmp_g1 = lmp_val(lmss_g1) 
    lms_g1 = lms_val(lmss_g1)
    Lmv_g1 = Lmv_val(lmss_g1)
    m_g1 = m_val(lmss_g1) 

! Loop over non-occupied states ( outgoing )
    do lmss_f1 = 1,nlmss

      iso_f1 = iso_val(lmss_f1) 
      isp_f1 = isp_val(lmss_f1) 
      ispm_f1 = ispm_val(lmss_f1)
      L_f1 = l_val(lmss_f1)
      Lm_f1 = lm_val(lmss_f1)
      lmp_f1 = lmp_val(lmss_f1) 
      lms_f1 = lms_val(lmss_f1) 
      Lmv_f1 = Lmv_val(lmss_f1)
      m_f1 = m_val(lmss_f1) 

! spin-flip is possible only for magnetic dipole
      if( isp_f1 /= isp_g1 .and. index_rank_s /= 0 ) cycle

      if( ( Lplus1 .and. L_f1 < L_g1 + 1 ) .or. ( Lmoins1 .and. L_f1 > L_g1 )  ) cycle
      if( L_f1 > L_g1 + index_rank_s .or. L_f1 < L_g1 - index_rank_s .or. mod(L_f1,2) /= mod(L_g1+index_rank_s,2) ) cycle

! < g1 m_s f1 >< f2 m_i g2 >
      if( index_rank_s == 0 ) then
!        Gaunt_s = PseudoGaunt_M1(L_f1,m_f1,isp_f1,m_s,L_g1,m_g1,isp_g1,Ylm_comp,Ylm_comp)
        Gaunt_s = conjg( PseudoGaunt_M1(L_g1,m_g1,isp_g1,m_s,L_f1,m_f1,isp_f1,Ylm_comp,Ylm_comp) )
      else
        Gaunt_s = Gaunt_xrx(L_f1,m_f1,L_s,m_s,L_g1,m_g1,Ylm_comp)
      endif
      if( abs(Gaunt_s) < eps10 ) cycle

! Loop over non-occupied states ( incoming )
      do lmss_f2 = 1,nlmss

        iso_f2 = iso_val(lmss_f2) 
        isp_f2 = isp_val(lmss_f2) 
        ispm_f2 = ispm_val(lmss_f2)
        L_f2 = l_val(lmss_f2)
        Lm_f2 = lm_val(lmss_f2)
        lmp_f2 = lmp_val(lmss_f2) 
        lms_f2 = lms_val(lmss_f2) 
        Lmv_f2 = Lmv_val(lmss_f2)
        m_f2 = m_val(lmss_f2) 
    
! Loop over occupied states ( incoming )
        do lmss_g2 = 1,nlmss

          iso_g2 = iso_val(lmss_g2) 
          isp_g2 = isp_val(lmss_g2) 
          ispm_g2 = ispm_val(lmss_g2)
          L_g2 = l_val(lmss_g2)
          lm_g2 = lm_val(lmss_g2)
          lmp_g2 = lmp_val(lmss_g2) 
          lms_g2 = lms_val(lmss_g2) 
          Lmv_g2 = Lmv_val(lmss_g2)
          m_g2 = m_val(lmss_g2) 

          if( isp_f2 /= isp_g2 .and. index_rank_i /= 0 ) cycle
          
!          if( .not. ( Final_tddft .or. Spinorbite ) .and. ( isp_g2 /= isp_g1 .or. isp_f2 /= isp_f1 ) ) cycle
!          if( ( isp_g2 /= isp_g1 .or. isp_f2 /= isp_f1 ) ) cycle
          
          if( ( Lplus1 .and. L_f2 < L_g2 + 1 ) .or. ( Lmoins1 .and. L_f2 > L_g2 ) ) cycle
          if( L_f2 > L_g2 + index_rank_i .or. L_f2 < L_g2 - index_rank_i .or. mod(L_f2,2) /= mod(L_g2+index_rank_i,2) ) cycle

! < g1 m_s f1 >< f2 m_i g2 >
          if( index_rank_i == 0 ) then
            Gaunt_i = PseudoGaunt_M1(L_f2,m_f2,isp_f2,m_i,L_g2,m_g2,isp_g2,Ylm_comp,Ylm_comp)
          else
            Gaunt_i = Gaunt_xrx(L_f2,m_f2,L_i,m_i,L_g2,m_g2,Ylm_comp)
          endif
          if( abs(Gaunt_i) < eps10 ) cycle

          is1 = ispm_f1 + (ispm_g1 - 1) * nspinp
          is2 = ispm_f2 + (ispm_g2 - 1) * nspinp
          io1 = iso_f1 + (iso_g1 - 1) * nspino
          io2 = iso_f2 + (iso_g2 - 1) * nspino

          Ro = roff_rr(lm_g1,Lm_f1,lmp_g1,lmp_f1,is1,io1,index_rank_s) * roff_rr(lm_g2,Lm_f2,lmp_g2,lmp_f2,is2,io2,index_rank_i)

          if( Final_tddft ) then
            is12 = 2 * isp_g1 - 2 + isp_g2
            is21 = 2 * isp_g2 - 2 + isp_g1
     ! + conjg is because this Taull must be multiplied by - i to be classical tau like.
            Tau_rad = - 0.5 * Ro * ( Chi(lms_f1,lms_f2,lms_g1,lms_g2,isp_f1,isp_f2,is12) &
                               + conjg( Chi(lms_f2,lms_f1,lms_g2,lms_g1,isp_f2,isp_f1,is21) ) )
          else                       
            if( Spinorbite ) then                       
              Tau_rad = - Ro * Taull(lmv_f1,iso_f1,lmv_f2,iso_f2,2) * Taull(lmv_g2,iso_g2,lmv_g1,iso_g1,1)
            else
              Tau_rad = - Ro * Taull(lmv_f1,ispm_f1,lmv_f2,ispm_f2,2) * Taull(lmv_g2,ispm_g2,lmv_g1,ispm_g1,1)
            endif
          endif

! Evenso Tau appear 2 times, one devide only 1 time by pi.
! Normalisation by pi must then do not to be done in coabs for dafs amplitude calculation.
          cf = Tau_rad * Conjg( Gaunt_s ) * Gaunt_i / pi
            
          Ten(1) = Ten(1) + cf
          if( Core_resolved ) then
            do i = 2,ninitr
              if( L_g1 == i-2 .and. L_g2 == i-2 ) Ten(i) = Ten(i) + cf
            end do
          endif

          if( icheck < 2  .or. abs( cf ) < 1.e-15_db ) cycle
          if( nlm_p_fp == 1 .and. Spinorbite .and. Final_tddft ) then
            write(3,150) L_g1, m_g1, isp_g1, iso_g1, L_f1, m_f1, isp_f1, iso_f1, &
                L_f2, m_f2, isp_f2, iso_f2, L_g2, m_g2, isp_g2, iso_g2, &
                Ten(1), cf, Chi(lms_f2,lms_f1,lms_g2,lms_g1,isp_f1,isp_f2,is12), Gaunt_s, Gaunt_i, & 
                roff_rr(lm_g1,Lm_f1,lmp_g1,lmp_f1,is1,io1,index_rank_s), &
                roff_rr(lm_g2,Lm_f2,lmp_g2,lmp_f2,is2,io2,index_rank_i)
          elseif( nlm_p_fp == 1 .and. Final_tddft ) then
            write(3,155) L_g1, m_g1, isp_g1, L_f1, m_f1, isp_f1, &
                L_f2, m_f2, isp_f2, L_g2, m_g2, isp_g2, &
                Ten(1), cf, Chi(lms_f2,lms_f1,lms_g2,lms_g1,isp_f1,isp_f2,is12), Gaunt_s, Gaunt_i, & 
                roff_rr(lm_g1,Lm_f1,lmp_g1,lmp_f1,is1,io1,index_rank_s), &
                roff_rr(lm_g2,Lm_f2,lmp_g2,lmp_f2,is2,io2,index_rank_i)
          elseif( nlm_p_fp == 1 .and. Spinorbite ) then
            write(3,150) L_g1, m_g1, isp_g1, iso_g1, L_f1, m_f1, isp_f1, iso_f1, &
                L_f2, m_f2, isp_f2, iso_f2, L_g2, m_g2, isp_g2, iso_g2, &
                Ten(1), cf, Taull(lmv_g2,iso_g2,lmv_g1,iso_g1,1), &
                Taull(lmv_f1,iso_f1,lmv_f2,iso_f2,2), Gaunt_s, Gaunt_i, & 
                roff_rr(lm_g1,Lm_f1,lmp_g1,lmp_f1,is1,io1,index_rank_s), &
                roff_rr(lm_g2,Lm_f2,lmp_g2,lmp_f2,is2,io2,index_rank_i)
          elseif( nlm_p_fp == 1 ) then
            write(3,155) L_g1, m_g1, isp_g1, L_f1, m_f1, isp_f1, &
                L_f2, m_f2, isp_f2, L_g2, m_g2, isp_g2, &
                Ten(1), cf, Taull(lmv_g2,ispm_g2,lmv_g1,ispm_g1,1), &
                Taull(lmv_f1,ispm_f1,lmv_f2,ispm_f2,2), Gaunt_s, Gaunt_i, & 
                roff_rr(lm_g1,Lm_f1,lmp_g1,lmp_f1,is1,io1,index_rank_s), &
                roff_rr(lm_g2,Lm_f2,lmp_g2,lmp_f2,is2,io2,index_rank_i)
          else
            write(3,160) L_g1, m_g1, lmp_g1, isp_g1, iso_g1, L_f1, m_f1, lmp_f1, isp_f1, iso_f1,&
                L_f2, m_f2, lmp_f2, isp_f2, iso_f2,  L_g2, m_g2, lmp_g2, isp_g2, iso_g2, &
                Ten(1), cf, Taull(lmv_f1,ispm_f1,lmv_f2,ispm_f2,2), Gaunt_s, Gaunt_i, & 
                roff_rr(lm_g1,Lm_f1,lmp_g1,lmp_f1,is1,io1,index_rank_s), &
                roff_rr(lm_g2,Lm_f2,lmp_g2,lmp_f2,is2,io2,index_rank_i)
          endif

        end do    ! end of loop over non-occupied states, incoming
      end do    ! end of loop over non-occupied states, outgoing
    end do    ! end of loop over occupied states, incoming
  end do   ! end of loop over occupied states, outgoing

  deallocate( iso_val, isp_val, ispm_val, l_val, lm_val, lmp_val, lms_val, Lmv_val, m_val )

  return
  110 format(/' lg1 mg1 sg1 og1 lf1 mf1 sf1 of1 lf2 mf2 sf2 of2 lg2 mg2 sg2 og2',11x,'Tens',21x,'d_Tens', &
              17x,'Tau_f1f2g1g2',14x,'Gaunte_g1f1',15x,'Gaunts_g2f2',10x,'rof_f1g1',5x,'rof_f2g2')
  112 format(/' lg1 mg1 sg1 lf1 mf1 sf1 lf2 mf2 sf2 lg2 mg2 sg2',11x,'Tens',21x,'d_Tens', &
              17x,'Tau_f1f2g1g2',14x,'Gaunte_g1f1',15x,'Gaunts_g2f2',10x,'rof_f1g1',5x,'rof_f2g2')
  115 format(/' lg1 mg1 sg1 og1 lf1 mf1 sf1 of1 lf2 mf2 sf2 of2 lg2 mg2 sg2 og22',11x,'Tens',21x,'d_Tens', &
              19x,'Tau_g1g2',18x,'Tau_f1f2',16x,'Gaunte_g1f1',15x,'Gaunts_g2f2',10x,'rof_f1g1',5x,'rof_f2g2')
  120 format(/' lg1 mg1 sg1 lf1 mf1 sf1 lf2 mf2 sf2 lg2 mg2 sg2',11x,'Tens',21x,'d_Tens', &
              19x,'Tau_g1g2',18x,'Tau_f1f2',16x,'Gaunte_g1f1',15x,'Gaunts_g2f2',10x,'rof_f1g1',5x,'rof_f2g2')
  130 format(/' lg1 mg1 pg1 sg1 og1 lf1 mf1 pf1 sf1 of1 lf2 mf2 pf2 sf2 of2 lg2 mg2 pg2 sg2 og2',11x,'Tens',21x,'d_Tens', &
              17x,'Tau_f1f2g1g2',14x,'Gaunte_g1f1',15x,'Gaunts_g2f2',10x,'rof_f1g1',3x,'rof_f2g2')
  150 format(i3,15i4,1p,14e13.5)
  155 format(i3,11i4,1p,14e13.5)
  160 format(i3,19i4,1p,14e13.5)
end

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

! Gaunt coefficient calculation with complex harmonics = Int( Y(l1,m1)*Y(l2,m2)Y(l3,m3)dOmega )
! Formula from M. E. Rose p. 62

function Gauntcp(l1,m1,l2,m2,l3,m3)

  use declarations
  
  implicit none
  
  integer:: l1, l2, l3, m1, m2, m3
  
  real(kind=db):: cgc, fac, Gauntcp

  if( m1 /= m2 + m3 .or. mod(l1+l2+l3,2) == 1 .or. l1 < abs(l2-l3) .or. l1 > abs(l2+l3) ) then

    Gauntcp = 0._db

  else

    fac = sqrt( (2*l3 + 1._db) * (2*l2 + 1._db) / ( (2*l1 + 1._db) * quatre_pi ) )
    Gauntcp = fac * cgc(l3,l2,l1,m3,m2) * cgc(l3,l2,l1,0,0)

  endif

  return
end

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

! Clebsch-gordan coefficient  eq. 3.18, Rose

function cgc(L1,L2,L3,m1,m2)

  use declarations
  implicit none

  integer:: i, i1, id, inum, j, j0, jm, L1, L2, L3, m1, m2, m3, n, n1, nfac, nff, nmin, nmax, ntemp, nu
  
  integer, dimension(5):: nd, num
  
  real(kind=db):: cgc, cgcp, ff

  cgc = 0._db

  m3 = m1 + m2

! Arguments of factorials
  num( 1 ) = L3 + L1 - L2
  num( 2 ) = L3 - L1 + L2
  num( 3 ) = L1 + L2 - L3
  num( 4 ) = L3 + m3
  num( 5 ) = L3 - m3
  nd( 1 ) = L1 + L2 + L3 + 1
  nd( 2 ) = L1 - m1
  nd( 3 ) = L1 + m1
  nd( 4 ) = L2 - m2
  nd( 5 ) = L2 + m2

! Check triangle and projection conditions
  do i = 1, 5
    if( num(i) < 0 .or. nd(i) < 0 ) return
  end do

  nff = 0
  nmin = max( 0, L2 + m3 - L1 )
  ff = 1.0_db

! Two sets of factorial products
  n = 5
  do nfac = 1, 2
    n1 = n - 1
! Arrange arguments in descending order
    do i = 1, n1
      inum = i
      id = i
      i1 = i + 1
      do j = i1, n
        if( num( j ) > num( inum ) ) inum = j
        if( nd( j ) > nd( id ) ) id = j
      end do
      ntemp = num( i )
      num( i ) = num( inum )
      num( inum ) = ntemp
      ntemp = nd( i )
      nd( i ) = nd( id )
      nd( id ) = ntemp
    end do
! Compute factorial ratios
    do i = 1, n
      if( num( i ) < nd( i ) ) then
        jm = nd( i )
        if( jm == 1 ) cycle
        j0 = num( i ) + 1
        if( num( i ) == 0 ) j0 = 2
        do j = j0, jm
          if( abs ( ff ) < 1.0e-20_db ) then
            ff = ff * 1.e20_db
            nff = nff - 2
          endif
          ff = ff / j
        end do
      elseif( num( i ) > nd( i ) ) then
        jm = num( i )
        if( jm == 1 ) cycle
        j0 = nd( i ) + 1
        if( nd( i ) == 0 ) j0 = 2
        do j = j0, jm
          if( abs ( ff ) > 1.0e20_db ) then
            ff = ff / 1.0e20_db
            nff = nff + 2
          endif
          ff = j * ff
        end do
      endif
    end do

    if ( nfac == 2 ) exit

    nff = nff / 2
    ff = sqrt( ( 2 * L3 + 1 ) * ff )

! Second set of factorial arguments
    num( 1 ) = L2 + L3 + m1 - nmin
    num( 2 ) = L1 - m1 + nmin
    num( 3 ) = 0
    nd( 1 ) = nmin
    if( nmin == 0 ) nd( 1 ) = L1 - L2 - m3
    nd( 2 ) = L3 - L1 + L2 - nmin
    nd( 3 ) = L3 + m3 - nmin
    n = 3
  end do

  if( mod( nmin + L2 + m2, 2 ) /= 0 ) ff = - ff
  ff = ff * 1.0e10_db**nff
  cgcp = ff
  nmax = min( L3 - L1 + L2, L3 + m3 )

  do nu = nmin+1, nmax
    ff = - ff * ( (L1-m1+nu) * (L3-L1+L2-nu+1) * (L3+m3-nu+1) ) / real( nu * (nu+L1-L2-m3) * (L2+L3+m1-nu+1), db )
    cgcp = cgcp + ff
  end do

  cgc = cgcp

  return
end

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

! This subroutine calculates the Gaunt coefficients for real spherical harmonics

function Gaunt_r( l1, m1, l2, m2, l3, m3 )

  use declarations
  implicit none

  integer:: itr, ityfct1, ityfct2, ityfct3, itygam, ityp, itypi, l1, l2, l3, ltsti, ltst, m, m1, m2, m3, mgam, mi, nexp, ns, nsf
 
  real(kind=db):: ab, cgc, pre, Gaunt_r, gf, sqr2i
  
  if( mod( l3+l1+l2, 2 ) /= 0 .or. abs( l1 - l3 ) > l2 .or. ( l1 + l3 ) < l2 ) then
    Gaunt_r = 0._db
    return
  endif

  mi = abs( m1 )
  mgam = abs( m2 )
  m = abs( m3 )

  if( mi+mgam /= m .and. mi+m /= mgam .and. m+mgam /= mi ) then
    Gaunt_r = 0._db
    return
  endif

  sqr2i = 1 / sqrt( 2._db )

  ltsti = l1
  itr = l2
  ltst = l3

  pre = sqrt( ( 2*ltst+1 ) * ( 2*itr+1 ) / ( ( 2*ltsti +1 ) * quatre_pi ) ) * cgc(ltst,itr,ltsti,0,0)

  if( m1 == 0 ) then
    itypi = 0
  else
    itypi = mi / m1
  endif
  if( m2 == 0 ) then
    itygam = 0
  else
    itygam = mgam / m2
  endif
  if( m3 == 0 ) then
    ityp = 0
  else
    ityp = m / m3
  endif

  nexp = abs(itypi) + abs(itygam) + abs(ityp)
  ns = itypi + ityp + itygam
  if( ns == 0 ) ns = 1
  nsf = isign(1,ns)
  gf = nsf * (-1)**mi * sqr2i**nexp * pre

  if( m+mgam+mi == 0 ) then
    ab = cgc(ltst,itr,ltsti,0,0)
  else
    ityfct1 = itypi + itygam * ityp
    ityfct2 = itygam + itypi * ityp
    ityfct3 = ityp + itypi * itygam
    ab = 0._db
    if( m+mgam == mi .and. ityfct1 /= 0 ) ab = (-1)**mi * cgc(ltst,itr,ltsti,m,mgam) * ityfct1
    if( m+mi == mgam .and. ityfct2 /= 0 ) ab = ab + (-1)**mgam * cgc(ltst,itr,ltsti,-m,mgam) * ityfct2
    if( mgam+mi == m .and. ityfct3 /= 0 ) ab = ab + (-1)**m * cgc(ltst,itr,ltsti,m,-mgam) * ityfct3
  endif

  Gaunt_r = gf * ab

  return
end

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

! Calculation lof Gaunt coefficient with Y(l,m) complex or real, Y(lo,mo) real and Y(L_g,m_g) complex. 

function Gaunt_xrc(l,m,lo,mo,L_g,m_g,Ylm_comp)

  use declarations
  implicit none
  
  integer:: l, L_g, lo, m, m_g, mo
  
  complex(kind=db):: Gaunt_xrc

  logical:: Ylm_comp

  real(kind=db):: Gaunt_r, Gauntcp, gi, gr
  
  if( Ylm_comp ) then

    if( mo == 0 ) then
      gr = gauntcp(l,m,lo,mo,L_g,m_g)
      gi = 0._db
    elseif( mo > 0 ) then
      gr = ( (-1)**mo * gauntcp(l,m,lo,mo,L_g,m_g) + gauntcp(l,m,lo,-mo,L_g,m_g) ) / sqrt( 2._db )
      gi = 0._db
    else
      gi = ( gauntcp(l,m,lo,mo,L_g,m_g) - (-1)**mo * gauntcp(l,m,lo,-mo,L_g,m_g) ) / sqrt( 2._db )
      gr = 0._db
    endif

  else

    if( m_g == 0 ) then
      gr = Gaunt_r(l,m,lo,mo,L_g,m_g)
      gi = 0._db
    elseif( m_g > 0 ) then
      gr = (-1)**m_g * Gaunt_r(l,m,lo,mo,L_g,m_g) / sqrt( 2._db )
      gi = (-1)**m_g * Gaunt_r(l,m,lo,mo,L_g,-m_g) / sqrt( 2._db )
    else
      gr = Gaunt_r(l,m,lo,mo,L_g,-m_g) / sqrt( 2._db )
      gi = - Gaunt_r(l,m,lo,mo,L_g,m_g) / sqrt( 2._db )
    endif

  endif

  Gaunt_xrc = cmplx( gr, gi, db )

  return
end

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

! Calcule le coefficient de Gaunt avec Y(l1,m1) complexe ou reel, Y(l2,m2) reel et Y(l3,m3) complexe ou reel

function Gaunt_xrx(l1,m1,l2,m2,l3,m3,Ylm_comp)

  use declarations
  implicit none
  
  integer:: l1, l2, l3, m1, m2, m3
  
  complex(kind=db):: Gaunt_xrx

  logical Ylm_comp

  real(kind=db):: Gaunt_r, Gauntcp, Gtot

  if( Ylm_comp ) then

    if( m2 == 0 ) then
      Gaunt_xrx = cmplx( gauntcp(l1,m1,l2,m2,l3,m3), 0._db, db )
    elseif( m2 > 0 ) then
      Gtot = ( Gauntcp(l1,m1,l2,m2,l3,m3) + Gauntcp(l1,m1,l2,-m2,l3,m3) ) / sqrt( 2._db ) 
      Gaunt_xrx = cmplx( Gtot, 0._db, db )
    else
      Gtot = ( (-1)**m2 * Gauntcp(l1,m1,l2,m2,l3,m3) - (-1)**m2 * Gauntcp(l1,m1,l2,-m2,l3,m3) ) / sqrt( 2._db )
      Gaunt_xrx = cmplx( 0._db, Gtot, db )
    endif

  else

    Gaunt_xrx = cmplx( Gaunt_r(l1,m1,l2,m2,l3,m3), 0._db, db )

  endif

  return
end

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

! Calcule le coefficient de Gaunt modifie: Int ( Y_L1* Y_L2* Y_L3 dOmega )
! avec Y(l1,m1) complexe ou reel, Y(l2,m2) et Y(l3,m3) complexe
! Pour Y(l2,m2), on prend Y(l2,m2)* = (-1)^m2 Y(l2,-m2)  

function Gaunt_nrixs(l1,m1,l2,m2,l3,m3,Ylm_comp)

  use declarations
  implicit none
  
  integer:: l1, l2, l3, m1, m2, m3
  
  complex(kind=db):: Gaunt_nrixs

  logical Ylm_comp

  real(kind=db):: Gauntcp, Gtot

  if( Ylm_comp .or. m1 == 0 ) then
    Gaunt_nrixs = (-1)**m2 * cmplx( gauntcp(l1,m1,l2,-m2,l3,m3), 0._db, db )
  elseif( m1 > 0 ) then
    Gtot = (-1)**m2 * ( (-1)**m1 * Gauntcp(l1,m1,l2,-m2,l3,m3) + Gauntcp(l1,-m1,l2,-m2,l3,m3) ) / sqrt( 2._db ) 
    Gaunt_nrixs = cmplx( Gtot, 0._db, db )
  else
    Gtot = (-1)**m2 * ( Gauntcp(l1,m1,l2,-m2,l3,m3) - (-1)**m1 * Gauntcp(l1,-m1,l2,-m2,l3,m3) ) / sqrt( 2._db )
    Gaunt_nrixs = cmplx( 0._db, Gtot, db )
  endif

  return
end

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

! Calculation of the pseudo Gaunt coefficient for the magnetic transition M1
! G = 1/hbar < (l,m,isp_f) | L_o + 2 S_o | (l_g,m_g,isp_g) >

function PseudoGaunt_M1(L,m,isp_f,mo,L_g,m_g,isp_g,Ylm_comp,Ylm_comp_g)

  use declarations
  implicit none

  integer:: i, isp_f, isp_g, j, L, L_g, m, m_g, mo, m_cf, m_cg

  complex(kind=db):: PseudoGaunt_M1

  logical:: Ylm_comp, Ylm_comp_g

  real(kind=db):: f

  PseudoGaunt_M1 = ( 0._db, 0._db )

  if( L /= L_g ) return

! Loop for real harmonics, right side
  do j = 1,2

    if( j == 1 ) then
      m_cg = m_g 
    else
      m_cg = - m_g
    endif

! Loop for real harmonics, left side
    do i = 1,2

      if( i == 1 ) then
        m_cf = m 
      else
        m_cf = - m
      endif

      f = 0._db

      select case(mo)

        case(-1)     ! Ly + 2*Sy
          if( isp_g == isp_f ) then
            if( m_cf == m_cg + 1 ) then
              f = - 0.5_db * sqrt( L_g * (L_g + 1._db) - m_cg * (m_cg + 1) )
            elseif( m_cf == m_cg - 1 ) then
              f = 0.5_db * sqrt( L_g * (L_g + 1._db) - m_cg * (m_cg - 1) )
            endif
          elseif( m_cf == m_cg ) then
            if( isp_g == 1 ) then
              f = 1._db
            else
              f = - 1._db
            endif
          endif

        case(0)     ! Lz + 2*Sz
          if( isp_g == isp_f .and. m_cf == m_cg ) then
            if( isp_g == 1 ) then
              f = m_cg + 1._db
            else
              f = m_cg - 1._db
            endif
          endif

        case(1)     ! Lx + 2*Sx
          if( isp_g == isp_f ) then
            if( m_cf == m_cg + 1 ) then
              f = 0.5_db * sqrt( L_g*(L_g + 1._db) - m_cg * (m_cg + 1) )
            elseif( m_cf == m_cg - 1 ) then
              f = 0.5_db * sqrt( L_g*(L_g + 1._db) - m_cg * (m_cg - 1) )
            endif
          elseif( m_cf == m_cg ) then
            f = 1._db
          endif

      end select

! For the conversion from real to complex
      if( .not. ( Ylm_comp .or. m == 0 ) ) then
        f = f / sqrt( 2._db )
        if(  m > 0 ) then
          if( i == 1 .and. mod(abs(m),2) == 1 ) f = - f
        else
          if( i == 2 .and. mod(abs(m),2) == 0 ) f = - f
        endif
      endif

      if( .not. ( Ylm_comp_g .or. m_g == 0 ) ) then
        f = f / sqrt( 2._db )
        if(  m_g > 0 ) then
          if( j == 1 .and. mod(abs(m_g),2) == 1 ) f = - f
        else
          if( j == 2 .and. mod(abs(m_g),2) == 0 ) f = - f
        endif
      endif

      if( mo == -1 ) then
        PseudoGaunt_M1 = PseudoGaunt_M1 + cmplx( 0._db, f, db )
      else
        PseudoGaunt_M1 = PseudoGaunt_M1 + cmplx( f, 0._db, db )
      endif

      if( Ylm_comp .or. m == 0 ) exit

    end do

    if( Ylm_comp_g .or. m_g == 0 ) exit
    
  end do

! For the conversion back from complex to real  
! Minus sign is because it is in <bra, so complex conjugate
  if( .not. Ylm_comp ) then
    if( m < 0 ) PseudoGaunt_M1 = - img * PseudoGaunt_M1
!    if( mod(abs(m),2) == 1 ) PseudoGaunt_M1 = - PseudoGaunt_M1
  endif

  if( .not. Ylm_comp_g ) then
    if( m_g < 0 ) PseudoGaunt_M1 = img * PseudoGaunt_M1
!    if( mod(abs(m_g),2) == 1 ) PseudoGaunt_M1 = - PseudoGaunt_M1
  endif

  PseudoGaunt_M1 = 0.5_db * alfa_sf * PseudoGaunt_M1

  return
end

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

! Calculation of the pseudo Gaunt coefficient for the transition transition M2
! G = i * alfa * /(hbar k) < (L,m,isp_f) | ( 2/3L_o + 2 S_o ) k.r | (L_g,m_g,isp_g) >
! Note that the "i" in front of k is kept here

function PseudoGaunt_M2(L,m,isp_f,m_o,m_k,L_g,m_g,isp_g,Ylm_comp,Ylm_comp_g)

  use declarations
  implicit none

  integer, parameter:: Two_third = 2._db / 3._db
  
  integer:: i, isp_f, isp_g, j, L, L_g, L_o, m, m_g, m_o, m_cf, m_cg, m_g_m_1, m_g_p_1, m_k

  complex(kind=db):: Cfac, Gaunt_xrc, PseudoGaunt_M2

  logical:: Ylm_comp, Ylm_comp_g

  real(kind=db):: p_plus, p_minus

  PseudoGaunt_M2 = ( 0._db, 0._db )

  L_o = 1
  if( L /= L_g + L_o .and. L /= L_g - L_o ) return

! Loop for real harmonics, right side
  do j = 1,2

    if( j == 1 ) then
      m_cg = m_g 
    else
      m_cg = - m_g
    endif

! Loop for real harmonics, left side
    do i = 1,2

      if( i == 1 ) then
        m_cf = m 
      else
        m_cf = - m
      endif
      
      if( m_o /= 0 ) then
        p_plus = 0.5_db * sqrt( L_g * (L_g + 1._db) - m_cg * (m_cg + 1) )
        p_minus = 0.5_db * sqrt( L_g * (L_g + 1._db) - m_cg * (m_cg - 1) )  
      endif

      Cfac = ( 0._db, 0._db )
            
      select case(m_o)

        case(-1)     ! 2/3Ly + 2*Sy
          if( isp_g == isp_f ) then
            m_g_p_1 = m_cg + 1
            if( abs( m_g_p_1 ) <= L_g ) Cfac = - Two_third * p_plus * Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_g_p_1,.true.)
            m_g_m_1 = m_cg - 1
            if( abs( m_g_m_1 ) <= L_g ) Cfac = Cfac + Two_third * p_minus * Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_g_m_1,.true.)
          elseif( isp_g == 1 ) then
            Cfac = Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_cg,.true.)
          else
            Cfac = - Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_cg,.true.)
          endif
          Cfac = - 0.5_db * alfa_sf * Cfac
          
        case(0)     ! 2/3Lz + 2*Sz
          if( isp_g == isp_f ) then
            if( isp_g == 1 ) then
              Cfac = ( Two_third * m_cg + 1 ) * Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_cg,.true.)
            else
              Cfac = ( Two_third * m_cg - 1 ) * Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_cg,.true.)
            endif
          endif
          Cfac = img * 0.5_db * alfa_sf * Cfac

        case(1)     ! 2/3Lx + 2*Sx
          if( isp_g == isp_f ) then
            m_g_p_1 = m_cg + 1
            if( abs( m_g_p_1 ) <= L_g ) Cfac = Two_third * p_plus * Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_g_p_1,.true.)
            m_g_m_1 = m_cg - 1
            if( abs( m_g_m_1 ) <= L_g ) Cfac = Cfac + Two_third * p_minus * Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_g_m_1,.true.)
          else
            Cfac = Gaunt_xrc(L,m_cf,L_o,m_k,L_g,m_cg,.true.)
          endif
          Cfac = img * 0.5_db * alfa_sf * Cfac

      end select

! For the conversion from real to complex  
      if( .not. ( Ylm_comp .or. m == 0 ) ) then
        Cfac = Cfac / sqrt( 2._db )
        if(  m > 0 ) then   ! the (-1)^m
          if( i == 1 .and. mod(abs(m),2) == 1 ) Cfac = - Cfac
        else
          if( i == 2 .and. mod(abs(m),2) == 2 ) Cfac = - Cfac
  ! minus because in bra
          Cfac = - img * Cfac
        endif
      endif

      if( .not. ( Ylm_comp_g .or. m_g == 0 ) ) then
        Cfac = Cfac / sqrt( 2._db )
        if(  m_g > 0 ) then   ! the (-1)^m
          if( j == 1 .and. mod(abs(m_g),2) == 1 ) Cfac = - Cfac
        else
          if( j == 2 .and. mod(abs(m_g),2) == 2 ) Cfac = - Cfac
          Cfac = img * Cfac
        endif
      endif

      PseudoGaunt_M2 = PseudoGaunt_M2 + Cfac

      if( Ylm_comp .or. m == 0 ) exit

    end do

    if( Ylm_comp_g .or. m_g == 0 ) exit
    
  end do

  return
end

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

! Mat_out = Rot_cluster . Mat . Transpose(Rot_cluster)
! Rot_cluster is real

subroutine Rot_tensor_2( mat2, Rot_cluster )

  use declarations
  implicit none

  integer:: i, j, k
  complex(kind=db):: cmat
  complex(kind=db), dimension(3,3):: mat, mat2

  real(kind=db), dimension(3,3):: Rot_cluster

  mat(:,:) = (0._db, 0._db)
  do i = 1,3
    do j = 1,3
      do k = 1,3
        cmat = sum( Rot_cluster(j,:) * mat2(k,:) )
        mat(i,j) = mat(i,j) + Rot_cluster(i,k) * cmat
      end do
    end do
  end do
  mat2(:,:) = mat(:,:)

  return
end

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

subroutine Rot_tensor_3( mat3, Rot_cluster )

  use declarations
  implicit none

  integer:: i, j, k, l, m
  complex(kind=db):: cmas, cmat
  complex(kind=db), dimension(3,3,3):: mat, mat3

  real(kind=db), dimension(3,3):: Rot_cluster

  mat(:,:,:) = (0._db, 0._db)
  do i = 1,3
    do j = 1,3
      do k = 1,3
        do l = 1,3
          cmas = (0._db, 0._db)
          do m = 1,3
            cmat = sum( Rot_cluster(k,:) * mat3(l,m,:) )
            cmas = cmas + Rot_cluster(j,m) * cmat
          end do
          mat(i,j,k) = mat(i,j,k) + Rot_cluster(i,l) * cmas
        end do
      end do
    end do
  end do
  mat3(:,:,:) = mat(:,:,:)

  return
end

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

subroutine Rot_tensor_4( Mat4, Rot_cluster )

  use declarations
  implicit none

  integer:: i, j, k, l, m, n, n1
  complex(kind=db):: cmar, cmas, cmat
  complex(kind=db), dimension(3,3,3,3):: mat, mat4

  real(kind=db), dimension(3,3):: Rot_cluster

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

  do i = 1,3
    do j = 1,3
      do k = 1,3
        do l = 1,3
          do m = 1,3
            cmar = (0._db, 0._db)
            do n = 1,3
              cmas = (0._db, 0._db)
              do n1 = 1,3
                cmat = sum( Rot_cluster(l,:) * mat4(m,n,n1,:) )
                cmas = cmas + Rot_cluster(k,n1) * cmat
              end do
              cmar = cmar + Rot_cluster(j,n) * cmas
            end do
            mat(i,j,k,l) = mat(i,j,k,l) + Rot_cluster(i,m) * cmar
          end do
        end do
      end do
    end do
  end do

  mat4(:,:,:,:) = mat(:,:,:,:)

  return
end

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

subroutine Rot_tensor_6( Mat6, Rot_cluster )

  use declarations

  implicit none

  integer:: i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5

  complex(kind=db):: Cm1, Cm2, Cm3, Cm4, Cm5, Cm6
  complex(kind=db), dimension(3,3,3,3,3,3):: Mat, Mat6

  real(kind=db), dimension(3,3):: Rot_cluster

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

  do i1 = 1,3
    do i2 = 1,3
      do i3 = 1,3
        do i4 = 1,3
          do i5 = 1,3
            do i6 = 1,3

              Cm1 = (0._db, 0._db)
              do j1 = 1,3
                Cm2 = (0._db, 0._db)
                do j2 = 1,3
                  Cm3 = (0._db, 0._db)
                  do j3 = 1,3
                    Cm4 = (0._db, 0._db)
                    do j4 = 1,3
                      Cm5 = (0._db, 0._db)
                      do j5 = 1,3
                        Cm6 = sum( Rot_cluster(i6,:) * Mat6(j1,j2,j3,j4,j5,:) )
                        Cm5 = Cm5 + Rot_cluster(i5,j5) * Cm6
                      end do
                      Cm4 = Cm4 + Rot_cluster(i4,j4) * Cm5
                    end do
                    Cm3 = Cm3 + Rot_cluster(i3,j3) * Cm4
                  end do
                  Cm2 = Cm2 + Rot_cluster(i2,j2) * Cm3
                end do
                Cm1 = Cm1 + Rot_cluster(i1,j1) * Cm2
              end do

              Mat(i1,i2,i3,i4,i5,i6) = Mat(i1,i2,i3,i4,i5,i6) + Cm1

            end do
          end do
        end do
      end do
    end do
  end do

  Mat6(:,:,:,:,:,:) = Mat(:,:,:,:,:,:)

  return
end

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

! Calculation of S = Sum_if | <f| exp( iq.r ) |i> |^2
! exp(iqr) = 4 * pi * Sum_L ( i^l * j_l(qr) * Y_L^*(Omega_r) * Y_L^*(Omega_q) )

subroutine S_nrixs_cal( Chi, Coef_g, Core_resolved, Ecinetic, &
                  Eimag, Energ, Enervide, Eseuil, Final_tddft, Full_potential, Green, Green_int, Hubb_a, &
                  Hubb_d, icheck, L0_nrixs, Lmax_nrixs, is_g, Lmax, Lmax_pot, Lmoins1, Lplus1, Lseuil, m_g, &
                  m_hubb, mpinodes, mpirank, n_Ec, n_V, nb_test, nbseuil, ns_dipmag, ninit_tddft, &
                  ninit1, ninit, ninitr, ninitv, nlm_pot, nlm_probe, nlm_p_fp, nq_nrixs, nr, nrm, nspin, nspino, nspinp, &
                  Numat, Psii, q_nrixs, r, Relativiste, Renorm, Rmtg, Rmtsd, &
                  S_nrixs, S_nrixs_m, Solsing, Solsing_only, Spinorbite, Taull, Test_case, &
                  V_hubb, V_intmax, V0bd, Vrato, Ylm_comp )

  use declarations
  implicit none
  
  integer, parameter::  ip_max = 0
  integer, parameter::  ip0 = 0

  integer:: icheck, index_rank_i, index_rank_s, initr, initlt, iq, isp, k_i, k_s, L_i, L_s, L0_nrixs, Lmax, &
    Lmax_nrixs, Lmax_pot, Lm_i, Lm_s, Lseuil, m_hubb, m_i, m_k_i, m_k_s, m_s, mpinodes, mpirank, &
    n_Ec, n_V, nb_test, nbseuil, ninit_tddft, ninit1, ninit, ninitr, ninitv, nlm_pot, nlm_probe, &
    nlm_p_fp, nq_nrixs, nr, nrm, ns_dipmag, nspin, nspino, nspinp, Numat

  complex(kind=db):: dfac   
  complex(kind=db), dimension(ninitr):: Ten, Ten_m
  complex(kind=db), dimension(nlm_probe,nspinp,nlm_probe,nspinp):: Taull
  complex(kind=db), dimension(nlm_probe*nspino,nlm_probe*nspino,ninit_tddft,ninit_tddft,2,2,ns_dipmag):: Chi
  complex(kind=db), dimension(-m_hubb:m_hubb,nspinp,-m_hubb:m_hubb,nspinp):: V_hubb
  complex(kind=db), dimension(nq_nrixs,(Lmax_nrixs+1)**2,(Lmax_nrixs+1)**2,ninitr,0:mpinodes-1):: S_nrixs, S_nrixs_m
  
  complex(kind=db), dimension(:,:,:,:,:,:,:), allocatable:: Singul
  complex(kind=db), dimension(:,:,:,:,:,:), allocatable:: rof

  integer, dimension(ninit,2):: m_g
  integer, dimension(ninit):: is_g

  logical:: Core_resolved, Final_tddft, Full_potential, Green, Green_int, &
    Hubb_a, Hubb_d, Lmoins1, Lplus1, Monocrystal, NRIXS, Relativiste, Renorm, RIXS, Solsing, Solsing_only, Spinorbite, &
    SpinPos_i, SpinPos_s, Ylm_comp
    
  logical, dimension(nb_test):: Test_case

  real(kind=db):: Eimag, Energ, Enervide_t, Rmtg, Rmtsd, V_intmax
  real(kind=db), dimension(3):: q_vec
  real(kind=db), dimension(nspin):: Ecinetic_e, V0bd_e
  real(kind=db), dimension(nbseuil):: Eseuil, Omega
  real(kind=db), dimension(ninit,2):: Coef_g
  real(kind=db), dimension(nspin,n_Ec):: Ecinetic
  real(kind=db), dimension(n_Ec):: Enervide
  real(kind=db), dimension(nspin,n_V):: V0bd
  real(kind=db), dimension(nr):: r
  real(kind=db), dimension(nr,nlm_pot,nspin):: Vrato_e
  real(kind=db), dimension(nrm,nbseuil):: psii
  real(kind=db), dimension(nr,nlm_pot,nspin,n_V):: Vrato
  real(kind=db), dimension(nr,L0_nrixs:Lmax_nrixs):: bessel
  real(kind=db), dimension(4,nq_nrixs):: q_nrixs

  NRIXS = .true.
  RIXS = .false.
  SpinPos_i = .false. !No spin position in Raman X
  SpinPos_s = .false.
  k_i = 0   ! not used
  k_s = 0   ! not used
  Omega(:) = 0._db ! not used
  m_k_i = 0; m_k_s = 0 ! .not.used

  allocate( Singul(nlm_probe,nspinp,nlm_probe,nspinp,L0_nrixs:Lmax_nrixs,L0_nrixs:Lmax_nrixs,ninitv) )
  allocate( rof(nlm_probe,nlm_p_fp,nspinp,nspino,L0_nrixs:Lmax_nrixs,ninitv) )
  S_nrixs(:,:,:,:,mpirank) = (0._db,0._db)
  S_nrixs_m(:,:,:,:,mpirank) = (0._db,0._db)

  do iq = 1,nq_nrixs

    if( icheck > 1 ) write(3,100) iq, q_nrixs(4,iq) / bohr
    
    rof(:,:,:,:,:,:) = (0._db, 0._db)
    Singul(:,:,:,:,:,:,:) = (0._db, 0._db)

    call cbessel(bessel,L0_nrixs,Lmax_nrixs,nr,q_nrixs(4,iq),r)
  
    do initlt = 1,n_Ec
      if( Final_tddft ) then
        do isp = 1,nspin
          Ecinetic_e(isp) = Ecinetic(isp,initlt)
          V0bd_e(isp) = V0bd(isp,initlt)
          Vrato_e(1:nr,:,isp) = Vrato(1:nr,:,isp,initlt)
        end do
      else
        do isp = 1,nspin
          Ecinetic_e(isp) = Ecinetic(isp,1)
          V0bd_e(isp) = V0bd(isp,1)
          Vrato_e(1:nr,:,isp) = Vrato(1:nr,:,isp,1)
        end do
      endif
      Enervide_t = Enervide(initlt)

      call Radial_integral(Ecinetic_e,Eimag,Energ,Enervide_t,Eseuil,Final_tddft,Full_potential,Hubb_a,Hubb_d,icheck, &
        initlt,Lmax_nrixs,L0_nrixs,Lmax,Lmax_pot,m_hubb,nb_test,nbseuil,ninit1,ninitv,nlm_pot,nlm_probe,nlm_p_fp,nr,NRIXS,nrm, &
        nspin,nspino,nspinp,Numat,psii,r,bessel,Relativiste,Renorm,Rmtg,Rmtsd,rof,Singul,Solsing,Spinorbite,Test_case, &
        V_hubb,V_intmax,V0bd_e,Vrato_e,Ylm_comp)

    end do

    if( icheck > 1 ) write(3,110) 
    
    q_vec(:) = q_nrixs(1:3,iq)
    Monocrystal = sum( abs(q_vec(:)) ) > eps10

    Lm_s = 0   
    do L_s = L0_nrixs,Lmax_nrixs
      index_rank_s = L_s
      
      do m_s = -L_s,L_s
        Lm_s = Lm_s + 1

        Lm_i = 0
        do L_i = L0_nrixs,Lmax_nrixs
          index_rank_i = L_i
          
          if( L_s == L_i ) then
            dfac = cmplx( quatre_pi**2, 0._db, db )
          else 
            dfac = img**(L_s - L_i) * quatre_pi**2
          endif

          do m_i = -L_i,L_i
            Lm_i = Lm_i + 1
    
            if( ( L_s /= L_i .or. m_s /= m_i ) .and. .not. Monocrystal ) cycle

            if( icheck > 1 ) write(3,150) L_s, m_s, L_i, m_i

            call Tens_ab( Chi, Coef_g, Core_resolved, Final_tddft, Full_potential, &
                                Green, Green_int, icheck, index_rank_i, index_rank_s, Lmax_nrixs, L0_nrixs, is_g, k_i, k_s, L_s, &
                                m_k_i, m_k_s, m_s, L_i, m_g, m_i, Lmax, Lmoins1, Lplus1, Lseuil, nbseuil, ns_dipmag, ninit1, &
                                ninit, ninit_tddft, ninitv, ninitr, nlm_probe, nlm_p_fp, NRIXS, nspinp, nspino, Omega, rof, &
                                Singul, Solsing, Solsing_only, Spinorbite, SpinPos_i, SpinPos_s, Taull, Ten, Ten_m, Ylm_comp )

            S_nrixs(iq,Lm_s,Lm_i,:,mpirank) = Ten(:) * dfac
            if( Green_int ) S_nrixs_m(iq,Lm_s,Lm_i,:,mpirank) = Ten_m(:) * dfac

            if( icheck > 1 ) write(3,160) L_s, m_s, L_i, m_i, dfac, &
                                         ( S_nrixs(iq,Lm_s,Lm_i,initr,mpirank), initr = 1,ninitr )  

          end do
        end do
      end do
    end do ! end of loop over L_s

  end do

  deallocate( rof, Singul )
  
  return
  100 format(/' ---- S_nrixs_cal Radial ---------',100('-'),//'  iq =',i3,', q =',f6.3)
  110 format(/' ---- S_nrixs_cal Tens_ab --------',100('-'))
  120 format(/'  iq =',i3,', q =',f6.3,', q_vec =',3f10.5)
  130 format(/' Normalized q_vec in the absorbing atom local basis =',3f10.5)
  140 format(/'  iq =',i3,', q =',f6.3)
  150 format(/' L_s, m_s =',2i3,',  L_i, m_i =',2i3)
  160 format(/' L_s, m_s =',2i3,',  L_i, m_i =',2i3,', dfac =',1p,2e15.7,', S_nrixs(1..ninitr) =',28e15.7)

end


