subroutine Egrad(E,Geom,Grad,NCoord,IGeom,IOpt,GeomCart)

  ! This routines calculates the energy E and the gradient Grad of 
  ! a molecule with Geometry Geom (may be in internal coordinates),
  ! using for now, either gaussian or Ext, more general later.

  use Path_module, only : Nat,AtName,Coord,dzdc,indzmat,Nom,Atome,massat,unit, &
                          prog,NCart,XyzGeomF,IntCoordF,XyzGeomI, renum, OrderInv
              ! IntCoordF(NGeomF,NCoord),GeomOld_all(NGeomF,NCoord)
              ! allocated in Path.f90
  use Io_module
  IMPLICIT NONE

  ! Energy (calculated if F300K=.F., else estimated)
  REAL(KREAL), INTENT (OUT) :: E
  ! NCoord: Number of the degrees of freedom
  ! IGeom: index of the geometry.
  INTEGER(KINT), INTENT (IN) :: NCoord, IGeom, IOpt
  ! Geometry at which gradient is calculated (cf Factual also):
  REAL(KREAL), INTENT (INOUT) :: Geom(NCoord)
  ! Gradient calculated at Geom geometry:
  REAL(KREAL), INTENT (OUT) :: Grad(NCoord)
  ! Cartesian geometry corresponding to (Internal Geometry) Geom:
  REAL(KREAL), INTENT (OUT) :: GeomCart(Nat,3)

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

  CHARACTER(LCHARS) :: LINE
  LOGICAL :: debug
  LOGICAL, SAVE :: first=.true.
  REAL(KREAL), SAVE :: Eold=1.e6
  REAL(KREAL), ALLOCATABLE :: GradTmp(:)
  REAL(KREAL), ALLOCATABLE :: GradCart(:)
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
  REAL(KREAL) :: d, a_val, Pi

  INTEGER(KINT) :: iat,jat,kat,i,j,IBeg,Idx
  INTEGER(KINT), PARAMETER :: IOLOG=12, IOCOM=11

  CHARACTER(LCHARS) :: ListName, CH32SVAR1
  CHARACTER(VLCHARS), SAVE :: RstrtCopy, RunCommand
  LOGICAL, SAVE :: FCopyRstrt=.False., FOrderChecked=.False.
  LOGICAL   :: FTmp
  INTEGER(kint) :: NbLists,NStep
  INTEGER(KINT), PARAMETER :: NbExtName=4


  ! ======================================================================
  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID
  END INTERFACE
  ! ======================================================================
  
  INTEGER(KINT) :: OptGeom
  Namelist/path/OptGeom
  
  Pi=dacos(-1.0d0)

  debug=valid('EGRAD')
  if (debug) Call header("Entering Egrad")

  Print *, 'EGrad.f90, L73, IOpt=', IOpt

  Grad=0.
  E=0.

  ALLOCATE(GradCart(3*Nat))
  ALLOCATE(x(Nat),y(Nat),z(Nat))
  SELECT CASE (COORD)
  CASE ('ZMAT')
     ! In order to avoid problem with numbering and linear angles/diehedral, we convert the
     ! zmat into cartesian coordinates
     ! A remplacer par Int2Cart :-)
     Call Int2Cart(nat,indzmat,Geom,GeomCart)
  CASE ('BAKER')
     WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop"
  CASE ('MIXED')
!  write(*,*) "Coucou 4-mixed"
     Call Mixed2Cart(Nat,indzmat,geom,GeomCart)
  CASE ('CART')
!  write(*,*) "Coucou 4-cart"
     if (debug) WRITE(*,*) "Coord=CART... in egrad"
     GeomCart=reshape(Geom,(/Nat,3/))
     DO Iat=1,Nat
        WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),Geom(3*Iat-2:3*Iat)
     END DO
  CASE ('HYBRID')
!  write(*,*) "Coucou 4-hybrid"
        GeomCart=reshape(Geom,(/Nat,3/))
  CASE DEFAULT
   WRITE (*,*) "Coord=",TRIM(COORD)," not recognized. Stop"
   STOP
  END SELECT

  SELECT CASE (Prog)
     CASE ('GAUSSIAN')
! we call the Gaussian routine.
! it will return the energy and the gradient in cartesian coordinates.
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT)
        Call egrad_gaussian(E,GeomCart,GradCart)
     CASE ('MOPAC')
! we call the Mopac routine.
! it will return the energy and the gradient in cartesian coordinates.
! Gradient calculated at GeomCart geometry and stored in GradCart, which has INTENT(OUT)
        Call egrad_mopac(E,GeomCart,GradCart)
     CASE ('VASP')
        Call egrad_vasp(E,Geomcart,GradCart)
     CASE ('TURBOMOLE')
        Call egrad_turbomole(E,Geomcart,GradCart)
     CASE ('EXT')
        Call egrad_ext(E,Geomcart,GradCart)
     CASE ('TEST')
        Call egrad_test(Nat,E,Geomcart,GradCart)
     CASE ('CHAMFRE')
        Call egrad_chamfre(Nat,E,Geomcart,GradCart)
  END SELECT
  if (debug) THEN
     WRITE(*,*) 'DBG EGRAD, GradCart read'
     DO I=1,Nat
        Iat=I
        if (renum) Iat=orderInv(I)
        WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GradCart(3*I-2:3*I)
     END DO
  END IF

!  If (PROG=="VASP") STOP

  ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
  !  but we have to convert it into Zmat if COORD=Zmat
  SELECT CASE (COORD)
     CASE ("ZMAT")
        ALLOCATE(GradTmp(3*Nat))
        if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
        CALL CalcBMat_int (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome)

        if (debug) WRITE(*,*) "DBG EGRAD Calling ConvGrad_Cart2Int"
        CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)

        if (debug) WRITE(*,*) "DBG EGRAD Storing Grad"
     Grad(1)=GradTmp(4)
     Grad(2)=GradTmp(7)
! We might have troubles whith angles that are not in the [0:pi] range because,
! when converted to catesian coord, they do appear in the [0:Pi] range to gaussian...
! so that the derivative is not good, and a multiplicative factor should be added...
!
! This in fact should be taken care of in B mat calculation... 
!
     IF ((Geom(3).LE.0).OR.(Geom(3).GE.Pi)) THEN
        Grad(3)=-1.0d0*GradTmp(8)
     ELSE
        Grad(3)=GradTmp(8)
     END IF
     Idx=4
     DO I=4,Nat
        Grad(Idx)=GradTmp(3*i-2)
        Grad(Idx+1)=GradTmp(3*i-1)
        Grad(Idx+2)=GradTmp(3*i)
        Idx=Idx+3
     END DO
     DEALLOCATE(GradTmp)
   CASE ('BAKER')
        WRITE (*,*) "Coord=",TRIM(COORD),": use Egrad_baker.f90, instead of Egrad.f90. Stop"
     CASE ('MIXED')
        ALLOCATE(GradTmp(3*Nat))
        if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
        CALL CalcBMat_mixed (nat, reshape(GeomCart,(/3,Nat/),ORDER=(/2,1/)), indzmat, dzdc, massat,atome)

        if (Debug) THEN
           WRITE(*,*) "DzDc"
           DO I=1,Nat
              DO J=1,3
                 WRITE(*,'(50(1X,F12.6))') dzdc(J,I,:,:)
              END DO
           END DO
        END IF
       

        CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)
        DO I=1,Nat
           WRITE(*,'(1X,3(1X,F10.5))') GradTmp(3*I-2:3*I)
        END DO
        SELECT CASE (NCART)
           CASE (1)
              Grad(1:3)=GradTmp(1:3)
              Grad(4)=GradTmp(4)
              Grad(5)=GradTmp(7)
              Grad(6)=GradTmp(8)
              Idx=7
              IBeg=4
           CASE(2)
              Grad(1:3)=GradTmp(1:3)
              Grad(4:6)=GradTmp(4:6)
              Grad(7)=GradTmp(7)
              Grad(8)=GradTmp(8)
              Idx=9
              IBeg=4
            CASE DEFAULT
               Idx=1
               IBeg=1
            END SELECT
     DO I=IBeg,Nat
        Grad(Idx)=GradTmp(3*i-2)
        Grad(Idx+1)=GradTmp(3*i-1)
        Grad(Idx+2)=GradTmp(3*i)
        Idx=Idx+3
     END DO
     DEALLOCATE(GradTmp)
  CASE ("CART","HYBRID")
     Grad=GradCart
  CASE DEFAULT
     WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP"
     STOP
  END SELECT

  DEALLOCATE(GradCart)
  DEALLOCATE(x,y,z)

  IF (Debug) WRITE(*,*) 'DBG EGRAD grad',grad(1:NCoord)
  if (debug) Call Header("Egrad Over")

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

