SUBROUTINE EGradPath(Iopt,Flag_Opt_Geom)

  ! This subroutine calculates the energy+gradient for all points of
  ! the path

  use Path_module, only : ProgExe, NCoord, PROG, Nat, NGeomF, Coord, IntCoordF, &
       IndZmat, XyzGeomF, Atome, Order, OrderInv, Latr, FrozAtoms, &
       Grad, Energies, renum, dzdc, massat, Pi, NCART, OptReac, &
       OptProd, IntCoordI, GeomOld_all, AtName, Unit, CalcEReac, CalcEProd
  ! GeomOld_all(NGeomF,NCoord), BTransInv (3*Nat,NCoord): used for Baker case
  use Io_module

  IMPLICIT NONE

  INTEGER(KINT), INTENT(IN) :: Iopt
  LOGICAL, INTENT(IN) :: Flag_Opt_Geom

  INTEGER(KINT) :: IGeom,IGeom0,IGeomF
  CHARACTER(LCHARS) :: Line,FileInp,RunCommand
  CHARACTER(SCHARS) :: Line2

  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:) ! NCoord
  REAL(KREAL), ALLOCATABLE :: GradCart(:) ! 3*Nat
  REAL(KREAL), ALLOCATABLE :: PathCart(:,:) ! (3*Nat,NGeomF)
  REAL(KREAL), ALLOCATABLE :: x(:), y(:), z(:)
  REAL(KREAL), ALLOCATABLE :: x_k(:), y_k(:), z_k(:)
  REAL(KREAL), ALLOCATABLE :: GeomCart(:,:) !(Nat,3)
  REAL(KREAL), ALLOCATABLE :: GeomOld_dummy(:),GeomCart_old_dummy(:,:)
  REAL(KREAL) :: E
  LOGICAL :: Debug
  LOGICAL, SAVE :: First=.TRUE.

  INTEGER(KINT) :: I, Iat, IBeg, Idx
  REAL(KREAL) :: RTmp1, RTmp2, RTmp3
  INTEGER(KINT) :: Itmp1, ITmp2

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

  debug=valid('EGradPath')

  if (debug) Call header("Entering EgradPath")

  ALLOCATE(GeomTmp(NCoord))
  ALLOCATE(x(Nat),y(Nat),z(Nat))
  ALLOCATE(x_k(Nat),y_k(Nat),z_k(Nat))
  ALLOCATE(GeomOld_dummy(NCoord))
  ALLOCATE(GeomCart_old_dummy(Nat,3))


  IF (RunMode=="PARA") THEN ! matches at L315

     ALLOCATE(GradCart(3*Nat),PathCart(3*Nat,NGeomF))


     SELECT CASE(Prog)
     CASE ("VASP")
        ! We will use the NEB routine of VASP to get all forces at once on Para8 queue

        ! First, we create all the POSCAR into the 0x directories

        ! With RunMode="Para" one cannot optimize reactants or products (due to Vasp constraints)
        ! so this has been done beforehand :-)
        IGeom0=2
        IGeomF=NGeomF-1
        PathCart=0.d0

        IF (First) THEN
           ! For the first iteration, we do write a POSCAR into 00 and 0Final
           ! Vasp needs it !
           First=.FALSE.
           IGeom0=1
           IGeomF=NGeomF
        END IF

        Line2="00"
        DO IGeom=IGeom0,IGeomF

           ! NOTE THAT HERE IGeom0 (IGeomF) IS NOT 1 (NGeomF) FOR NON-FIRST OPTIMIZATION.
           ! THIS WILL HAVE EFFECT ON XyzGeomF UPDATE IN BAKER CASE.
           WRITE(Line,'(I5)') IGeom-1
           Idx=2-Len_TRIM(ADJUSTL(Line))
           FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/POSCAR"
           if (debug) WRITE(*,*) "Creating ",TRIM(FileInp)
           OPEN(IOTMP,File=TRIM(FileInp))

           WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt
           if (debug) WRITE(*,*) Line
           SELECT CASE(Coord)
           CASE ('ZMAT')
              GeomTmp=IntCoordF(IGeom,:)
              Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
           CASE ('BAKER')
              GeomTmp=IntCoordF(IGeom,:)
              !
              IF (IOpt .EQ. 0) THEN
                 ! EgradPath(...) is called only after the call of PathCreate(...)
                 ! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always 
                 ! true if IOpt==0), THEN ExtraPol_baker(...) is called.
                 ! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting
                 ! the interpolated internal coordinates into cartesian coordinates.
                 ! Thus here we don't need to reconvert the internal coordinates again
                 ! to the cartesian coordinates. 
                 DO I=1,Nat
                    PathCart(3*I-2,IGeom) = XyzGeomF(IGeom,1,I)
                    PathCart(3*I-1,IGeom) = XyzGeomF(IGeom,2,I)
                    PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I)
                 END DO
              ELSE
                 ! XyzGeomF(NGeomF,3,Nat).
                 x_k(:) = XyzGeomF(IGeom,1,:)
                 y_k(:) = XyzGeomF(IGeom,2,:)
                 z_k(:) = XyzGeomF(IGeom,3,:)

                 !Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), &
                 !                               x_k,y_k,z_k,x,y,z)
                 ! PathCart(3*Nat,NGeomF)
                 DO I=1,Nat
                    PathCart(3*I-2,IGeom) = x(I)
                    PathCart(3*I-1,IGeom) = y(I)
                    PathCart(3*I,IGeom) = z(I)
                 END DO

                 XyzGeomF(IGeom,1,:)=x(:)
                 XyzGeomF(IGeom,2,:)=y(:)
                 XyzGeomF(IGeom,3,:)=z(:)
              END IF
           CASE ('MIXED')
              GeomTmp=IntCoordF(IGeom,:)
              Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
           CASE ('CART','HYBRID')
!!! CAUTION : PFL 29.AUG.2008 ->
              ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
              ! This should be made  consistent !!!!!!!!!!!!!!!
              GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
!!! <- CAUTION : PFL 29.AUG.2008 
            PathCart(:,IGeom)=GeomTmp
          CASE DEFAULT
            WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L134. STOP"
            STOP
        END SELECT

           IF (debug) THEN
              WRITE(*,*) "Calling PrintGeomVasp for geom",Igeom
              Call PrintGeom(Line,Nat,NCoord,PathCart(1,Igeom),Coord,IoOut,Atome,Order,OrderInv,IndZmat)
           END IF
           Call PrintGeomVasp(Line,PathCart(1,IGeom),Latr,FrozAtoms,IoTmp)

           !           IF (debug) THEN 
           !              WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD), ' Igeom=',IGeom
           !              WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
           !              WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
           !              DO Iat=1,Nat
           !                 WRITE(*,*) XyzGeomF(IGeom,1:3,Iat)
           !              END DO
           !           END IF

           IF (debug) THEN 
              WRITE(*,*) 'DBG MAin, COORD=',TRIM(COORD),' Igeom=',IGeom,' Iopt=',Iopt
              WRITE(*,*) GeomTmp(1:min(NCoord,9))
           END IF

           CLOSE(IOTMP)
        END DO ! matches DO IGeom=IGeom0,IGeomF

!!!!!!!!!!!!!!!!!!!
!
        ! We calculate the energies and gradients
!

        RunCommand=Trim(ProgExe)
        if (debug) WRITE(*,*) "Launching:", TRIM(RunCommand)
        call system(RunCommand)
        if (debug) WRITE(*,*) "Back from:", TRIM(RunCommand)

!!!!!!!!!!!!!
!
        ! We read the gradients and energies
!
        IGeom0=2
        IGeomF=NGeomF-1
        Grad=0.d0

        Line2="00"
        DO IGeom=IGeom0,IGeomF
           WRITE(Line,'(I5)') IGeom-1
           Idx=2-Len_TRIM(ADJUSTL(Line))
           FileInp=Line2(1:Idx) // TRIM(AdjustL(Line)) // "/OUTCAR"
           if (debug) WRITE(*,*) "Reading E and forces from ", TRIM(FileInp)
           OPEN(IOTMP,File=TRIM(FileInp))

           ! 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,I5,A,F20.6)') "DBG Egrad_VASP E(",Igeom,")=",E
           Energies(Igeom)=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

           Close(IoTmp)

           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

           ! We convert it into the right coordinate system
           ! This should be put into a subroutine (because it is also in egradient)
           SELECT CASE (COORD)
           CASE ("ZMAT")
              ALLOCATE(GradTmp(3*Nat))
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
              CALL CalcBmat_int(nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome)
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)

              Grad(IGeom,1)=GradTmp(4)
              Grad(IGeom,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 calc Bmat ...
              !
              IF ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN
                 Grad(IGeom,3)=-1.0d0*GradTmp(8)
              ELSE
                 Grad(IGeom,3)=GradTmp(8)
              END IF
              Idx=4
              DO I=4,Nat
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
                 Idx=Idx+3
              END DO
              DEALLOCATE(GradTmp)
              ! 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      
           CASE ('BAKER')
              ALLOCATE(GradTmp(3*Nat))
              GradTmp=0.d0
              ! Below is to be corrected.
              ! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN'
              !DO I=1, 3*Nat
              ! here GradTmp and Grad are gradients in Baker coordinates
              ! Not implemented  IF ((PROG=="VASP").AND.(RunMode=="PARA"))
              ! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I)
              !END DO
              Grad(IGeom,:) = GradTmp
              DEALLOCATE(GradTmp)
           CASE ('MIXED')
              ALLOCATE(GradTmp(3*Nat))
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
              CALL CalcBMat_mixed (nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome)
              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(IGeom,1:3)=GradTmp(1:3)
                 Grad(IGeom,4)=GradTmp(4)
                 Grad(IGeom,5)=GradTmp(7)
                 Grad(IGeom,6)=GradTmp(8)
                 Idx=7
                 IBeg=4
              CASE(2)
                 Grad(IGeom,1:3)=GradTmp(1:3)
                 Grad(IGeom,4:6)=GradTmp(4:6)
                 Grad(IGeom,7)=GradTmp(7)
                 Grad(IGeom,8)=GradTmp(8)
                 Idx=9
                 IBeg=4
              CASE DEFAULT
                 Idx=1
                 IBeg=1
              END SELECT
              DO I=IBeg,Nat
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
                 Idx=Idx+3
              END DO
              DEALLOCATE(GradTmp)
           CASE ("CART","HYBRID")
              Grad(IGeom,:)=GradCart
           CASE DEFAULT
              WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP"
              STOP
           END SELECT

           IF (debug) THEN
              Line='DBG Path, GradTmp'
              Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat)
           END IF
        END DO ! matches  DO IGeom=IGeom0,IGeomF

     CASE ("GAUSSIAN")
        ! We create the output files.
        IGeom0=2
        IGeomF=NGeomF-1
        PathCart=0.d0

        If (OptReac) IGeom0=1
        If (OptProd) IGeomF=NGeomF

        OPEN(IOTMP2,File="ListJob",STATUS="REPLACE")

        DO Igeom=IGeom0,IGeomF
           ! We first convert the geometry into a Cartesian one
           SELECT CASE(Coord)
           CASE ('ZMAT')
              GeomTmp=IntCoordF(IGeom,:)
              Call Int2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
           CASE ('BAKER')
              GeomTmp=IntCoordF(IGeom,:)
              !
              IF (IOpt .EQ. 0) THEN
                 ! EgradPath(...) is called only after the call of PathCreate(...)
                 ! and in PathCreate, IF (MOD(IOpt,IReparam).EQ.0) (which is always 
                 ! true if IOpt==0), THEN ExtraPol_baker(...) is called.
                 ! In ExtraPol_baker(...), XyzGeomF(...) is calculated upon converting
                 ! the interpolated internal coordinates into cartesian coordinates.
                 ! Thus here we don't need to reconvert the internal coordinates again
                 ! to the cartesian coordinates. 
                 DO I=1,Nat
                    PathCart(3*I-2,IGeom) = XyzGeomF(IGeom,1,I)
                    PathCart(3*I-1,IGeom) = XyzGeomF(IGeom,2,I)
                    PathCart(3*I,IGeom) = XyzGeomF(IGeom,3,I)
                 END DO
              ELSE

                 ! XyzGeomF(NGeomF,3,Nat).
                 x_k(:) = XyzGeomF(IGeom,1,:)
                 y_k(:) = XyzGeomF(IGeom,2,:)
                 z_k(:) = XyzGeomF(IGeom,3,:)

                 !Call ConvertBakerInternal_cart(GeomOld_all(IGeom,:),IntCoordI(IGeom,:), &
                 !                               x_k,y_k,z_k,x,y,z)
                 ! PathCart(3*Nat,NGeomF)
                 DO I=1,Nat
                    PathCart(3*I-2,IGeom) = x(I)
                    PathCart(3*I-1,IGeom) = y(I)
                    PathCart(3*I,IGeom) = z(I)
                 END DO

                 XyzGeomF(IGeom,1,:)=x(:)
                 XyzGeomF(IGeom,2,:)=y(:)
                 XyzGeomF(IGeom,3,:)=z(:)
              END IF
           CASE ('MIXED')
              GeomTmp=IntCoordF(IGeom,:)
              Call Mixed2cart(nat,indzmat,geomTmp,PathCart(1,IGeom))
           CASE ('CART','HYBRID')
!!! CAUTION : PFL 29.AUG.2008 ->
              ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
              ! This should be made  consistent !!!!!!!!!!!!!!!
              GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
!!! <- CAUTION : PFL 29.AUG.2008 
              PathCart(:,IGeom)=GeomTmp
           CASE DEFAULT
              WRITE(*,*) "Coord=",TRIM(Coord)," not recognized. EgradPath L75.STOP"
              STOP
           END SELECT

           WRITE(Line,'(I5)') IGeom
           FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".com"
           WRITE(IOTMP2,'(1X,A)') TRIM(FileInp)

           OPEN(IOTMP,File=TRIM(FileInp))
           Current => Gauss_root
           DO WHILE (ASSOCIATED(Current%next))
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
              WRITE(*,'(1X,A)') Trim(current%line)
              Current => current%next
           END DO

           WRITE(IOTMP,*) 

           Current => Gauss_comment
           DO WHILE (ASSOCIATED(Current%next))
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
              Current => current%next
           END DO
           !  WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)"),Igeom,NgeomF,Iopt

           WRITE(IOTMP,*) 
           WRITE(IOTMP,*) Trim(Gauss_charge)

           DO I=1,Nat
              If (renum) THEN
                 Iat=Order(I)
                 WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)) &
                      ,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom)  &
                      ,TRIM(Gauss_Paste(I))
              ELSE
                 Iat=OrderInv(I)
                 WRITE(IOTMP,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)) &
                 ,PathCart(Iat,IGeom),PathCart(Iat+Nat,IGeom),PathCart(Iat+2*Nat,IGeom)  &
                      , TRIM(Gauss_Paste(Iat))
              END IF
           END DO

           WRITE(IOTMP,*) 

           Current => Gauss_End
           DO WHILE (ASSOCIATED(Current%next))
              WRITE(IOTMP,'(1X,A)') Trim(current%line)
              Current => current%next
           END DO

           WRITE(IOTMP,*) 
           CLOSE(IOTMP)

        END DO

        Close(IOTMP2)

! We launch the parallel calculations

        Line=Trim(ProgExe) 
        IF (DEBUG) WRITE(*,*)'RunCommand:',TRIM(Line)
        call system(Line)
        IF (DEBUG) WRITE(*,*)'Back from RunCommand:'

! We read energies and gradients from the log files

        DO IGeom=IGeom0,IGeomF
           WRITE(Line,'(I5)') IGeom
           FileInp="PathPar_" // TRIM(ADJUSTL(Line)) // ".log"
           OPEN(IOTMP,File=TRIM(FileInp))

  ! We first search for the forces
  READ(IOTMP,'(A)') LINE
  DO WHILE (INDEX(LINE,'Forces (Hartrees/Bohr)')==0)
     ! here, we have a problem, because there are so many ways to write E for Gaussian :-)
     ! but we do not really need E except if we want to weight the different points
     ! on the path... that will be done later :-)
     ! And I should use ConvGaussXmol ...
! PFL 3rd Jan 2008 ->
! I have finally upgraded the energy reading phase so that it should be able to read 
! many formats corresponding to HF, DFT, MP2, AM1, PM3, ONIOM with or w/o PCM
!
     
! For MP2, energy is after the last =
     IF ((Line(2:9)=="E2 =    ")) THEN
      Itmp1=Index(LINE,"=",BACK=.TRUE.)+1
    READ(LINE(Itmp1:),*) e
      END IF
! For other methods, it is after the first =
     IF ((LINE(2:9)=='Energy= ').OR.(Line(2:9)=="SCF Done").OR.(Line(2:9)=="ONIOM: e") &
   .OR.(Line(2:9)==" with al") &
         .OR.(Line(2:31)=="Counterpoise: corrected energy")) THEN
      Itmp1=Index(LINE,"=")+1
    READ(LINE(Itmp1:),*) e
  END IF
! <- PFL 3 Jan 2008
     READ(IOTMP,'(A)') LINE
  END DO
  READ(IOTMP,'(A)') LINE
  READ(IOTMP,'(A)') LINE
  DO I=1,Nat
     Iat=I
     IF (renum) Iat=Order(I)
     READ(IOTMP,*)  Itmp1,ITmp2, GradCart(3*Iat-2:3*Iat)
  END DO

! Gaussian gives the Forces in ua/a0, so we convert it into the 
! gradient in ua/Angstrom 
  gradCart=-1.0d0*Unit*gradCart

  CLOSE(IOTMP)
  Energies(IGeom)=E
  
           ! We convert it into the right coordinate system
           ! This should be put into a subroutine (because it is also in egradient)
           SELECT CASE (COORD)
           CASE ("ZMAT")
              ALLOCATE(GradTmp(3*Nat))
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_int"
              CALL CalcBMat_int (nat, PathCart(1:3*Nat,Igeom), indzmat, dzdc, massat,atome)
              CALL ConvGrad_Cart2Int(Nat,dzdc,IndZmat,GradCart,GradTmp)

              Grad(IGeom,1)=GradTmp(4)
              Grad(IGeom,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 ((IntCoordF(Igeom,3).LE.0).OR.(IntCoordF(Igeom,3).GE.Pi)) THEN
                 Grad(IGeom,3)=-1.0d0*GradTmp(8)
              ELSE
                 Grad(IGeom,3)=GradTmp(8)
              END IF
              Idx=4
              DO I=4,Nat
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
                 Idx=Idx+3
              END DO
              DEALLOCATE(GradTmp)
              ! 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      
           CASE ('BAKER')
              ALLOCATE(GradTmp(3*Nat))
              GradTmp=0.d0
              ! Below is to be corrected.
              ! This is inside 'IF ((PROG=="VASP").AND.(RunMode=="PARA")) THEN'
              !DO I=1, 3*Nat
              ! here GradTmp and Grad are gradients in Baker coordinates
              ! Not implemented  IF ((PROG=="VASP").AND.(RunMode=="PARA"))
              ! GradTmp(:) = GradTmp(:) + BTransInv(:,I)*GradCart(I)
              !END DO
              Grad(IGeom,:) = GradTmp
              DEALLOCATE(GradTmp)
           CASE ('MIXED')
              ALLOCATE(GradTmp(3*Nat))
              if (debug) WRITE(*,*) "DBG EGRAD Calling CalcBMat_mixed"
              CALL CalcBMat_mixed(nat,PathCart(1:3*Nat,IGeom), indzmat, dzdc, massat,atome)
              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(IGeom,1:3)=GradTmp(1:3)
                 Grad(IGeom,4)=GradTmp(4)
                 Grad(IGeom,5)=GradTmp(7)
                 Grad(IGeom,6)=GradTmp(8)
                 Idx=7
                 IBeg=4
              CASE(2)
                 Grad(IGeom,1:3)=GradTmp(1:3)
                 Grad(IGeom,4:6)=GradTmp(4:6)
                 Grad(IGeom,7)=GradTmp(7)
                 Grad(IGeom,8)=GradTmp(8)
                 Idx=9
                 IBeg=4
              CASE DEFAULT
                 Idx=1
                 IBeg=1
              END SELECT
              DO I=IBeg,Nat
                 Grad(IGeom,Idx)=GradTmp(3*i-2)
                 Grad(IGeom,Idx+1)=GradTmp(3*i-1)
                 Grad(IGeom,Idx+2)=GradTmp(3*i)
                 Idx=Idx+3
              END DO
              DEALLOCATE(GradTmp)
           CASE ("CART","HYBRID")
              Grad(IGeom,:)=GradCart
           CASE DEFAULT
              WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in EgradPath.L257.STOP"
              STOP
           END SELECT

           IF (debug) THEN
              Line='DBG Path, GradTmp'
              Call PrintGeom(Line,Nat,NCoord,reshape(GRad(Igeom,1:NCoord),(/Ncoord/)),Coord,6,Atome,Order,OrderInv,IndZmat)
           END IF

        END DO

     CASE DEFAULT
        WRITE(*,*) "Prog=",TRIM(Prog)," not yet supported for Para mode."
        STOP
     END SELECT

     DEALLOCATE(GradCart,PathCart)

  ELSE ! matches IF (RunMode=="PARA") THEN
 ! We will launch all calculations sequentially
     ALLOCATE(GradTmp(NCoord))
     GeomOld_dummy=0.d0 ! Internal coordinates
     GeomCart_old_dummy=0.d0
     ! We have the new path, we calculate its energy and gradient
     IGeom0=2
     IGeomF=NGeomF-1
     IF (OptReac.OR.CalcEReac) IGeom0=1
     IF (OptProd.OR.CalcEProd) IGeomF=NGeomF
     CalcEReac=.FALSE.
     CalcEProd=.FALSE.
     ALLOCATE(GeomCart(Nat,3))

     DO IGeom=IGeom0,IGeomF
        WRITE(Line,"('Geometry ',I3,'/',I3,' for iteration ',I3)") Igeom,NgeomF,Iopt
        if (debug) WRITE(*,*) Line
        SELECT CASE(Coord)
        CASE ('ZMAT','MIXED')
           GeomTmp=IntCoordF(IGeom,:)
        CASE ('BAKER')
           GeomTmp=IntCoordF(IGeom,:)
        CASE ('CART','HYBRID')
!!! CAUTION : PFL 29.AUG.2008 ->
           ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
           ! This should be made  consistent !!!!!!!!!!!!!!!
           GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
           !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
!!! <- CAUTION : PFL 29.AUG.2008
        CASE DEFAULT
           WRITE(*,*) "Coord=",TRIM(Coord)," not recognized in EgradPath. L299"
           STOP
        END SELECT

        !if (debug) WRITE(*,*) "Calling PrintGeom"
        !Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)

        !           IF (debug) THEN 
        !              WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD), ' Igeom=',IGeom
        !              WRITE(*,'(3(1X,F9.4))') Reshape(Reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
        !              WRITE(*,'(3(1X,F9.4))') Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
        !              DO Iat=1,Nat
        !                 WRITE(*,*) XyzGeomF(IGeom,1:3,Iat)
        !              END DO
        !           END IF

        IF (debug) THEN 
           WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD),' IGeom=',IGeom,' Iopt=',Iopt
           WRITE(*,*) GeomTmp(1:min(NCoord,12))
        END IF

        GeomCart=0.d0
        IF(IOpt .EQ. 1) Then
           Print *, 'GeomTmp='
           WRITE(*,'(12(1X,F6.3))') GeomTmp
           Print *, 'GeomCart='
           WRITE(*,'(12(1X,F6.3))') GeomCart
      END IF
    
    IF (COORD.EQ.'BAKER') THEN
       ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
       ! GeomCart has INTENT(OUT)
       Call EGrad_baker(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom, &
                  GeomOld_dummy,GeomCart_old_dummy)
    ELSE
       Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart,Flag_Opt_Geom) !GeomTmp=IntCoordF(IGeom,:)
        END IF

    ! Egrad calls ConvertBakerInternal_cart to convert GeomTmp=IntCoordF(IGeom,:) into 
    ! GeomCart (Cartesian Coordinates) so that energy and gradients can be calculated 
    ! for each geometry.
        XyzGeomF(IGeom,1,:)=GeomCart(:,1)
        XyzGeomF(IGeom,2,:)=GeomCart(:,2)
        XyzGeomF(IGeom,3,:)=GeomCart(:,3)
      
        Energies(IGeom)=E

        if (debug) THEN
           Line='DBG Path, GradTmp'
           Call PrintGeom(Line,Nat,NCoord,GRadTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
        END IF
        Grad(IGeom,:)=GradTmp
     END DO ! matches DO IGeom=IGeom0,IGeomF
     DEALLOCATE(GradTmp,GeomCart)
  END IF ! matches IF (RunMode=="PARA") THEN AFTER ELSE

  DEALLOCATE(GeomTmp)
  DEALLOCATE(x,y,z)
  DEALLOCATE(x_k,y_k,z_k)
  DEALLOCATE(GeomOld_dummy,GeomCart_old_dummy)
  if (debug) Call header('Exiting EgradPath')
  RETURN
999 WRITE(*,*) "EgradPath : We should not be here !"
  STOP

END SUBROUTINE EGRADPATH
