! Subroutines for FDMNES
! Contains different tools

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

subroutine ad_number(ib,nomfich,Length)

  implicit none

  integer:: i, ib, in, ipuis, iu, iumax, L, Length

  character(len=Length):: nomfich

  L = len_trim(nomfich)

  if( ib < 0 ) then
    L = L + 1
    if( L <= Length ) nomfich(L:L) = '-'
  endif

  i = abs(ib)

  do iu = 2,10
    if( i / 10**(iu-1) < 1 ) exit
  end do
  iumax = iu - 1

  do iu = iumax,1,-1

    ipuis = 10**(iu-1)

    in = i / ipuis

! When toolong one writes only the last numbers
    if( L + iu <= Length ) then
      L = L + 1
      nomfich(L:L) = achar(in+48)
    elseif( L + iu == Length-1 ) then
      if( nomfich(L:L) == '_' ) nomfich(L:L) = achar(in+48)
    endif

    i = i - ipuis * in

  end do

  return
end

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

subroutine ad_number_r(r,nom)

  use declarations
  implicit none

  integer:: i, l, l2, l3

  character(len=1) a
  character(len=Length_word) mot, nom

  real(kind=db):: r
  
  if( Length_word >= 15 ) then 
    write(mot,'(f15.4)') r
  elseif( Length_word >= 13 ) then 
    write(mot,'(f13.3)') r
  elseif( Length_word >= 10 ) then 
    write(mot,'(f10.2)') r
  else
    write(mot,'(f5.2)') r
  endif
  
  mot = adjustl( mot )
  
  l2 = len_trim(mot) 
  l3 = l2

  do i = l3,2,-1
    a = mot(i:i)
    if( a /= '0' .and. a /= '.' ) exit
    mot(i:i) = ' '
    l2 = l2 - 1
    if( a == '.' ) exit
  end do
 
  l = len_trim( nom )

  do i = l+1,min(l+l2,Length_word) 
    nom(i:i) = mot(i-l:i-l)
  end do
      
  return
end

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

subroutine ad_sufix(Core_resolved,initr,jseuil,Length,mot,nseuil)

  use declarations
  implicit none

  integer:: initr, jseuil, l, Length, nseuil
    
  character(len=Length):: mot
  
  Logical:: Core_resolved 

  l = len_trim(mot)

  if( Core_resolved ) then
    if( l < Length - 2 ) mot(l+1:l+1) = '_'
    call ad_number(initr,mot,Length)
  else  
    if( l < Length - 3 ) then
      l = l + 1
      mot(l:l) = '_'
    endif
    if( l < Length - 2 ) then
      l = l + 1
      mot(l:l) = achar(nseuil+74)
    endif
    call ad_number(jseuil+initr-1,mot,Length)
  endif

  return

end

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

subroutine center_word( mot, Length_word )

  implicit none
  
  integer:: L, Length_word, Lm, Lshift
  
  character(len=*) mot
  character(len=Length_word) mot2

  mot2 = ' '
  mot = adjustl( mot )
  L = len_trim( mot )
  Lshift = ( Length_word - L + 1 ) / 2
  Lshift = max( 0, lshift )
  Lm = min( L, Length_word )
  mot2(1+Lshift:Lm+Lshift) = mot(1:Lm)
  mot = mot2

  return
end

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

! Fonction giving the number of number in the next non empty line et setting at the biginning of this line.
! With line starting with a character, nnombre = 0

function nnombre(irec,Length)

  use declarations
  implicit none

  integer, parameter:: nlet = 53

  integer:: eof, i, icol, irec, j, k, l, Length, ligne, nmots, n, nnombre

  character(len=Length):: mot, test
  character(len=1), dimension(nlet):: let

  real(kind=db):: x

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

  do i = 1,Length
    mot(i:i) = 'P'
  end do

  nmots = 0

  boucle_ligne: do ligne = 1,1000000

    read(irec,'(A)',iostat=eof) mot

! An "END=" condition leaves the file position after the endfile record,
! thus backspace to be just before the endfile record. This way a
! subsequent READ with another "END=" clause can execute correctly.
    if( eof /= 0 ) then
      nnombre = 0
      backspace(irec)
      return
    endif

    do i = 1,Length
      if( mot(i:i) /= ' ' .and. mot(i:i) /= char(9) ) exit boucle_ligne
    end do

  end do boucle_ligne

  backspace(irec)

  Open( 16, status='SCRATCH' )

  do i = 1,Length
    if( mot(i:i) == char(9) ) mot(i:i) = ' '
  end do

  n = 0
  i = 0
  do icol = 1,Length
    i = i + 1
    if( i > Length ) exit
    if( mot(i:i) == ' ' ) cycle
    do j = i+1,Length
      if( mot(j:j) == ' ' ) exit
    end do
    j = j - 1
    n = n + 1
    write(16,'(A)') mot(i:j)
    i = j
  end do

  Rewind(16)

  boucle_i: do i = 1,n
    read(16,*,iostat=eof) x
    if( eof /= 0 ) exit boucle_i
    backspace(16)
    read(16,'(A)',iostat=eof) test
    l = len_trim(test)
    do j = 1,l
      do k = 5,nlet
        if( test(j:j) == let(k) ) exit boucle_i
      end do
    end do
    do k = 1,4
      if( test(1:1) == let(k) .or. test(l:l)==let(k) ) exit boucle_i
    end do
    nmots = nmots + 1
  end do boucle_i
  Close(16)


  nnombre = nmots
  return

end

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

! Column name of type "Name(L,m,isp)" for n_index = 3
!                     "Name(L,isp)" for n_index = 2 and nspin = 2 
!                     "Name(L,m)" for n_index = 2 and nspin = 1
!                     "Name(L)" for n_index = 1
!             L rage from L_0 to L_max
!             m range from -L up to +L if m_max < m_0
!             m range from m_0 up to m_max in the other case

subroutine Create_col_name(Col_name,L_0,L_max,m_0,m_max,Length_col,Length_n,Name,n_index,n_dim,nspin)

  use declarations
  implicit none
  
  integer:: i, L_0, L_max, isp, m, m_0, m_1, m_2, m_max, k, Length, Length_col, Length_n, n_index, n_dim, nspin
  
  character(len=Length_n):: Name, Name_b 
  character(len=Length_col):: Word 
  character(len=Length_col), dimension(n_dim):: Col_name 

  k = 0

  Col_name(:) = ' '
    
  Name_b = adjustl(Name)

  do isp = 1,nspin

    do i = L_0, L_max    
      if( n_index == 1 .and. nspin == 2 .and. i > L_0 ) exit 

      if( m_max < m_0 ) then
        m_1 = - i; m_2 = i
      else
        m_1 = m_0; m_2 = m_max
      endif
       
      do m = m_1, m_2
        if( n_index == 2 .and. nspin == 2 .and. m > m_0 ) exit 

        k = k + 1
          
        Word = ' '
        Length = Len_trim(Name_b)
        Length = min( Length, Length_col-2 )
        Word(1:Length) = Name(1:Length)
        
        Length = len_trim(Word) + 1
        if( Length < Length_col - 2 ) Word(Length:Length) = '('
        
        if( n_index > 1 .or. nspin == 1 ) then
          call ad_number(i,Word,Length_col)
          Length = len_trim(Word) + 1
        endif

        if( n_index > 2 .or. ( n_index == 2 .and. nspin == 1 ) ) then
          if( Length < Length_col-1) Word(Length:Length) = ','
          call ad_number(m,Word,Length_col)
          Length = len_trim(Word) + 1
        endif

        if( n_index == 3 .or. nspin == 2 ) then
          if( Length+1 < Length_col ) then
            if( isp == 1 ) then
              Word(Length:Length+1) = ',u'
            else
              Word(Length:Length+1) = ',d'
            endif
          elseif( Length < Length_col ) then
            if( isp == 1 ) then
              Word(Length:Length+1) = 'u'
            else
              Word(Length:Length+1) = 'd'
            endif
          endif
          Length = len_trim(Word) + 1
        endif

        if( Length < Length_col )  Word(Length:Length) = ')'
        
        call center_word( Word, Length_col ) 
        Col_name(k) = Word

      end do
    end do
  end do
    
  return
end
  
!***********************************************************************

! Calculation of the integral of f(r), from 0 to Radius
! Integral is calculated using a second order polynomia
! nr: number of points

  function f_integr(r,fct,nr,ir0,nrm,Radius)

  use declarations
  implicit none

  integer:: i, ir0, nr, nrm

  real(kind=db):: a, b, c, d0m, dp0, dpm, f_integr, f0, fm, fp, r0, Radius, rm, rp, tiers, x1, x2
  real(kind=db), dimension(ir0:nrm):: fct, r

  tiers = 1._db / 3

  f_integr = 0._db

  do i = 2,nr-1
    rm = r(i-1)
    r0 = r(i)
    rp = r(i+1)
    if( i == 2 ) then
      x1 = 0._db
    else
      x1 = 0.5_db * ( rm + r0 )
    endif
    if( rp > Radius ) then
      x2 = Radius
    else
      if( i == nr - 1 ) then
        x2 = rp
      else
        x2 = 0.5_db * ( rp + r0 )
      endif
      fm = fct(i-1)
      f0 = fct(i)
      fp = fct(i+1)
      if( abs(fp) < 1e-20_db .and. abs(f0) < 1e-20_db .and. abs(fm) < 1e-20_db ) cycle

      dp0 = rp - r0
      dpm = rp - rm
      d0m = r0 - rm
      a = ( fm * dp0 - f0 * dpm + fp * d0m ) / ( d0m * dp0 * dpm )
      b = ( f0 - fm ) / d0m - a * ( r0 + rm )
      c = f0 - a * r0**2 - b * r0
    endif

    f_integr = f_integr + ( tiers * a * ( x1**2 + x1 * x2  + x2**2 ) + 0.5_db * b * ( x1 + x2 ) + c ) * ( x2 - x1 )

    if( rp > Radius ) exit

  end do

  return
end

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

! Calcul de L'integrale de 0 a Radius de f.
! L'integrale est calculee avec un polynome d'interpolation d'ordre 3
! r: contient les valeurs des rayons pour les points qu'on integre
! On considere que seule les valeurs jusqu'à L'indice n tel que
! r(i) > Radius > Radius(i-1) sont sures.

  real(kind=db) function f_integr3(r,fct,ir0,nrm,Radius)

  use declarations
  implicit none

  integer:: i, ipr, ir0, nrm

  logical:: This_is_the_end

  real(kind=db):: a, b, c, d, d0m, dp0, dpm, f0, f1, f2, f3, f4, fac1, fac2, fm, fp, r0, Radius, r1, r2, r3, r4, rap14, rap24, &
                  rap34, rm, rp, Tiers
  real(kind=db), dimension(ir0:nrm):: fct, r

  f_integr3 = 0._db

  if( Radius < eps10 ) return

  if( nrm < 4 ) then
    call write_error
    do ipr = 6,9,3
      write(ipr,110) nrm
    end do
    stop
  endif

  Tiers = 1._db / 3
  This_is_the_end = .false.

  do i = 1,nrm-1
    if( i == 1 ) then
      if( ir0 == 1 ) then
 ! Construction de L'ordonnée par extrapolation du second degre
        rm = r(1)
        r0 = r(2)
        rp = r(3)
        fm = fct(1)
        f0 = fct(2)
        fp = fct(3)
        if( abs(fp) < 1e-20_db .and. abs(f0) < 1e-20_db .and. abs(fm) < 1e-20_db ) cycle
        dp0 = rp - r0
        dpm = rp - rm
        d0m = r0 - rm
        a = ( fm * dp0 - f0 * dpm + fp * d0m ) / ( d0m * dp0 * dpm )
        b = ( f0 - fm ) / d0m - a * ( r0 + rm )
        c = f0 - a * r0**2 - b * r0
        f1 = c
      else
        f1 = fct(i-1)
      endif
      r1 = 0._db
      r2 = r(i)
      f2 = fct(i)
      r3 = r(i+1)
      f3 = fct(i+1)
      r4 = r(i+2)
      f4 = fct(i+2)
      rm = 0._db
      rp = r3
    elseif( i == nrm-1 .or. r(i+1) > Radius ) then
      rm = r(i)
      rp = r(nrm)
      if( i < 3 ) then
        r1 = r(i-1)
        f1 = fct(i-1)
        r2 = r(i)
        f2 = fct(i)
        r3 = r(i+1)
        f3 = fct(i+1)
        r4 = r(i+2)
        f4 = fct(i+2)
      endif
    else
      r1 = r(i-1)
      f1 = fct(i-1)
      r2 = r(i)
      f2 = fct(i)
      r3 = r(i+1)
      f3 = fct(i+1)
      r4 = r(i+2)
      f4 = fct(i+2)
      rm = r2
      rp = r3
    endif
    if( rp > Radius ) then
      rp = Radius
      This_is_the_end = .true.
    endif
    if( abs(f2) < 1e-20_db .and. abs(f3) < 1e-20_db ) cycle

    rap14 = (f1 - f4) / (r1 - r4)
    rap24 = (f2 - f4) / (r2 - r4)
    rap34 = (f3 - f4) / (r3 - r4)
    fac1 = ( rap14 - rap34 ) / (r1 - r3)
    fac2 = ( rap24 - rap34 ) / (r2 - r3)
    a = ( fac1 - fac2 ) / ( r1 - r2 )
    b = fac1 - a * ( r1 + r3 + r4 )
    c = rap14 - a * ( r1**2 + r1 * r4 + r4**2 ) - b * ( r1 + r4 )
    d = f1 - ( a * r1**2 + b * r1 + c ) * r1

    f_integr3 = f_integr3 + ( 0.25_db * ( rm + rp ) * ( rm**2 + rp**2 ) * a + tiers * ( rm**2 + rm*rp + rp**2 ) * b &
              + 0.5_db * ( rm + rp ) * c + d ) * ( rp - rm )

    if( This_is_the_end ) exit

  end do

  return
  110 format(/' Call to function f_integr3 needs a number of point higher than 3, here it is =',i2,' !')
end

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

! Calcul de L'integrale de 0 a Radius de f. f(0) = 0.
! L'integrale est calculee avec un polynome d'interpolation d'ordre 3
! r: contient les valeurs des rayons pour les points qu'on integre
! nr: nombre de points qu'on integre

  real(kind=db) function f_integr32(r,fct,ir0,nr,R0,Radius)

  use declarations
  implicit none

  integer:: i, ir0, ir1, nr

  logical This_is_the_end

  real(kind=db):: a, b, c, d, f1, f2, f3, f4, fac1, fac2, R0, r1, r2, r3, r4, rap14, rap24, rap34, rm, rp, Radius, tiers
  real(kind=db), dimension(ir0:nr):: fct, r

  tiers = 1._db / 3
  This_is_the_end = .false.
  f_integr32 = 0._db
  ir1 = ir0 + 1

  do i = ir1,nr-2
    if( r(i) < R0 ) cycle
    r1 = r(i-1)
    f1 = fct(i-1)
    r2 = r(i)
    f2 = fct(i)
    r3 = r(i+1)
    f3 = fct(i+1)
    r4 = r(i+2)
    f4 = fct(i+2)
    if( i == ir1 .or. r1 < R0 ) then
      rm = R0
    else
      rm = r2
    endif
    if( i == nr-2 .or. r4 > Radius ) then
      rp = Radius
      This_is_the_end = .true.
    else
      rp = r3
    endif

    if( abs(f2) < 1e-20_db .and. abs(f3) < 1e-20_db ) cycle

    rap14 = (f1 - f4) / (r1 - r4)
    rap24 = (f2 - f4) / (r2 - r4)
    rap34 = (f3 - f4) / (r3 - r4)
    fac1 = ( rap14 - rap34 ) / (r1 - r3)
    fac2 = ( rap24 - rap34 ) / (r2 - r3)
    a = ( fac1 - fac2 ) / ( r1 - r2 )
    b = fac1 - a * ( r1 + r3 + r4 )
    c = rap14 - a * ( r1**2 + r1 * r4 + r4**2 ) - b * ( r1 + r4 )
    d = f1 - ( a * r1**2 + b * r1 + c ) * r1

    f_integr32 = f_integr32 + ( 0.25_db * ( rm + rp ) * ( rm**2 + rp**2 ) * a + tiers * ( rm**2 + rm*rp + rp**2 ) * b &
               + 0.5_db * ( rm + rp ) * c + d ) * ( rp - rm )

    if( This_is_the_end ) exit

  end do

  return
end

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

  real(kind=db) function f_interp1(x,x1,x2,y1,y2)

  use declarations
  implicit none

  real(kind=db),intent(in):: x, x1, x2, y1, y2
  real(kind=db):: p

  p = (x - x1) / (x2 - x1)
  f_interp1 = p*y2 + (1-p)*y1

  return
end

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

  real(kind=db) function f_interp2(r,rm,r0,rp,fm,f0,fp)

  use declarations
  implicit none

  real(kind=db):: a, b, c, d0m, dp0, dpm, f0, fm, fp, r0, rm, rp, r

  dp0 = rp - r0
  dpm = rp - rm
  d0m = r0 - rm

  a = ( fm * dp0 - f0 * dpm + fp * d0m ) / ( d0m * dp0 * dpm )
  b = ( f0 - fm ) / d0m - a * ( r0 + rm )
  c = f0 - a * r0**2 - b * r0

  f_interp2 = a * r**2 + b * r + c

  return
end

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

  function f_interp2_cp(r,rm,r0,rp,fm,f0,fp)

  use declarations
  implicit none

  complex(kind=db):: a, b, c, f_interp2_cp, f0, fm, fp
  real(kind=db):: d0m, dp0, dpm, r0, rm, rp, r

  dp0 = rp - r0
  dpm = rp - rm
  d0m = r0 - rm

  a = ( fm * dp0 - f0 * dpm + fp * d0m ) / ( d0m * dp0 * dpm )
  b = ( f0 - fm ) / d0m - a * ( r0 + rm )
  c = f0 - a * r0**2 - b * r0

  f_interp2_cp = a * r**2 + b * r + c

  return
end

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

! Interpolation d'ordre 3, pour obtenir fonction et derivee

subroutine interp3(f,fp,r,r1,r2,r3,r4,f1,f2,f3,f4)

  use declarations
  implicit none

  real(kind=db):: a, b, c, d, f, f1, f2, f3, f4, fac1, fac2, fp, r, r1, r2, r3, r4, rap14, rap24, rap34

  rap14 = (f1 - f4) / (r1 - r4)
  rap24 = (f2 - f4) / (r2 - r4)
  rap34 = (f3 - f4) / (r3 - r4)
  fac1 = ( rap14 - rap34 ) / (r1 - r3)
  fac2 = ( rap24 - rap34 ) / (r2 - r3)
  a = ( fac1 - fac2 ) / ( r1 - r2 )
  b = fac1 - a * ( r1 + r3 + r4 )
  c = rap14 - a * ( r1**2 + r1 * r4 + r4**2 ) - b * ( r1 + r4 )
  d = f1 - ( a * r1**2 + b * r1 + c ) * r1

  f = a * r**3 + b * r**2 + c * r + d
  fp = 3 * a * r**2 + 2 * b * r + c

  return
end

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

! Interpolation d'ordre 3, pour obtenir fonction et derivee

subroutine interp3_pp(f,fp,fpp,r,r1,r2,r3,r4,f1,f2,f3,f4)

  use declarations
  implicit none

  real(kind=db):: a, b, c, d, f, f1, f2, f3, f4, fac1, fac2, fp, fpp, r, r1, r2, r3, r4, rap14, rap24, rap34

  rap14 = (f1 - f4) / (r1 - r4)
  rap24 = (f2 - f4) / (r2 - r4)
  rap34 = (f3 - f4) / (r3 - r4)
  fac1 = ( rap14 - rap34 ) / (r1 - r3)
  fac2 = ( rap24 - rap34 ) / (r2 - r3)
  a = ( fac1 - fac2 ) / ( r1 - r2 )
  b = fac1 - a * ( r1 + r3 + r4 )
  c = rap14 - a * ( r1**2 + r1 * r4 + r4**2 ) - b * ( r1 + r4 )
  d = f1 - ( a * r1**2 + b * r1 + c ) * r1

  f = a * r**3 + b * r**2 + c * r + d
  fp = 3 * a * r**2 + 2 * b * r + c
  fpp = 6 * a * r + 2 * b

  return
end

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

! Interpolation d'ordre 3

  real(kind=db) function f_interp3(r,r1,r2,r3,r4,f1,f2,f3,f4)

  use declarations
  implicit none

  real(kind=db):: a, b, c, d, f1, f2, f3, f4, fac1, fac2, r, r1, r2, r3, r4, rap14, rap24, rap34

  rap14 = (f1 - f4) / (r1 - r4)
  rap24 = (f2 - f4) / (r2 - r4)
  rap34 = (f3 - f4) / (r3 - r4)
  fac1 = ( rap14 - rap34 ) / (r1 - r3)
  fac2 = ( rap24 - rap34 ) / (r2 - r3)
  a = ( fac1 - fac2 ) / ( r1 - r2 )
  b = fac1 - a * ( r1 + r3 + r4 )
  c = rap14 - a * ( r1**2 + r1 * r4 + r4**2 ) - b * ( r1 + r4 )
  d = f1 - ( a * r1**2 + b * r1 + c ) * r1

  f_interp3 = a * r**3 + b * r**2 + c * r + d

  return
end

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

! Calculation of the derivative of a function
! Second order polynomia approximation
  
  subroutine Derivation(df,f,n,r)

  use declarations
  implicit none
  
  integer:: i, n
  
  real(kind=db):: hm, hp
  real(kind=db):: df(n), f(n), r(n)

  if( n < 2 ) then
    if( n == 1 ) df(1) = 0._db
    return
  endif
  
  df(1) = ( f(2) - f(1) ) / ( r(2) - r(1) )
  
  do i = 2,n-1
    hp = r(i+1) - r(i)
    hm = r(i) - r(i-1)
    df(i) = ( ( f(i+1) - f(i) ) * hm / hp - ( f(i-1) - f(i) ) * hp / hm ) / ( hm + hp )
  end do

  df(n) = ( f(n) - f(n-1) ) / ( r(n) - r(n-1) )
     
  return
end

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

subroutine mult_cell(itape,File_out)

  use declarations
  implicit none

  integer, parameter:: nam = 100000
  integer, parameter:: ntypem = 100

  integer:: eof, eoff, i, igrdat, ipr, itape, ix, iy, iz, j, l, n, na, nnombre, ntype

  real(kind=db):: Thickness
  real(kind=db), dimension(3):: A, abc, abc_new, Angle, Angle_new, B, C, CosAngle, q, V
  real(kind=db), dimension(3,3):: Mat, Mat_i
  real(kind=db), dimension(3,nam):: p

  integer, dimension(3):: nm
  integer, dimension(nam):: itype, Z
  integer, dimension(ntypem):: Numat

  Logical:: Ang, Mat_mul, Surface, Tr_cub_hexa, Typ

  character(len=2):: Chemical_Symbol
  character(len=9):: keyword
  character(len=132):: identmot, mot
  character(len=Length_name):: File_out

  Ang = .false.
  Mat_mul = .false.
  Tr_cub_hexa = .false.
  Surface = .false.
  Typ = .false.
  Thickness = 0._db
  Numat(:) = 0
  nm(:) = 1

  l = len_trim(File_out)
  if( File_out(l-3:l) /= '.txt' ) File_out(l+1:l+4) = '.txt'
  open(2,file = File_out)

! Lecture

  Rewind(itape)

  do igrdat = 1,100000

    read(itape,'(A)',iostat=eoff) mot
    if( eoff /= 0 ) exit

    keyword = identmot(mot,9)
    if( keyword(1:1) == '!' ) cycle

    select case(keyword)

      case('end')
        exit

      case('mult_cell')
        n = nnombre(itape,132)
        if( n == 0 ) call write_err_form(itape,keyword)
        n = min(n,3)
        read(itape,*) nm(1:n)

      case('mat_mul')
        Mat_mul = .true.
        do i  = 1,3
          n = nnombre(itape,132)
          if( n < 3 ) call write_err_form(itape,keyword)
          read(itape,*) Mat(:,i)
        end do

      case('cub_hexa')
        Tr_cub_hexa = .true.
        
      case('surf_cell')
        Surface = .true.
        n = nnombre(itape,132)
        if( n == 0 ) call write_err_form(itape,keyword)
        n = min(n,3)
        if( n == 1 ) then
          read(itape,*) Thickness
        else
          read(itape,*) nm(1:2), Thickness
        endif

      case('atomic_nu')
        Typ = .true.
        ntype = nnombre(itape,132)
        if( ntype > ntypem ) then
          call write_error
            do ipr = 6,9,3
              write(ipr,110) ntype, ntypem
            end do
          stop
        endif
        read(itape,*,iostat=eoff)  Numat(1:ntype)
        if( eoff /= 0 ) call write_err_form(itape,keyword)

      case('unit_cell')
        n = nnombre(itape,132)
        if( n == 3 .or. n == 4 .or. n == 5 ) then
          read(itape,*) abc(:)
        elseif( n >= 6 ) then
          Ang = .true.
          read(itape,*) abc(:), Angle(:)
        else
          call write_err_form(itape,keyword)
        endif

        i = 0
        do j = 1,100000
          if( j <= nam ) i = i + 1
          read(itape,*,iostat=eof) itype(i), p(:,i)
          if( eof /= 0 ) exit
          if( Typ ) then
            Z(i) = Numat(itype(i))
          else
            Z(i) = itype(i)
          endif
        end do

! ESRF changes start
! An "END=" condition leaves the file position after the endfile record,
! thus backspace to be just before the endfile record. This way a
! subsequent READ with another "END=" clause can execute correctly.
        backspace(itape)
! ESRF changes end

        if( j > nam ) then
          call write_error
          do ipr = 6,9,3
            write(ipr,130) j, nam
          end do
          stop
        endif

        na = i-1

      case default

        if( keyword(1:1) /= ' ' ) then
          call write_error
          do ipr = 6,9,3
            write(ipr,140)
            write(ipr,150) mot
          end do
          stop
        endif

    end select

  end do

  if( Mat_mul ) then
    do i = 1,3
      A(i) = Mat(1,i) * abc(i) 
      B(i) = Mat(2,i) * abc(i) 
      C(i) = Mat(3,i) * abc(i)
    end do
    if( Ang ) then
      CosAngle(:) = cos( Angle(:) * pi / 180._db )
      abc_new(1) = sqrt( sum( A(:)**2 ) + 2 * A(1)*A(2) * cosAngle(3) + 2 * A(1)*A(3) * cosAngle(2) &
                                        + 2 * A(2)*A(3) * cosAngle(1) ) 
      abc_new(2) = sqrt( sum( B(:)**2 ) + 2 * B(1)*B(2) * cosAngle(3) + 2 * B(1)*B(3) * cosAngle(2) &
                                        + 2 * B(2)*B(3) * cosAngle(1) ) 
      abc_new(3) = sqrt( sum( C(:)**2 ) + 2 * C(1)*C(2) * cosAngle(3) + 2 * C(1)*C(3) * cosAngle(2) &
                                        + 2 * C(2)*C(3) * cosAngle(1) )
      Angle_new(1) = sum( B(:)*C(:) ) + ( B(1)*C(2) + B(2)*C(1) ) * cosAngle(3) + ( B(1)*C(3) + B(3)*C(1) ) * cosAngle(2) &
                                      + ( B(2)*C(3) + B(3)*C(2) ) * cosAngle(1)
      Angle_new(1) = Angle_new(1) / ( abc_new(2) * abc_new(3) )                                           

      Angle_new(2) = sum( A(:)*C(:) ) + ( A(1)*C(2) + A(2)*C(1) ) * cosAngle(3) + ( A(1)*C(3) + A(3)*C(1) ) * cosAngle(2) &
                                      + ( A(2)*C(3) + A(3)*C(2) ) * cosAngle(1)
      Angle_new(2) = Angle_new(2) / ( abc_new(1) * abc_new(3) )                                           

      Angle_new(3) = sum( A(:)*B(:) ) + ( A(1)*B(2) + A(2)*B(1) ) * cosAngle(3) + ( A(1)*B(3) + A(3)*B(1) ) * cosAngle(2) &
                                      + ( A(2)*B(3) + A(3)*B(2) ) * cosAngle(1)
      Angle_new(3) = Angle_new(3) / ( abc_new(1) * abc_new(2) )
      
      Angle_new(:) = acos( Angle_new(:) ) * 180._db / pi                                           
    else
      abc_new(1) = sqrt( sum( A(:)**2 ) ) 
      abc_new(2) = sqrt( sum( B(:)**2 ) )
      abc_new(3) = sqrt( sum( C(:)**2 ) )
    endif
  
  elseif( Tr_cub_hexa ) then
! cubic to hexa
    abc_new(1:2) = abc(1:2) * sqrt(2._db)
    abc_new(3) = abc(3) * sqrt(3._db)
    Angle_new(1:2) = Angle(1:2); 
    Angle_new(3) = 120._db  

    A(1) =  1._db; A(2) = 0._db; A(3) =-1._db
    B(1) = -1._db; B(2) = 1._db; B(3) = 0._db
    C(1:3) = 1._db
    Mat(:,1) = A(:)
    Mat(:,2) = B(:)
    Mat(:,3) = C(:)

  else
  
    if( Surface ) nm(3) = 1
    abc_new(:) = abc(:) * nm(:)
    Angle_new(:) = Angle(:)
 
  endif
  
  if( Ang ) then
    write(2,120) abc_new(:), Angle_new(:)
    write(6,120) abc_new(:), Angle_new(:)
  else
    write(2,120) abc_new(:)
    write(6,120) abc_new(:)
  endif
    
  if( Mat_mul .or. Tr_cub_hexa ) then

    call invermat(Mat,Mat_i)
    
    j = 0
    do i = 1,na

      do ix = -10,10
        v(1) = p(1,i) + ix
        do iy = -10,10
          v(2) = p(2,i) + iy
          do iz = -10,10
            v(3) = p(3,i) + iz
          
            q = Matmul( Mat_i, v )
            
            if( q(1) > 1._db - eps10 .or. q(2) > 1._db - eps10 .or. q(3) > 1._db - eps10 &
             .or. q(1) < - eps10 .or. q(2) < - eps10 .or. q(3) < - eps10 ) cycle
         
            j = j + 1
               
            write(2,160) itype(i), q(:), j, Chemical_Symbol(Z(i))
            write(6,160) itype(i), q(:), j, Chemical_Symbol(Z(i))
            
          end do
        end do
      end do
    end do
    
  else
  
    if( Surface ) nm(3) = int( Thickness / abc(3) ) + 1
    j = 0
    do ix = 0,nm(1)-1
      do iy = 0,nm(2)-1
        do iz = 0,nm(3)-1
          write(2,*)
          do i = 1,na
            j = j + 1
            q(1) = ( p(1,i) + ix ) / nm(1)
            q(2) = ( p(2,i) + iy ) / nm(2)
            if( Surface ) then
              q(3) = p(3,i) + iz
            else
              q(3) = ( p(3,i) + iz ) / nm(3)
            endif
            if( Surface .and. q(3) > Thickness ) cycle
            write(2,160) itype(i), q(:), j, Chemical_Symbol(Z(i))
            write(6,160) itype(i), q(:), j, Chemical_Symbol(Z(i))
          end do

        end do
      end do
    end do

  endif
  
  Close(itape)
  Close(2)

  return

  110 format(//' number of type =',i5,' > ntypem =',i4,// ' Change the parameter ntypem in the code !'//)
  120 format(/' Crystal',/5x,3f16.12,3f12.5)
  130 format(//' number of atoms =',i6,' > nam =',i6,// ' Change the parameter nam in the code !'//)
  140 format(//'  Error in the indata file :')
  150 format(//' The following line is not understood :',/A,// &
          ' If it is a keyword, check the spelling.'/, &
          ' If the line is not supposed to be a keyword but contains numbers, check:'/ &
          5x,' - How many numbers must be in the line ?'/, &
          5x,' - Are there spaces between the numbers ?'/, &
          5x,' - Tabulations are forbidden !'//)
  160 format(i5,3f16.12,'  ! ',i5,3x,a2)
end

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

! Direct calculation of the inverse of a matrix 3 x 3

subroutine invermat(a,b)

  use declarations

  integer:: i, i1, i2, j, j1, j2, n 
  
  real(kind=db):: det
  real(kind=db), dimension(3,3):: a, b, cofac
  real(kind=db), dimension(2,2):: c

  n = 3

  do i = 1,n
    do j = 1,n

      do i1 = 1,n
        do j1 = 1,n
          if( i1 == i .or. j1 == j) cycle
          if( i1 < i ) then
            i2 = i1
          else
            i2 = i1 - 1
          endif
          if( j1 < j ) then
            j2 = j1
          else
            j2 = j1 - 1
          endif
          c(i2,j2) = a(i1,j1)
        end do
      end do
      cofac(i,j) = ( c(1,1)*c(2,2) - c(1,2)*c(2,1) ) * (-1)**(i+j)

    end do
  end do

  det = sum( a(1,:) * cofac(1,:) )

  b = transpose( cofac ) / det

  return
end

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

subroutine prodvec(u,v,w)

  use declarations
  real(kind=db), dimension(3):: u, v, w

  u(1) = v(2) * w(3) - v(3) * w(2)
  u(2) = v(3) * w(1) - v(1) * w(3)
  u(3) = v(1) * w(2) - v(2) * w(1)

  return
end

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

subroutine prodvec_cp(u,v,w)

  use declarations
  complex(kind=db), dimension(3):: u, v, w

  u(1) = v(2) * w(3) - v(3) * w(2)
  u(2) = v(3) * w(1) - v(1) * w(3)
  u(3) = v(1) * w(2) - v(2) * w(1)

  return
end

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

  real(kind=db) function Detmat(Mat)

  use declarations
  implicit none

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

  Detmat = Mat(1,1) * ( Mat(2,2)*Mat(3,3) - Mat(3,2)*Mat(2,3) ) - Mat(2,1) * ( Mat(1,2)*Mat(3,3) - Mat(3,2)*Mat(1,3) ) &
         + Mat(3,1) * ( Mat(1,2)*Mat(2,3) - Mat(2,2)*Mat(1,3) )

  return
end

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

! Calculation of the rotation matrix versus the Euler angles
! First rotation about z
! second rotation about y' (new y)  ( and not x' )
! third rotation about z' (new z)

subroutine Mat_Euler(Ang,Rot)

  use declarations
  implicit none

  real(kind=db):: c1, c2, c3, s1, s2, s3
  real(kind=db), dimension(3):: Ang
  real(kind=db), dimension(3,3):: Rot

  c1 = cos( Ang(1) )
  s1 = sin( Ang(1) )
  c2 = cos( Ang(2) )
  s2 = sin( Ang(2) )
  c3 = cos( Ang(3) )
  s3 = sin( Ang(3) )

! Euler notes Duc and Wikipedia (English)
  Rot(1,1) = c1*c2*c3 - s1*s3; Rot(1,2) = - c1*c2*s3 - s1*c3; Rot(1,3) = c1*s2
  Rot(2,1) = s1*c2*c3 + c1*s3; Rot(2,2) = - s1*c2*s3 + c1*c3; Rot(2,3) = s1*s2
  Rot(3,1) = - s2*c3;          Rot(3,2) = s2*s3;              Rot(3,3) = c2

  return
end

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

! Calculation of Euler angles from the rotation matrix  (Z1 Y2 Z3)

subroutine Euler_mat(icheck,Mat,alfa,beta,gamma)

  use declarations
  implicit none

  integer:: icheck

  real(kind=db):: alfa, beta, cos_b, gamma, sin_b
  real(kind=db), dimension(3,3):: Mat
  
  cos_b = Mat(3,3)
  sin_b = sqrt( abs(1 - cos_b**2) )    ! beta between 0 and pi
  beta = atan2( sin_b, cos_b )

  if( sin_b > 0.00001_db ) then
    alfa = datan2( Mat(2,3), Mat(1,3) )
    gamma = datan2( Mat(3,2), - Mat(3,1) )
  else
    gamma = 0._db
    if( abs( Mat(1,2) - 1._db ) < 0.0001_db ) then
      alfa = - pi / 2
    else
      alfa = datan2( - Mat(1,2), Mat(2,2) )
    endif
    if( cos_b < 0 ) alfa = - alfa 
  endif 

  if( icheck > 2 ) write(3,100) alfa / radian, beta / radian, gamma / radian

  return
  100 format(/' alfa, beta, gamma = ',3f10.5)
end

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

! Calculation of the Wigner matrices Dlmsms(alfa,beta,gamma) for j half integer, in the (L,m,s) basis

Subroutine Dlmsms_half_integer(Dlmsms,icheck,Inverse_rotation,Lmax,Mat_rot)

  use declarations
  implicit none

  integer:: icheck, index_J, iso, iso_p, L, Lmax, m, m_p, Two_j, Two_j_max, Two_j_p, Two_j_z, Two_j_z_p, Two_s, Two_s_p

  complex(kind=db), dimension(:,:,:), allocatable:: Dlmm_hi
  complex(kind=db), dimension(-Lmax:Lmax,2,-Lmax:Lmax,2,0:Lmax):: Dlmsms
  
  logical:: Inverse_rotation

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

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

  Two_j_max = 2 * Lmax + 1
 
  allocate( Dlmm_hi(-Two_j_max:Two_j_max,-Two_j_max:Two_j_max,(Two_j_max+1)/2) )

  call Djmm_half_integer(Dlmm_hi,icheck,Inverse_rotation,Mat_rot,Two_j_max)

  do L = 0,Lmax
    do iso = 1,2
 !     s = 1.5_db - real( iso, db )
      Two_s = 3 - 2 * iso
      do m = -L,L
        Two_j_z = 2 * m + Two_s
        if( iso == 1 .or. m == -L ) then
 ! case j+
          Two_j = 2 * L + 1
        else
          Two_j = 2 * L - 1
        endif

        do iso_p = 1,2
          Two_s_p = 3 - 2 * iso_p
          do m_p = -L,L
            Two_j_z_p = 2 * m_p + Two_s_p
            if( iso_p == 1 .or. m_p == -L ) then
 ! case j+
              Two_j_p = 2 * L + 1
            else
              Two_j_p = 2 * L - 1
            endif
          
            if( Two_j /= Two_j_p ) cycle

            index_J = ( Two_j + 1 ) / 2
            Dlmsms(m,iso,m_p,iso_p,L) = Dlmm_hi(Two_j_z,Two_j_z_p,index_j) 

          end do
        end do
      end do
    end do
  end do 

  deallocate( Dlmm_hi )
 
  if( icheck < 3 ) return
  
  do L = 0,Lmax
    write(3,110) L
    write(3,120) ( ( m_p, iso_p, m_p = -L,L ), iso_p = 1,2 )    
    do iso = 1,2
      do m = -L,L
        write(3,'(2i3,1p,100(1x,2e13.5))') m, iso, ( ( Dlmsms(m,iso,m_p,iso_p,L), m_p = -L,L ), iso_p = 1,2 )    
      end do
    end do
    
  end do
  
  
  return
  110 format(/' Wigner matrix in the L,m,s basis, L =',i3)
  120 format('  m is',100(11x,2i3,10x))
end

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

! Calculation of the Wigner matrices Djmm_hi(alfa,beta,gamma) for j half integer
! Storage with indexes from -2j to 2j (even index not used).

Subroutine Djmm_half_integer(Dlmm_hi,icheck,Inverse_rotation,Mat_rot,Two_j_max)

  use declarations
  implicit none

  integer:: icheck, index, Two_j, Two_j_max, Two_m, Two_m_p

  complex(kind=db), dimension(-Two_j_max:Two_j_max,-Two_j_max:Two_j_max,(Two_j_max+1)/2):: Dlmm_hi

  logical:: Inverse_rotation
   
  real(kind=db):: a, alfa, beta, gamma

  real(kind=db), dimension(3,3):: Mat_rot
  real(kind=db), dimension(-Two_j_max:Two_j_max,-Two_j_max:Two_j_max,(Two_j_max+1)/2):: dj

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

  call Euler_mat(icheck,Mat_rot,alfa,beta,gamma)
  
  if( Inverse_rotation ) then
    a = alfa
    alfa = pi - gamma
    gamma = - pi - a
  endif
   
  call djmmp_half_integer(beta,dj,Two_j_max)

  index = 0
  do Two_j = 1, Two_j_max, 2
    index = index + 1
  
    do Two_m = - Two_j, Two_j, 2
      do Two_m_p = - Two_j, Two_j, 2
        Dlmm_hi(Two_m,Two_m_p,index) = exp( - img * 0.5_db * ( Two_m * alfa + Two_m_p * gamma)  ) * dj(Two_m,Two_m_p,index)
      end do
    end do
  end do

  if( icheck < 4 ) return

  write(3,100) alfa*180/pi, beta*180/pi, gamma*180/pi
  
  index = 0
  do Two_j = 1, Two_j_max, 2
    index = index + 1
    write(3,'(/a16,f4.1)') ' djmm(beta), j =', index - 0.5_db
    do Two_m = - Two_j, Two_j, 2
      write(3,'(f4.1,1p,10e13.5)') 0.5*Two_m, ( dj(Two_m,Two_m_p,index), Two_m_p = - Two_j, Two_j, 2 ) 
    end do
  end do

  if( abs( alfa ) < eps10 .and. abs( gamma ) < eps10 ) return 

  index = 0
  do Two_j = 1, Two_j_max, 2
    index = index + 1
    write(3,'(/a13,f4.1)') ' Dlmm_hi, j =', index - 0.5_db
    do Two_m = - Two_j, Two_j, 2
      write(3,'(f4.1,1p,10(1x,2e13.5))') 0.5*Two_m, ( Dlmm_hi(Two_m,Two_m_p,index), Two_m_p = - Two_j, Two_j, 2 ) 
    end do
  end do

  return
  100 format(/' alfa, beta, gamma =',3f10.4)
end

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

! Calculation of djmm'(beta) for j half integer, for the Wigner matrices

Subroutine djmmp_half_integer_o(beta,dj,Two_j_max)

  use declarations
  implicit none

  integer, parameter:: Nfac = 125

  integer:: i, index, is, n, Two_j, S, S_min, S_max, Two_j_max, Two_m, Two_m_p
  
  integer, dimension(0:Two_j_max):: Fact 

  real(kind=db):: beta, cos_beta, j, m, m_p, Prefac, sin_beta, Sum
  
  real(kind=db), dimension(0:Two_j_max):: cos_b, sin_b 
  real(kind=db), dimension(-Two_j_max:Two_j_max,-Two_j_max:Two_j_max,(Two_j_max+1)/2):: dj

  dj(:,:,:) = 0._db
    
  cos_beta = cos( beta / 2 )
  sin_beta = sin( beta / 2 )
  
  cos_b(0) = 1._db
  Sin_b(0) = 1._db
  
  do i = 1,Two_j_max
    cos_b(i) = cos_b(i-1) * cos_beta  
    sin_b(i) = sin_b(i-1) * sin_beta  
  end do
   
  Fact(0) = 1

  do i = 1,Two_j_max
    Fact(i) = i * Fact(i-1)
  end do

  index = 0
  do Two_j = 1, Two_j_max, 2
    index = index + 1
  
    j = 0.5_db * Two_j
    
    do Two_m = - Two_j, Two_j, 2

      m = 0.5_db * Two_m 

      do Two_m_p = Two_m, Two_j, 2

        if( Two_m_p == Two_m .and. Two_m_p < 0 ) cycle

        m_p = 0.5_db * Two_m_p 

        S_min = nint( max( 0._db, m - m_p ) )
        S_max = nint ( min( j + m, j - m_p ) )
  
        Sum = 0._db      
        do s = S_min,S_max
        
          n = nint( m_p - m + s ) 
          is = (-1)**n
          
          Sum = Sum + is * cos_b( nint( 2*j + m - m_p - 2*s ) ) * sin_b( nint( m_p - m + + 2*s ) ) &
                    / ( Fact( nint(j+m-s) ) * Fact( s ) * Fact( nint(m_p-m+s) ) * Fact( nint(j-m_p-s) ) )    
        end do
 
        Prefac = sqrt( real( Fact( nint(j+m_p) ) * Fact( nint(j-m_p) ) * Fact( nint(j+m) ) * Fact( nint(j-m) ),db ) )
                
        dj(Two_m_p,Two_m,index) = Prefac * Sum
        if( Two_m_p == Two_m ) then
          dj(-Two_m_p,-Two_m,index) = dj(Two_m_p,Two_m,index)
        else
          dj(Two_m,Two_m_p,index) = (-1)**nint(m_p-m) * dj(Two_m_p,Two_m,index)
        endif
            
      end do
          
    end do 
  
  end do

  return
end

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

Subroutine djmmp_half_integer(beta,dj,Two_j_max)

  use declarations
  implicit none

  integer, parameter:: Nfac = 125

  integer:: i, index, is, n, Two_j, S, S_min, S_max, Two_j_max, Two_m, Two_m_p
  
  integer, dimension(0:Two_j_max):: Fact 

  real(kind=db):: beta, cos_beta, Fac, j, m, m_p, Prefac, sin_beta, Sum
  
  real(kind=db), dimension(0:Two_j_max):: cos_b, sin_b, sqrt_Fact 
  real(kind=db), dimension(-Two_j_max:Two_j_max,-Two_j_max:Two_j_max,(Two_j_max+1)/2):: dj

  dj(:,:,:) = 0._db
    
  cos_beta = cos( beta / 2 )
  sin_beta = sin( beta / 2 )
  
  cos_b(0) = 1._db
  Sin_b(0) = 1._db
  
  do i = 1,Two_j_max
    cos_b(i) = cos_b(i-1) * cos_beta  
    sin_b(i) = sin_b(i-1) * sin_beta  
  end do
   
  Fact(0) = 1

  do i = 1,Two_j_max
    Fact(i) = i * Fact(i-1)
  end do
  do i = 0,Two_j_max
    sqrt_Fact(i) = sqrt( real( Fact(i), db ) )
  end do

  index = 0
  do Two_j = 1, Two_j_max, 2
    index = index + 1
  
    j = 0.5_db * Two_j
    
    do Two_m = - Two_j, Two_j, 2

      m = 0.5_db * Two_m 

      do Two_m_p = Two_m, Two_j, 2

        if( Two_m_p == Two_m .and. Two_m_p < 0 ) cycle

        m_p = 0.5_db * Two_m_p 

        S_min = nint( max( 0._db, m - m_p ) )
        S_max = nint ( min( j + m, j - m_p ) )
  
        Sum = 0._db      
        do s = S_min,S_max
        
          n = nint( m_p - m + s ) 
          is = (-1)**n
          
          Fac = is * cos_b( nint( 2*j + m - m_p - 2*s ) ) * sin_b( nint( m_p - m + 2*s ) ) / Fact( nint(j+m-s) )
          Fac = Fac / Fact( s )
          Fac = Fac / Fact( nint(m_p-m+s) )
          Fac = Fac / Fact( nint(j-m_p-s) )
            
          Sum = Sum + Fac    
        end do
 
        Prefac = sqrt_Fact( nint(j+m_p) ) * sqrt_Fact( nint(j-m_p) ) * sqrt_Fact( nint(j+m) ) * sqrt_Fact( nint(j-m) )
                
        dj(Two_m_p,Two_m,index) = Prefac * Sum
        if( Two_m_p == Two_m ) then
          dj(-Two_m_p,-Two_m,index) = dj(Two_m_p,Two_m,index)
        else
          dj(Two_m,Two_m_p,index) = (-1)**nint(m_p-m) * dj(Two_m_p,Two_m,index)
        endif
            
      end do
          
    end do 
  
  end do

  return
end





