subroutine egrad_siesta(e,geomcart,gradcart)

  ! This routines calculates the energy and the gradient of 
  ! a molecule, using 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 Path_module, only : Nat, renum,Order,OrderInv, Coord, ProgExe,v_direct,latr
  use Io_module

  IMPLICIT NONE

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


    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


  ! Energy 
  REAL(KREAL), INTENT (OUT) :: e
  ! Geometry in cartesian coordinates
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
  ! Gradient calculated at Geom geometry
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)

  ! ======================================================================

  character(LCHARS) :: LINE

  logical           :: debug, TChk
  LOGICAL, SAVE :: first=.true.

  REAL(KREAL), ALLOCATABLE :: X(:),Y(:),Z(:) ! Nat
  REAL(KREAL), ALLOCATABLE :: GeomCartLoc(:,:) !Nat,3

  INTEGER(KINT) :: iat, i, J,Idx

  CHARACTER(LCHARS) :: FileIn, FileOut, FileForces,FileE

  CHARACTER(3) :: SepIn=' < ', SepOut=' > '
  CHARACTER(VLCHARS), SAVE :: RunCommand

  ! ======================================================================

  debug=valid('EGRAD').OR.valid('egrad_siesta')

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

  if (first) THEN
     First=.FALSE.
     call system("uname > Path.tmp")
     OPEN(IOTMP,FILE='Path.tmp')
     READ(IOTMP,'(A)') Line
     CLOSE(IOTMP)
     IF (index(Line,'CYGWIN').NE.0) THEN
        SepIn='   '
        SepOut='   '
     END IF
  END IF

  ALLOCATE(GeomCartLoc(Nat,3))

  gradcart=0.
  E=0.
  FileIn=Trim(CalcName) // "." // Trim(ISuffix)
  FileOut=Trim(CalcName) // "." // Trim(OSuffix)
  FileForces=Trim(Siesta_Label) // ".FA"
  FileE='BASIS_ENTHALPY'

  RunCommand=Trim(ProgExe) // SepIn // Trim(FileIn) // SepOut // Trim(FileOut)

  IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(RunCommand)

  ! we create the input file

  OPEN(IOTMP,File=FileIn)

  Call WriteList(Siesta_input,Unit=IOTMP)

  WRITE(IOTMP,*) 

  IF (V_Direct=='DIRECT') THEN
     ALLOCATE(X(Nat),Y(Nat),Z(Nat))
     X=0.
     Y=0.
     Z=0.
     DO I=1,Nat
        DO J=1,3
           X(I)=X(I)+GeomCart(I,J)*Latr(1,J)
           Y(I)=Y(I)+GeomCart(I,J)*Latr(2,J)
           Z(I)=Z(I)+GeomCart(I,J)*Latr(3,J)
        END DO
     END DO
     GeomCartLoc(1:Nat,1)=X(1:Nat)
     GeomCartLoc(1:Nat,2)=Y(1:Nat)
     GeomCartLoc(1:Nat,3)=Z(1:Nat)
     DEALLOCATE(X,Y,Z)
  ELSE
     GeomCartLoc=GeomCart
  END IF

  WRITE(IOTMP,'(1X,A)')  '%block AtomicCoordinatesAndAtomicSpecies' 

  DO I=1,Nat
     If (renum) THEN
        Iat=Order(I)
        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(Iat,:)/Siesta_Unit_Read, IdxSpecies(Iat),TRIM(Siesta_Paste(I))
     ELSE
        Iat=OrderInv(I)
        WRITE(IOTMP,'(1X,3(1X,F15.8),1X,I5,1X,A)') GeomCart(I,:)/Siesta_Unit_Read, IdxSpecies(Iat), TRIM(Siesta_Paste(Iat))
     END IF
  END DO

  WRITE(IOTMP,'(1X,A)')  '%endblock AtomicCoordinatesAndAtomicSpecies'
  WRITE(IOTMP,*) 

  CLOSE(IOTMP)

  IF (debug) THEN
     WRITE(*,*) 'DBG EGRAD ' //  Trim(FileIn),' COORD=',TRIM(COORD)
     call system('cat ' // Trim(FileIn))
  END IF

  call system(RunCommand)

  if (debug) WRITE(*,*) 'DBG EGRAD, back from calculation'

  ! Whatever the coordinate system, we read the cartesian gradient
  ! and we convert it to Zmat or other if needed

! in Siesta, the forces are in the file 'Label'.FA so this is easy to read :)
  OPEN(IOTMP,FILE=FileForces, STATUS='old')
! firt line is the atom number
  READ(IOTMP,'(A)') LINE

  if (debug) WRITE(*,*) "GradCart reading - renum=",Renum
  DO I=1,Nat
     Iat=I
     IF (renum) Iat=Order(I)
     READ(IOTMP,'(6X,3(F12.6))')  GradCart(3*Iat-2:3*Iat)
     if (debug) WRITE(*,*) "I,Iat,GradCart:",I,Iat,GradCart(3*Iat-2:3*Iat)
  END DO

! Siesta gives the Forces in eV/Ang, so we convert it into the 
! gradient in au/Angstrom 
  gradCart=-1.0d0*ev2au*gradCart

  CLOSE(IOTMP)

! We now have to read the energy

  INQUIRE(File=FileE,EXIST=TChk)
  If (Tchk) THEN
! We have a BASIS_ENTHALPY file. Great. 
! Energy is the last number on the second line
     OPEN(IOTMP,File=FileE,Status='Old')
     READ(IOTMP,*) Line
     READ(IOTMP,'(A)') Line
     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)
     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)
     Line=Line(Idx+1:)
     if (debug) WRITE(*,*) "Line E:",TRIM(LINE)
     Read(Line,*) E
  ELSE
! We do not have this file... we have to look for it in the .out file
     FileE=FileOut
     OPEN(IOTMP,File=FileE,Status='Old')
! We search for  "siesta:         Total ="
     Line=""
     DO WHILE (INDEX(Line,"siesta:         Total =")/=1)
        READ(IOTMP,*) Line
     END DO

     Idx=INDEX(TRIM(LINE),' ',BACK=.TRUE.)
     Line=Line(Idx+1:)
     Read(Line,*) E
  END IF
  CLOSE (IOTMP)

  WRITE(*,*) 'E=',E
  WRITE(*,*) 'Gradcart='
  DO I=1,Nat
     WRITE(*,'(1X,3(1X,F15.8))') GradCart(3*I-2:3*I)
  END DO

  DeAllocate(GeomCartLoc)

!  Call Die('Egrad_Siesta','Ah ah.',UNIT=IOOUT)
  if (debug) Call header("Egrad_siesta Over")



  RETURN

  ! ======================================================================
end subroutine egrad_siesta
