subroutine egrad_vasp(e,geomcart,gradcart)

  ! This routines calculates the energy and the gradient of 
  ! a molecule, using VASP


  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, dzdc, &
       indzmat,Nom,Atome, massat, unit,Nfroz,Frozen,ProgExe,Latr,FrozAtoms
  use Io_module

  !
  IMPLICIT NONE

  ! Energy (calculated if F300K=.F., else estimated)
  REAL(KREAL), INTENT (OUT) :: e
  ! Nb degree of freedom
  REAL(KREAL), INTENT (IN) :: geomcart(Nat,3)
  ! Gradient calculated at Geom geometry
  REAL(KREAL), INTENT (OUT) :: gradcart(3*NAt)

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

  character(LCHARS) :: LINE

  REAL(KREAL) :: d, a_val, Pi

  INTEGER(KINT) :: iat, jat, kat, i, j, n3at, Nfr
  REAL(KREAL) :: RTmp1, RTmp2, RTmp3

  !
  CHARACTER(132), PARAMETER :: FileIn='POSCAR', FileOut='OUTCAR'

  CHARACTER(VLCHARS), SAVE :: RunCommand
  INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit

  INTEGER(KINT) :: IGeom, Idx
!  REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:)
  REAL(KREAL) :: MRot(3,3), Rmsd

  REAL(KREAL) :: X,Y,Z
  LOGICAL :: Debug

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

  Pi=dacos(-1.0d0)
  debug=valid('EGRAD')
  if (debug) Call header('Entering Egrad_VASP')

  n3at=3*nat

!  ALLOCATE(GeomTmpC(n3at))


  gradcart=0.
  E=0.
  RunCommand=Trim(ProgExe) 

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

  ! we create the input file

!   write(*,*) "Coucou 2"

  OPEN(IOTMP,File=FileIn)

  ! we convert the coordinates into Cartesian coordinates
!     GeomTmpC=Reshape(GeomCart,SHAPE=(/3*Nat/))

   Call PrintGeomVasp(Vasp_Title,geomcart,Latr,FrozAtoms,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_VASP, back from calculation'
!  if (debug) WRITE(*,*) "DBG EGRAD_VASP Order(I)",Order(1:Nat)

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

  OPEN(IOTMP,FILE=FileOut, STATUS='old')
  ! We first search for the forces
  READ(IOTMP,'(A)',END=999,ERR=999) LINE
! We search for the  last part of the OUTCAR file, after wavefunction is converged
  DO WHILE (INDEX(LINE,'EDIFF is reached')==0)
     READ(IOTMP,'(A)',END=999,ERR=999) LINE
  END DO
  READ(IOTMP,'(A)',END=999,ERR=999) LINE
  DO WHILE (INDEX(LINE,'TOTEN')==0)
     READ(IOTMP,'(A)',END=999,ERR=999) LINE
  END DO
  Line=Line(26:)
  READ(LINE,*) E
  if (debug) WRITE(*,'(1X,A,F20.6)') "DBG Egrad_VASP E=",E

! We search for the forces
  DO WHILE (INDEX(LINE,'TOTAL-FORCE')==0)
     READ(IOTMP,'(A)',END=999,ERR=999) LINE
  END DO
  READ(IOTMP,'(A)',END=999,ERR=999) LINE
  DO I=1,Nat
     Iat=I
     IF (renum) Iat=Order(I)
     READ(IOTMP,*,END=999,ERR=999)  Rtmp1,RTmp2,RTmp3,GradCart(3*Iat-2:3*Iat)
  END DO


  if (debug) THEN
     WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - original order"
     DO I=1,Nat
        Iat=Order(I)
        WRITE(*,'(1X,I5,3(1X,F20.6))') I,GradCart(3*Iat-2:3*Iat)
     END DO
     WRITE(*,*) "DBG EGRAD_VASP - GradCart=FORCES - internal order"
     DO I=1,Nat
        Iat=Order(I)
        WRITE(*,'(3(1X,I5),3(1X,F20.6))') I,Iat,OrderInv(I),GradCart(3*I-2:3*I)
     END DO
  END IF

  ! VASP gives the Forces in eV/Angstrom, so we convert it into the 
  ! gradient in ua/Angstrom 
    gradCart=-1.0d0*ev2Au*gradCart

  !  if (debug.AND.(COORD.EQ.'ZMAT')) THEN
  !     DO WHILE (INDEX(LINE,'Internal Coordinate Forces')==0)
  !        READ(IOTMP,'(A)') LINE
  !!        WRITE(*,*) 'Pas bon:' // TRIM(LINE)
  !     END DO
  !     WRITE(*,*) TRIM(LINE)
  !     DO I=1,Nat+2
  !        READ(IOTMP,'(A)') LINE
  !        WRITE(*,*) TRIM(LINE)
  !     END DO
  !  END IF

  CLOSE(IOTMP)

!  DeALLOCATE(GeomTmpC)


  if (debug) Call header('Egrad_VASP Over')

  RETURN

999 CONTINUE
  WRITE(*,*) 'We should not be here !!!!'
  STOP
  ! ======================================================================
end subroutine egrad_vasp




