SUBROUTINE EGradPath(Iopt)

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

!----------------------------------------------------------------------
!  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 : ProgExe, NCoord, PROG, Nat, NGeomF, Coord, IntCoordF, &
       IndZmat, XyzGeomF, Atome, Order, OrderInv, Latr, FrozAtoms, &
       Grad, Energies, renum, dzdc, massat, Pi, NCART, OptReac, &
       OptProd,  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

  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) :: E
! Flag to see if we should print debug info
  LOGICAL :: Debug
! Flags to check if we have found the information we need
  LOGICAL, SAVE :: FReadE, FReadGrad
  LOGICAL, SAVE :: First=.TRUE.
! Flag to check if a file/directory exists
  LOGICAL :: FExists

  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

     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, &
       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

  ! 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

 END subroutine Egrad

  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))


  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

        WRITE(IOOUT,*) " Creating the POSCAR files "


        ! 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
! PFL 2013 OCt ->
! We might have to create the directories
           Line2="00"
           DO IGeom=IGeom0,IGeomF
              WRITE(Line,'(I5)') IGeom-1
              Idx=2-Len_TRIM(ADJUSTL(Line))
              FileInp=Line2(1:Idx) // TRIM(AdjustL(Line))
              INQUIRE(File=FileInp,Exist=FExists)
              If (.NOT.Fexists) THEN
                 FileInp="mkdir " // TRIM(FileInp)
                 if (debug) WRITE(*,'(1X,A)') "Creating dir: " // TRIM(FileInp)
                 CALL SYSTEM(FileInp)
              ELSE IF (DEBUG) THEN
                 WRITE(*,'(1X,A)') "Directory " // TRIM(FileInp) // " already exists"
              END IF
           END DO
        END IF   ! If (First)

        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="." // FileDelim // Line2(1:Idx) // TRIM(AdjustL(Line)) // FileDelim // "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)
        WRITE(IOOUT,*) " Launching VASP to get the Forces "
        call system(RunCommand)
        if (debug) WRITE(*,*) "Back from:", TRIM(RunCommand)

!!!!!!!!!!!!!
!
        ! We read the gradients and energies
!
        WRITE(IOOUT,*) " Reading Forces "

        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="." // FileDelim // Line2(1:Idx) // TRIM(AdjustL(Line)) // FileDelim // "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=9991) LINE
           END DO

! PFL oct 2013 ->
! The order of the things that we search might change
! so that we look for TOTEN and TOTAL-FORCE in parallel :)
           FReadE=.FALSE.
           FReadGrad=.FALSE.
           DO WHILE (.NOT.(FReadE.AND.FReadGrad))
              READ(IOTMP,'(A)',END=999,ERR=9992) LINE

              IF (INDEX(LINE,'TOTEN')>0) THEN
                 Line=Line(26:)
                 READ(LINE,*) E
                 if (debug) WRITE(*,'(1X,A,I5,A,F20.6)') "DBG Egrad_VASP E(",Igeom,")=",E
                 Energies(Igeom)=E
                 FReadE=.TRUE.
              END IF

              IF (INDEX(LINE,'TOTAL-FORCE')>0) THEN
                 READ(IOTMP,'(A)',END=999,ERR=9993) LINE
                 DO I=1,Nat
                    Iat=I
                    IF (renum) Iat=Order(I)
                    READ(IOTMP,*,END=999,ERR=9994)  Rtmp1,RTmp2,RTmp3,GradCart(3*Iat-2:3*Iat)
                 END DO
                 FReadGrad=.TRUE.
              END IF
           END DO

           Close(IoTmp)

           WRITE(IOOUT,*) " Forces read for geometry ",IGeom

           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")

        WRITE(IOOUT,*) " Creating the Input files and ListJob"
        ! 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)
        WRITE(IOOUT,*) " Launching Gaussian to compute the forces"
        call system(Line)
        IF (DEBUG) WRITE(*,*)'Back from RunCommand:'

! We read energies and gradients from the log files

        WRITE(*,*) "Gaussian over. Reading the Forces"

        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

  WRITE(*,*) "Forces and energy read for geometry ",IGeom
           ! 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)


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Serial executions
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ELSE ! matches IF (RunMode=="PARA") THEN
 ! We will launch all calculations sequentially
     ALLOCATE(GradTmp(NCoord))
     ! 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
        WRITE(IOOUT,*) "Creating Input for ", TRIM(Line)
!        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) THEN 
           WRITE(*,*) 'DBG Main, COORD=',TRIM(COORD),' IGeom=',IGeom,' Iopt=',Iopt
           WRITE(*,*) GeomTmp(1:min(NCoord,12))
        END IF

        GeomCart=0.d0
    
       ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
       ! GeomCart has INTENT(OUT)

       WRITE(IOOUT,*) "Computing and reading forces for ", TRIM(Line)
       Call EGrad(E,GeomTmp,GradTmp,NCoord,IGeom,IOpt,GeomCart) !GeomTmp=IntCoordF(IGeom,:)

    ! 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

  WRITE(IOOUT,*) " Forces all read. Moving on"


  DEALLOCATE(GeomTmp)
  DEALLOCATE(x,y,z)
  DEALLOCATE(x_k,y_k,z_k)
  if (debug) Call header('Exiting EgradPath')
  RETURN
999 WRITE(*,*) "EgradPath : We should not be here !"
   WRITE(*,*) "End of file happened too soon"
  STOP
9991 WRITE(*,*) "EgradPath : We should not be here !"
  WRITE(*,*) "EgradPath was Looking for 'EDIFF is reached'"
  STOP
9992 WRITE(*,*) "EgradPath : We should not be here !"
  WRITE(*,*) "EgradPath was Looking for 'TOTEN'"
  STOP
9993 WRITE(*,*) "EgradPath : We should not be here !"
  WRITE(*,*) "EgradPath was Looking for 'TOTAL-FORCE'"
  STOP
9994 WRITE(*,*) "EgradPath : We should not be here !"
  WRITE(*,*) "EgradPath was reading the forces "
  STOP

END SUBROUTINE EGRADPATH
