Chimie Théorique » OpenPath

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

  ! 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) IGeom0=1

     IF (OptProd) IGeomF=NGeomF

     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