subroutine egrad_vasp(e,geomcart,gradcart)

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

!----------------------------------------------------------------------
!  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,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) :: Pi

  INTEGER(KINT) :: iat, i, n3at
  REAL(KREAL) :: RTmp1, RTmp2, RTmp3

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

  CHARACTER(VLCHARS), SAVE :: RunCommand

!  REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:)

  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




