 SUBROUTINE ReadInput_siesta

! This routine reads an input template for Siesta

!----------------------------------------------------------------------
!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 
!  Centre National de la Recherche Scientifique,
!  Université Claude Bernard Lyon 1. All rights reserved.
!
!  This work is registered with the Agency for the Protection of Programs 
!  as IDDN.FR.001.100009.000.S.P.2014.000.30625
!
!  Authors: P. Fleurat-Lessard, P. Dayal
!  Contact: optnpath@gmail.com
!
! This file is part of "Opt'n Path".
!
!  "Opt'n Path" is free software: you can redistribute it and/or modify
!  it under the terms of the GNU Affero General Public License as
!  published by the Free Software Foundation, either version 3 of the License,
!  or (at your option) any later version.
!
!  "Opt'n Path" is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU Affero General Public License for more details.
!
!  You should have received a copy of the GNU Affero General Public License
!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.
!
! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,
! for commercial licensing opportunities.
!----------------------------------------------------------------------

  use VarTypes
  use Path_module
  use Io_module

  IMPLICIT NONE

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
  INTERFACE
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

    FUNCTION SearchInput(Input,String,Line,Clean)  Result (Found)

      Use Vartypes
      use io_module
! Input
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
      CHARACTER(*), INTENT(IN) :: String
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean

! Output
      TYPE (Input_line), POINTER, INTENT(OUT) :: Line

      LOGICAL :: Found

    END FUNCTION SearchInput

    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)

      Use VarTypes

      implicit none
! Input
      CHARACTER(*), INTENT(IN) :: Line
      CHARACTER(*), INTENT(IN) :: String
      LOGICAL, OPTIONAL, INTENT(IN) :: Case
      LOGICAL, OPTIONAL, INTENT(IN) :: Back
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean

! Output
! the position of String in Line (the first one) unless Back is present
      INTEGER(KINT) :: Pos
    END FUNCTION InString

    SUBROUTINE die(routine, msg, file, line, unit)

      Use VarTypes
      Use io_module

      implicit none
!--------------------------------------------------------------- Input Variables
      character(len=*), intent(in)           :: routine, msg
      character(len=*), intent(in), optional :: file
      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die

    SUBROUTINE Warning(routine, msg, file, line, unit)

      Use VarTypes
      Use io_module

      implicit none

      character(len=*), intent(in)           :: routine, msg
      character(len=*), intent(in), optional :: file
      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE Warning


 SUBROUTINE WriteList(Input,Unit)

! This routine reads an input template for Siesta

  use VarTypes
  use Io_module

  IMPLICIT NONE

! Input
      TYPE (Input_line), POINTER, INTENT(IN) :: Input
      INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit

    END SUBROUTINE WriteList

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
 END  INTERFACE
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  CHARACTER(132) ::  Line,Line2
  INTEGER(KINT) ::  Idx,LineL
  INTEGER(KINT) ::  IAt,I
  INTEGER(KINT) :: IoRead, ITmp
  REAL(KREAL) :: B(3),Xtmp, Ytmp, Ztmp
  REAL(KREAL) :: Lata,Latb,Latc,Alpha,Beta,Gamma
  REAL(KREAL), ALLOCATABLE :: GeomCartTmp(:,:) ! (Nat,3)
  REAL(KREAL) :: LatLoc(3,3)
  
  LOGICAL :: Debug
  LOGICAL :: FSpecies, FCoord

  TYPE(Input_Line), POINTER :: Search,Bla

  Debug=Valid("readinput").OR.Valid("readinput_siesta")

 if (debug) Call Header("Entering ReadInput_Siesta")

  ! We read the Siesta input file

     ALLOCATE(Siesta_Input)
     NULLIFY(Siesta_Input%next)
     NULLIFY(Siesta_Input%prev)
     Current => Siesta_Input

     FSpecies=.FALSE.
     FCoord=.FALSE.

     READ(IOIN,'(A)',iostat=IoRead) Line

     DO WHILE (IoRead==0)
        Line=AdjustL(Line)
        if (debug) WRITE(*,*) 'Line:', Line
        current%Line=TRIM(Line)
        ALLOCATE(current%next)
        Current%next%prev => Current
        Current => Current%next
        current%line="toto"
        Nullify(Current%next)

        READ(IOIN,'(A)',iostat=IoRead) Line
     END DO
! With this procedure, the last current is not valid and should be deleted.
       Bla => Current%prev
       nullify(Bla%next)
       deallocate(Current)
       current => Bla

     if (debug) THEN
        WRITE(*,*) 'Input just read'
        Call WriteList(Siesta_input,Unit=IOOUT)
     END IF

! We analyse the input

! We look for the NumberofAtoms
     If (SearchInput(Siesta_Input,"NUMBEROFATOMS",Search,Clean=".-_")) THEN
        Line=AdjustL(Search%line)
        Idx=Index(Line," ")
        Line2=Trim(AdjustL(Line(Idx+1:)))
        Read(Line2,*) IAt
        if (Iat/=Nat) THEN
           Call Die('ReadInput_siesta','Nat in FDF sample different from  Nat Path input', Unit=IOOUT)
        END IF
     ELSE
! There is no atom number defined !!!
           Call Warning('ReadInput_siesta','No NumberOfAtoms in FDF sample', Unit=IOOUT)
           WRITE(current%Line,'(1X,"NumberOfAtoms ",I5)') Nat
           ALLOCATE(current%next)
           Current%next%prev => Current
           Current => Current%next
           Nullify(Current%next)
         
        END IF

! We look for the NumberofSpecies
        IF (SearchInput(Siesta_Input,"NUMBEROFSPECIES",Search,Clean=".-_")) THEN
           Line=AdjustL(Search%line)
           Idx=Index(Line," ")
           Line2=Trim(AdjustL(Line(Idx+1:)))
           Read(Line2,*) Siesta_NbSpecies
           ALLOCATE(Siesta_SpeciesMass(Siesta_NbSpecies))
           ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies))
        END IF

! We look for SystemLabel
        If (SearchInput(Siesta_Input,"SYSTEMLABEL",Search,Clean=".-_")) THEN
           Line=AdjustL(Search%line)
           Idx=Index(Line," ")
           Siesta_Label=Trim(adjustl(Line(Idx+1:)))
        ELSE
           Siesta_label='siesta'
        END IF

! We look for the ChemicalSpeciesLabel block
        IF (SearchInput(Siesta_Input,"CHEMICALSPECIESLABEL",Search,Clean=".-_")) THEN
           if (.NOT.ALLOCATED(IdxSpecies)) ALLOCATE(IdxSpecies(NAt))
           I=0
           Search => Search%next
           DO WHILE (InString(Search%line,'ENDBLOCK',Case=.FALSE.,Clean=".-_")==0)
              I=I+1
              if (I>Siesta_NbSpecies) THEN
                 Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel','Found more line in this block than NbSpecies !')
              END IF
              Line=AdjustL(Search%Line)
              Read(Line,*) ITmp
              Idx=Index(Line,' ')
              Line=AdjustL(Line(Idx+1:))
              Read(Line,*) Idx
              Siesta_SpeciesMass(ITmp)=Idx
              Idx=Index(Line,' ')
              Line=AdjustL(Line(Idx+1:))
              Idx=Index(Line,' ')
              Siesta_SpeciesName(ITmp)=Line(1:Idx-1)
              Search => Search%next
           END DO
           IF (I/=Siesta_NbSpecies) Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel', &
                'Number of lines in this block different from NbSpecies')
        ELSE
           Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel', &
                'Block ChemicalSpeciesLabel is mandatory !')
        END IF


! We look for the AtomicCoordinatesAndAtomicSpecies  block
        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESANDATOMICSPECIES",Search,Clean=".-_")) THEN
           ALLOCATE(Siesta_Paste(Nat))
           Current=>Search
           DO I=1,NAt
              Search => Search%next
              Read(Search%line,*) XTmp,YTmp,ZTmp,IdxSpecies(I)
! We give a name to this atom
              IF (AtName(I)/="") THEN
!                    Write(Line,'(A,I5,1X,A2,1X,A32)') 'AtName & SpeciesName for atom ',I,AtName(I), &
!                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))
!                    WRITE(*,'(1X,A)') Line
                 If (InString(Atname(I),TRIM(Siesta_SpeciesName(IdxSpecies(I))))==0) THEN
                    Write(Line,'(A,I5,1X,A2,1X,A32)') 'AtName /= SpeciesName for atom ',I,AtName(I), &
                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))
                    Call Die('Readinput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
                 END IF
! We look for the Mass number also
                 IF (Atome(I)/=Siesta_SpeciesMass(IdxSpecies(I))) THEN
                    Write(Line,'(A,I5,1X,I5,1X,I5)') 'AtMass /= SpeciesMass for atom ',I,Atome(I), &
                         Siesta_SpeciesMass(IdxSpecies(I))
                    Call Die('ReadInput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
                 END IF
              ELSE
                 AtName(I)=AdjustL(Trim(Siesta_SpeciesName(IdxSpecies(I))))
                 Atome(I)=Siesta_SpeciesMass(IdxSpecies(I))
              END IF
! we look for something else at the end of the line
! We save everything but the x,y,z, and species description              
              Line=AdjustL(search%line)
! We skip x
              Idx=Index(Line,' ')
              Line=AdjustL(Line(Idx+1:))
! we skip y
              Idx=Index(Line,' ')
              Line=AdjustL(Line(Idx+1:))
! we skip z
              Idx=Index(Line,' ')
              Line=AdjustL(Line(Idx+1:))
! we skip species
              Idx=Index(Line,' ')
              Line=AdjustL(Line(Idx+1:))
              Siesta_Paste(I)=TRIM(Line)
           END DO
           if (debug) THEN
              WRITE(*,*) 'Input before deletion of coord block'
              Call WriteList(Siesta_input)
           END IF
! We will now delete this block from our input sample as it will then be
! written directly by Opt'n Path
! Search%next point on %endblock
           search => search%next
           IF (ASSOCIATED(Search%next)) THEN
! we are not at the end of the input file
              if (ASSOCIATED(Current%Prev)) THEN
                 Search%next%prev => current%prev
                 current%prev%next => search%next
              ELSE
! the coordinate block is at the begining of the input
                 siesta_input => Search%next
                 Nullify(Siesta_Input%Prev)
              END IF
           ELSE
! the coordinate block is the last part of the input
              nullify(current%prev%next)
           END IF
           if (debug) THEN
              WRITE(*,*) 'Input after deletion of coord block 1'
              Call WriteList(Siesta_input)
           END IF

           DO I=1,Nat+2
              bla=>current
              current=> current%next
              deallocate(bla)
           END DO
           if (debug) THEN
              WRITE(*,*) 'Input after deletion of coord block 2'
              Call WriteList(Siesta_input)
           END IF
        ELSE
           IF (SearchInput(Siesta_Input,"ZMATRIX",Search,Clean=".-_")) THEN
              call Die('ReadInput_Siesta','For now, I need the full block'//&
                   'AtomicCoordinatesAndAtomicSpecies, ZMatrix not yet handled.')
           ELSE
              call Die('ReadInput_Siesta','For now, I need the full block' //&
              'AtomicCoordinatesAndAtomicSpecies.')
           END IF
        END IF

        IF (SearchInput(Siesta_Input,"LATTICE",Search,Clean=".-_")) THEN
! This is a pbc calculation
           FPBC=.TRUE.
           IPer=3
           Kabeg=-1
           kaEnd=1
           kbBeg=-1
           kbEnd=1
           kcBeg=-1
           kcEnd=1
         

         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN
           Line=Adjustl(Search%Line)
! We discar the label
           Idx=Index(Line,' ')
           Line=AdjustL(Line(Idx+1:))
! we read the value
           Read(Line,*) Siesta_LatticeConstant
! We discard the value
           Idx=Index(Line,' ')
           Line=AdjustL(Line(Idx+1:))
           LineL=Len_Trim(Line)
! We read the unit
           If (LineL>0) THEN
              Call UpCase(Line)
              SELECT CASE (Line) 
                CASE ('ANG')
                   Siesta_Lat_Unit=1.d0
                CASE ('BOHR')
                   Siesta_Lat_Unit=a0
              END SELECT
           ELSE
              Siesta_Lat_Unit=1.d0
           END IF
         ELSE
! We have no LatticeConstant but we need one for our calculation
            Siesta_LatticeConstant=1.
            Siesta_Lat_Unit=1.d0
         END IF

         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN
           search => search%next
           Line=Adjustl(Search%Line)
           Read(Line,*) Lat_a(1:3)
           search => search%next
           Line=Adjustl(Search%Line)
           Read(Line,*) Lat_b(1:3)
           search => search%next
           Line=Adjustl(Search%Line)
           Read(Line,*) Lat_c(1:3)

         ELSEIF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN
           search => search%next
           Line=Adjustl(Search%Line)
           Read(Line,*) Lata,Latb,Latc,Alpha,Beta,Gamma

            Lat_a=0.
            Lat_b=0.
            Lat_c=0.
            Lat_a(1)=Lata
            Lat_b(1)=LatB*cos(Gamma*Pi/180d0)
            Lat_b(2)=LatB*sin(Gamma*Pi/180d0)
            Lat_c(1)=Latc*cos(Beta*Pi/180.)
            Lat_c(2)=(Latb*Latc*cos(Alpha*Pi/180.)-Lat_b(1)*Lat_c(1))/Lat_b(2)
            Lat_c(3)=sqrt(Latc*Latc-Lat_c(1)**2-Lat_c(2)**2)
         ELSE
! We have a lattice constant but nothing else. We issue a warning just in case,
! and we take the defaut values 1. 1. 1. 90. 90. 90.
            Line="LatticeConstant found, but not LatticeVectors nor LatticeParameters. Taking 1. 1. 1. 90. 90. 90."
            Call Warning("ReadInput_siesta",Line)
            Lat_a=(/1.,0.,0./)
            Lat_b=(/0.,1.,0./)
            Lat_c=(/0.,0.,1./)
         END IF

         IF (debug) THEN 
            WRITE(*,*) "Lattice vectors before scaling with LatticeConstant"
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c
         END IF
           Lat_a=Lat_a*Siesta_LatticeConstant*Siesta_Lat_Unit
           Lat_b=Lat_b*Siesta_LatticeConstant*Siesta_Lat_Unit
           Lat_c=Lat_c*Siesta_LatticeConstant*Siesta_Lat_Unit

         IF (debug) THEN 
            WRITE(*,*) "Lattice vectors After scaling with LatticeConstant"
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b
            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c
         END IF

      END IF

         IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN
           Line=Adjustl(Search%Line)
           Idx=Index(Line,' ')
           Line=AdjustL(Line(Idx+1:))
           Call UpCase(Line)
           SELECT CASE (Line) 
              CASE ('ANG','NOTSCALEDCARTESIANANG')
                 Siesta_Unit_Read=1.d0
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
                 IF (FPBC) THEN
                    Siesta_Unit_Read=1.d0
                    V_DIRECT="DIRECT"
                    Latr(1:3,1)=Lat_a
                    Latr(1:3,2)=Lat_b
                    Latr(1:3,3)=Lat_c
                    B=1.
! TO DO: replace by LAPACK
                    CALL Gaussj(Latr,3,3,B,1,1)
                    LatLoc(1:3,1)=Lat_a
                    LatLoc(1:3,2)=Lat_b
                    LatLoc(1:3,3)=Lat_c
                    Allocate(GeomCartTmp(Nat,3))
                    Do I=1,NGeomI
                       GeomCartTmp=Reshape(XyzGeomI(I,1:3,1:Nat),(/Nat,3/),ORDER=(/2,1/))
! TO DO: shall we replace MatMull by DGEMM ?
                       XyzGeomI(I,1:3,1:Nat)=Reshape(MatMul(GeomCartTmp,Latloc),(/3,Nat/),ORDER=(/2,1/))
                   END DO
                 ELSE
                    WRITE(Line2,*) "AtomicCoordinatesFormat=" // Trim(Line) &
               //", but there is no information about Lattice parameters"
                    Call Die("ReadInput_siesta",Line2)
                 END IF
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
                 Siesta_Unit_Read=a0
                 IF (INPUT=='SIESTA') THEN
! We have read the coordinates, but not the unit. This is corrected here
                    XyZGeomI=XyzGeomI*a0
                 END IF
              CASE ('SCALEDCARTESIAN')
                 Siesta_Unit_Read=Siesta_LatticeConstant
                 If (INPUT=='SIESTA') THEN
                    XyzGeomI=XyzGeomI*Siesta_Unit_Read
                 END IF
            END SELECT
         ELSE
            Siesta_Unit_Read=1.d0
         END IF

        IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN
           Line=Adjustl(Search%Line)
           Idx=Index(Line,' ')
           Line=AdjustL(Line(Idx+1:))
           Call UpCase(Line)
           SELECT CASE (Line) 
              CASE ('ANG','NOTSCALEDCARTESIANANG')
                 Siesta_Unit_Write=1.d0
              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')
                 Siesta_Unit_Write=1.d0
                 V_DIRECT_Write='DIRECT'
              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')
                 Siesta_Unit_Write=Unit ! 1/a0
            END SELECT
         END IF
         IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN
            Line=Adjustl(search%line)
            Call Upcase(Line)
            Call CleanString(Line,".-_")
            If (Index(Line,"%BLOCK")/=0) THEN
               Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")
            END IF
         END IF

!         Call Die('Readinput_siesta:fin',' Lecture finie')

     if (debug) Call Header("Exiting ReadInput_Siesta")

   END SUBROUTINE READINPUT_SIESTA
