Chimie Théorique » Opt'n Path

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