 subroutine egrad_ext(e,geomcart,gradcart)

  ! This routines calculates the energy and the gradient of 
  ! a molecule, using an external code


  use Path_module, only : Nat, renum,Order,OrderInv,AtName, Coord, dzdc, indzmat,Nom,Atome, massat, unit,ProgExe
  use Io_module

  !
  IMPLICIT NONE

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

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

  character(LCHARS) :: LINE

  logical           :: debug
  LOGICAL :: fexist, FSim
  LOGICAL, SAVE :: first=.true.
  LOGICAL, ALLOCATABLE :: Done(:)

  REAL(KREAL), SAVE :: Eold=1.e6
  REAL(KREAL) :: d, a_val, Pi

  REAL(KREAL) :: coef,x
  INTEGER(KINT) :: iat, jat, kat, i, j, n3at, absidg, absidg2
  INTEGER(KINT) :: kTmp, Istart,ITmp1,ITmp2, Idx
  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11

  !
  CHARACTER(132) :: FileIn,FileOut

  CHARACTER(LCHARS) :: ListName, TitleTmp, CH32SVAR1
  CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand
  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
  LOGICAL   :: FRdyn,FStopNucl, FTmp, Tchk,Tchk1
  INTEGER(kint) :: NbLists,LastIt,Nat4,NStep, Firstit
  INTEGER(KINT) :: NStepAd,NStepTmp
  INTEGER(KINT) :: IShowTmp
  REAL(KREAL)   :: FricPsiMin,FricPsiT, FricNuclTmp


  INTEGER(KINT), PARAMETER :: NbExtName=4

  INTEGER(KINT) :: ICouc

  INTEGER(KINT) :: ILine
  INTEGER(KINT) :: NPulses, PulseLen, NStepWarm,NStepAv,NStepEq,NStepTot
  REAL(KREAL) :: TempWarm,Mean,Slope,Dev,RTmp1,RTmp2

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

  LOGICAL, EXTERNAL :: valid

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


  Pi=dacos(-1.0d0)
  n3at=3*nat

  debug=valid('EGRAD')
  if (debug) WRITE(*,*) '================ Entering Egrad_ext ===================='

  RunCommand=Trim(Adjustl(ProgExe))
  FileIn=Trim(CalcName) // Trim(ISuffix)
  FileOut=Trim(CalcName) // Trim(OSuffix)

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

  ! we create the input file

  OPEN(IOTMP,File=FileIn)

  WRITE(IOTMP,'(1X,I10)') NAt
  WRITE(IOTMP,'(1X,A)') Coord

     DO I=1,Nat
        If (renum) THEN
           Iat=Order(I)
           WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),GeomCart(Iat,:)
        ELSE
           Iat=OrderInv(I)
           WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomCart(I,:)
        END IF
     END DO

  call system(RunCommand)

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

  OPEN(IOTMP,FILE=FileOut, STATUS='old')
  ! We first search for the forces
  READ(IOTMP,*) e
  DO I=1,Nat
     Iat=I
     IF (renum) Iat=Order(I)
     READ(IOTMP,*)  GradCart(3*Iat-2:3*Iat)
  END DO

  CLOSE(IOTMP)


  if (debug) WRITE(*,*) '================  Egrad_ext Over ===================='

  RETURN

999 CONTINUE
  if (.NOT.Ftmp) WRITE(*,*) 'We should not be here !!!!'
  STOP
  ! ======================================================================
end subroutine egrad_ext
