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

  ! 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,BTransInv, XyzGeomI, renum, &
       GeomOld_all,BTransInvF,BTransInv_local,UMatF,UMat_local &
       , BprimT,a0,ScanCoord, Coordinate,NPrim,BMat_BakerT,BBT,BBT_inv &
      , Order,OrderInv, XPrimitiveF

  ! 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)
!!! Optional, just for geometry optimization with Baker coordinates
  REAL(KREAL), INTENT (IN), OPTIONAL :: GeomCart_old(Nat,3)
  REAL(KREAL), INTENT (INOUT), OPTIONAL :: GeomOld(NCoord)
! FOptGeom is a flag indicating if we are doing a geom optimization
! it can be omitted so that we use a local flag: Flag_Opt_Geom
  LOGICAL,INTENT (IN), OPTIONAL :: FOptGeom
! if FOptGeom is given Flag_Opt_Geom=FOptGeom, else Flag_Opt_Geom=F
  LOGICAL  :: Flag_Opt_Geom

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

  LOGICAL :: debug
  REAL(KREAL), ALLOCATABLE :: GradTmp(:)
  REAL(KREAL), ALLOCATABLE :: GradCart(:)
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
  REAL(KREAL) ::  Pi

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

  INTEGER(KINT), PARAMETER :: NbExtName=4

  REAL(KREAL), ALLOCATABLE :: XPrimRef(:), XPrim(:) ! NPrim

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!! PFL Debug
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:,:) !3,Nat
  LOGICAL ::  FALOBBT,FALOBPrimt,FAloBBTInv
  LOGICAL :: DebugPFL

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

  Pi=dacos(-1.0d0)

  if (present(FOptGeom)) THEN
     Flag_Opt_Geom=FOptGeom
  ELSE
     Flag_Opt_Geom=.FALSE.
  END IF

  debug=valid('EGRAD')
  debugPFL=valid('bakerPFL')
  if (debug) Call Header("Entering Egrad")

  if (debug) Print *, 'EGrad.f90, L73, IOpt=', IOpt, ' Flag_Opt_Geom=',Flag_Opt_Geom

  Grad=0.
  E=0.


  ALLOCATE(GradCart(3*Nat))
  ALLOCATE(x(Nat),y(Nat),z(Nat))
  ALLOCATE(XPrimRef(NPrim),XPrim(NPrim))

  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')
     XPrimRef=XPrimitiveF(IGeom,:)
     IF (Flag_Opt_Geom) THEN
        IF (IOpt .EQ. 1) THEN
           GeomCart(:,1) = XyzGeomI(IGeom,1,:)
           GeomCart(:,2) = XyzGeomI(IGeom,2,:)
           GeomCart(:,3) = XyzGeomI(IGeom,3,:)
        ELSE
           if (debug) WRITE(*,*) 'Before Call ConvertBakerInternal_cart; L125, Egrad.f90'
           if (.NOT.present(GeomOld)) THEN
              WRITE(*,*) "ERROR: GeomOld not passed as an argument"
              STOP
           END IF
           WRITE(*,*) 'GeomOld='
           WRITE(*,'(12(1X,F6.3))') GeomOld
           WRITE(*,*) 'Geom='
           WRITE(*,'(12(1X,F6.3))') Geom
           if (.NOT.present(GeomCart_Old)) THEN
              WRITE(*,*) "ERROR: GeomCart_Old not passed as an argument"
              STOP
           END IF

           WRITE(*,*) 'GeomCart_old='
           WRITE(*,'(20(1X,F6.3))') GeomCart_old(1:Nat,1:3)
           Call ConvertBakerInternal_cart(GeomOld,Geom,GeomCart_old(:,1), &
                GeomCart_old(:,2),GeomCart_old(:,3),x,y,z,XPrim,XPrimRef)
           GeomCart(:,1) = x(:)
           GeomCart(:,2) = y(:)
           GeomCart(:,3) = z(:)
        END IF
     ELSE
        IF (IOpt .EQ. 0) THEN  
           GeomCart(:,1) = XyzGeomF(IGeom,1,:)
           GeomCart(:,2) = XyzGeomF(IGeom,2,:)
           GeomCart(:,3) = XyzGeomF(IGeom,3,:)
        ELSE
           ! IntCoordF(NGeomF,NCoord), XyzGeomF(NGeomF,3,Nat)
           ! Geom has to be converted into x,y,z
           ! GeomOld_all and XyzGeomF are old geometry in internal and cartesian coords resp.
           DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart()
              BTransInv_local(I,:) = BTransInvF(IGeom,I,:)
              UMat_local(:,I) = UMatF(IGeom,:,I)
           END DO
           if (debug) Print *, 'I am before Call ConvertBakerInternal_cart; L160, Egrad.f90'
           WRITE(*,*) 'GeomOld='
           WRITE(*,'(12(1X,F6.3))') GeomOld_all(IGeom,:)
           WRITE(*,*) 'Geom='
           WRITE(*,'(12(1X,F6.3))') Geom
           WRITE(*,*) 'DBG Egrad L165 GeomCart_old='
           DO I=1,Nat
              WRITE(*,'(3(1X,F12.8))') XyzGeomF(IGeom,:,I)
           END DO

           Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),Geom,XyzGeomF(IGeom,1,:), &
          XyzGeomF(IGeom,2,:),XyzGeomF(IGeom,3,:),x,y,z,XPrim,XPrimRef)
           if (debug) Print *, 'I am after Call ConvertBakerInternal_cart; L111, Egrad.f90'
           DO I=1,NCoord ! these variables are used in ConvertBakerInternal_cart()
              BTransInvF(IGeom,I,:) = BTransInv_local(I,:)
           END DO
           ! GeomCart is actually XyzCoord and has INTENT(OUT) in this subroutine.
           GeomCart(:,1) = x(:)
           GeomCart(:,2) = y(:)
           GeomCart(:,3) = z(:)
        END IF ! matches IF (IOpt .EQ. 0) THEN
     END IF
  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
  !WRITE(*,*) Nat
  !WRITE(*,*) 'GeomCart in Egrad_baker'
  !DO I=1,Nat
  !  WRITE(*,'(1X,A,3(1X,F12.8))') AtName(I),GeomCart(I,1:3)
  !END DO

  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')
     ! We have the gradient in Cartesian coordinates... which is ok for COORD=Cart, Hybrid
     ! but we have to convert it into internal coordinates if COORD=Baker	 
!!!!!!!!!!!!!!!!!!!!
     !
     !   PFL Debug
     !
!!!!!!!!!!!!!!!!!!!!
     if (debugPFL) THEN
        WRITE(*,*) "DBG PFL Egrad_Baker... computing dzdc"
     if (.not.ALLOCATED(indzmat)) THEN
        ALLOCATE(indzmat(Nat,5))
     IndZmat=0
     Indzmat(1,1)=1
     IndZmat(2,1)=2
     IndZmat(2,2)=1
     IndZmat(3,1)=3
     IndZmat(3,2)=2
     IndZmat(3,3)=1
     IndZmat(4,1)=4
     IndZmat(4,2)=3
     IndZmat(4,3)=2
     IndZmat(4,4)=1
     END IF
     IF (.NOT.ALLOCATED(DzDc)) THEN
        ALLOCATE(DzDc(3,nat,3,Nat))
     END IF
     if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat))                                            
     DO I=1,Nat
     GeomTmp(1,I)=GeomCart(OrderInv(I),1)
     GeomTmp(2,I)=GeomCart(OrderInv(i),2)
     GeomTmp(3,I)=GeomCart(OrderInv(i),3)
     END DO

     DzDc=0.d0
     CALL CalcBMat_int (nat, GeomTmp, indzmat, dzdc, massat,atome)     

     END IF ! if (debugPFL)
!!!!!!!!!!!!!!!!!!!!!!!!!
     ! Debugging PFL
!!!!!!!!!!!!!!!!!!!!!!!!

     WRITE(*,*) "BTransInv_local Trans that is used originally"
     DO J=1,3*Nat
        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
     END DO

     WRITE(*,*) "Calculating actual values using GeomCart"
     if (.NOT.Allocated(GeomTmp)) Allocate(GeomTmp(3,Nat))
     GeomTmp(1,:)=GeomCart(1:Nat,1)
     GeomTmp(2,:)=GeomCart(1:Nat,2)
     GeomTmp(3,:)=GeomCart(1:Nat,3)

     ! Computing B^prim:
     FAloBBT=.NOT.ALLOCATED(BBT)
     FAloBBTInv=.NOT.ALLOCATED(BBT_inv)
     FAloBPrimT=.NOT.ALLOCATED(BprimT)
     if (FAloBPrimT) ALLOCATE(Bprimt(3*Nat,NPrim))
     if (FAloBBT) ALLOCATE(BBt(NCoord,NCoord))
     if (FAloBBTInv) ALLOCATE(BBt_inv(NCoord,NCoord))
     BprimT=0.d0
     ScanCoord=>Coordinate
     I=0
     DO WHILE (Associated(ScanCoord%next))
        I=I+1
        SELECT CASE (ScanCoord%Type)
        CASE ('BOND') 
           CALL CONSTRAINTS_BONDLENGTH_DER(Nat,ScanCoord%at1,ScanCoord%AT2,GeomTmp,BprimT(1,I))
        CASE ('ANGLE')
           CALL CONSTRAINTS_BONDANGLE_DER(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,GeomTmp,BprimT(1,I))
        CASE ('DIHEDRAL')
           CALL CONSTRAINTS_TORSION_DER2(Nat,ScanCoord%At1,ScanCoord%AT2,ScanCoord%At3,ScanCoord%At4,GeomTmp,BprimT(1,I))
        END SELECT
        ScanCoord => ScanCoord%next
     END DO

     if (debug) THEN
        WRITE(*,*) "Bprim "
        DO J=1,3*Nat
           WRITE(*,'(50(1X,F12.6))') BprimT(J,:)
        END DO
     END IF

     BMat_BakerT = 0.d0
     DO I=1,NCoord
        DO J=1,NPrim
           !BprimT is transpose of B^prim.
           !B = UMat^T B^prim, B^T = (B^prim)^T UMat
           BMat_BakerT(:,I)=BMat_BakerT(:,I)+BprimT(:,J)*UMat_local(J,I)
        END DO
     END DO

     BBT = 0.d0
     DO I=1,NCoord
        DO J=1,3*Nat   ! BBT(:,I) forms BB^T
           BBT(:,I) = BBT(:,I) + BMat_BakerT(J,:)*BMat_BakerT(J,I)
        END DO
     END DO

     BBT_inv = 0.d0

     Call GenInv(NCoord,BBT,BBT_inv,NCoord)

!     ALLOCATE(V(NCoord,NCoord))
!     tol = 1e-10
!     ! BBT is destroyed by GINVSE
!     Call GINVSE(BBT,NCoord,NCoord,BBT_inv,NCoord,NCoord,NCoord,tol,V,NCoord,FAIL,IOOUT)
!     DEALLOCATE(V)
!     IF(Fail) Then
!        WRITE (*,*) "Failed in generalized inverse, L220, ConvertBaker...f90"
!        STOP
!     END IF

     !Print *, 'BBT_inv='
     DO J=1,NCoord
        ! WRITE(*,'(50(1X,F10.2))') BBT_inv(:,J)
        ! Print *, BBT_inv(:,J)
     END DO
     !Print *, 'End of BBT_inv'

     ! Calculation of (B^T)^-1 = (BB^T)^-1B:
     BTransInv_local = 0.d0
     DO I=1,3*Nat
        DO J=1,NCoord
           BTransInv_local(:,I)=BTransInv_local(:,I)+BBT_inv(:,J)*BMat_BakerT(I,J)
        END DO
     END DO

     Print *, 'BMat_BakerT='
     DO J=1,3*Nat
        WRITE(*,'(50(1X,F12.8))') BMat_BakerT(:,J)
     END DO

     WRITE(*,*) "BTransInv_local Trans corresponding to GeomCart"
     DO J=1,3*Nat
        WRITE(*,'(50(1X,F12.8))') BTransInv_local(:,J)
     END DO

     if (FAloBPrimT) DEALLOCATE(Bprimt)
     if (FAloBBT) DEALLOCATE(BBt)
     if (FAloBBTInv) DEALLOCATE(BBt_inv)

     Grad=0.d0
     DO I=1, 3*Nat
        ! Here Grad is gradient in Baker coordinates. GradCart(3*Nat)
        Grad(:) = Grad(:) + BTransInv_local(:,I)*GradCart(I) !BTransInv_local = BTransInv in Opt_Geom.f90
     END DO
     !WRITE(*,*) 'In Egrad.f90, L263, Gradient:' 
     !WRITE(*,'(12(1X,F8.4))') Grad(1:NCoord)


     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

