! This programs generates and optimizes a reaction path
! between at least two endpoints.
!
! It uses NAMELIST for input, see below !
!
! v3.0
! First version entirely in Fortran.
! It uses routines from Cart (path creation and parameterisation) 
! and from IRC06, especially routines for point optimization.
!
! TO DO:
! 1) Possibility to read and/or calculate some hessian structures
! 2) Use of internal coordinate to estimate the hessian for stable
! species in a much better way...
!
!!!!!!!!!!!!!!!
! v3.1
! The Hessian update includes neighbour points.
!
! v3.2
! We add the option Hinv that uses the BFGS update of the Hessian inverse
!
! v3.3
! We add the full option for ZMAT coordinates, where all is done using
! internal coordinates.
! v3.31
! The step size is controlled using the step norm and not the maximum step
! component.
! v3.5
! Big modifications of Calc_zmat_const_frag in order to be able
! to treat VASP geometries. 
! For now, only XYZ geometries are read and written. VASP interface
! has to be written.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! v3.6
! Some minor modif in Calc_zmat_const_frag to have nicer programming
! and a new mode added: Coord=MIXED where part of the system
! is extrapolated in cartesian coordinates and the rest of it in zmat.
! it might be Useful for VASP, or when lots of atoms are frozen
! by default, when Coord=MIXED, all frozen atoms are described in cartesian
! coordinates.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! v3.61
! We move further and further to VASP program. New variables added to the 'path' namelist:
! input: what is the program used for the geometries ?
!     Possible values: xyz (or cart) for Xmol format; VASP. Xyz should be used for Gaussian
! prog: what is the program used for the calculations ?
!     Possible values for now: gaussian, vasp.
! For now (02/2007), when prog='vasp' is used, only the initial path is created.
! That means that we have only added input/output subroutines :-)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! v3.7
! We try to improve the optimization steps. For that, we first modify many routines
! so that the new step is calculated in the 'free' degrees of freedom, as it is done
! in IRC07/ADF for example. That will allow us to impose constraints in an easier way.
!
! v3.701
! Add a new option for the initial extrapolation of the Hessian
!
!v3.71
! The optimization step has been improved... but the vfree part has not yet been included.
! This has still to be done :-)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! v3.8
! I have included part of the vfree procedure: for now freemv returns the Id
! matrix, that is used to project out the tangent vector.
! This improves the optimization a lot, but I will have to implement
! a real freemw subroutine to have constraints there !
! v3.81 Technical version
! I have added a new program 'Ext'. When this one is used, Path creates a Ext.in file,
! calls Ext.exe file to get a Ext.out file that contains the Energy and the gradient.
! Format:
! Ext.in : Xmol format. The second line (comment in Xmol) contains the COORD name.
! Ext.out: Energy on firt line (1X,D25.15), then gradient in CARTESIAN coordinates
!  (3(1X,D25.15)) format. Natoms lines. 
! TO DO:  Gradient could be given in COORD format, ie cartesian for COORD=CART, XYZ or HYBRID
! ZMAT for CORD=ZMAT (in that case, 6 zeros are there at the begining !)
! For MIXED, it is equivalent to CART :-)
! v3.811
! I have added a prog="TEST" option, where egradient calls egrad_test subroutine.
! It then calculates the energy and gradient as it wants and it returns the cartesian 
! gradient. The purpose is to have an analytical approximation of the PES of the system
! under study to gain some time in testing the program. This is not documented :-)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  v3.9
! In order to be fully interfaced with VASP, I have to change the architecture
! of the program. In particular, with VASP, one can use a Para8+NEB routine
! of VASP to calculate energy+forces for all points of the path.
! v3.91
! One small modification: in the cart list, one can freeze a whole sequence by
! entering a negative number for the last atom of the sequence
! for example : cart=1 -160 165 -170 will define atoms from 1 to 160 and from
! 165 to 170 (included) as cartesian atoms.  
! Of course, Same applies to Frozen
! v3.92 / v3.93
! Slight modifications. One noticeable: when using Vasp, the program analyses the
! displacements to suggest which atoms could be described in CART.
! v3.94 
! VaspMode=Para now implemented. One needs to be careful in defining 'ProgExe'.
! In particular, one must put the mpirun command in ProgExe !!! There is NO test.
! v3.95
! When using internal coordinates (zmat, hybrid or mixes), Path calculates the number
! of fragments for each geometry. The reference one (ie the one used to compute the internal
! coordinates) is now the one with the least fragments.
! user can also choose one geometry by specifying IGeomRef in the Path namelis.
! v3.96
! I have added a test for the geometry step: it does not allow valence angle to 
! be negative or bigger than 180. I hope that this might help geometry optimization
! when dealing with nearly linear molecules.
! v3.961
! The options cart and frozen are now LOGICAL:
! Fcart=T indicates that a &cartlist namelist should be read !
! Ffrozen=T indicates that a &frozenlist namelist should be read !
! TO DO: Autocart=T indicates that the atoms described in cartesian coordiantes should
! be found automatically
! v3.97
! Baker has been implemented by P. Dayal: we are interpolating Xint for now.
! TO do: interpolate Umat ?
! -----
! v3.971
! AutoCart=T is being implemented. Goal: having a 'black box' version for Vasp users.
! Vmd: True if user wants to watch the path using VMD. Only used for VASP for now.
!------------
! v3.98
! Calc_zmat has been replaced by Calc_zmat_frag which is more robust.
! Still missing: linear molecules and some tests.
!
! v3.981
! * Added a non document flag in .valid file that switch on the numerical
! calculation of the hessian instead of the update.
! * Added HesUpd variable in the namelist. Default value is MS for path
!   optimization because it does not imposes a positive hessian, and BFGS
!  for geometry optimization
! v 4.0
! * This version has been made public within the CARTE project.
! * Added:
! - Step_method to choose the optimization method: RFO or Diis
! - DynMaxStep: if TRUE, the maximum size of a step is updated at each step. 
!             If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. 
!             It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]
!
!!!!!!!!!!!!!!!!!!!!!!
! v4.1
! * Para mode has been partly implemented for Gaussian.
!  VaspMode is thus now 'RunMode' and can be Serial or Para
! Only Gaussian and VASP are supported for Para mode now.
! v4.12 
! Some bugs in unconstrained zmat construction are corrected.
! v4.13
! Prakash has implemented the GEDIIS optimization procedure.
! v4.14
! There is a small problem for some inputs in VASP, where the
! last geometry is not the one given by the user. 
! It seems to come from the fact that some people try to freeze
! some unfrozen atoms
! So some tests have been added to check that frozen atoms do not move.
!!!
! v4.15
! There were some bugs when reading and expanding cartlist and frozenlist
! I have changed the way frozen atoms are compared by adding partial alignment.
! v4.16
! Some bugs were corrected.
! v4.17
! Added MOPAC as a program.
! v4.175
! Added Three more parameters for defining the input and output names for
! the calculations.
! CalcName: Prefix for the files used for the energy and gradient calculations
! ISuffix: Suffix for the input file
! OSuffix: suffix for the output file.
! This maters for Gaussian for example.
!
! v4.177 (P) 2010 Nov 22
! - We add TurboMole as a new external code.
! - Subtle change for Input: the default is XYZ whatever the code, except
! for VASP. 
!
! v4.178 (P) 2010 Nov 28
! - We add the possibility to call CHAMFRE reactive force field using 
! prog=TEST2 or CHAMFRE
!


PROGRAM PathOpt

  use VarTypes
  use Path_module
  use Io_module

  IMPLICIT NONE

  CHARACTER(132) :: FileIn, FileOut, Line,Line2
  CHARACTER(32) :: Input,poscar

  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF
  INTEGER(KINT) :: IOpt
  INTEGER(KINT) :: Idx, LineL,LenLine
  INTEGER(KINT) :: N3at, NFree, Nfr
  ! Temporary integer counter
  INTEGER(KINT) :: NTmp

  INTEGER(KINT) :: List(2*MaxFroz)

  REAL(KREAL) :: E,FactStep,B(3) !E=Energy
  REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad
  REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord
  REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord)
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord
  REAL(KREAL), ALLOCATABLE :: HessTmp(:,:)
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
  LOGICAL, SAVE :: allocation_flag

  ! Flag to see if frozen atoms are frozen (see below)
  LOGICAL :: FChkFrozen

  ! Energies for all points along the path at the previous iteration
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF, where it is deallocated: Prakash
  ! Maximum step allowed for this geometry
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF, where it is deallocated: Prakash 

  LOGICAL, ALLOCATABLE :: FTmp(:) ! Nat
  LOGICAL :: FFrozen,FCart

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! To analyse the frozen atoms
  !
  REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat
  REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat
  REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist
  LOGICAL, ALLOCATABLE :: ListAtFroz(:)
  INTEGER(KINT) :: Iat1, Iat2, Iat3

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! For VASP
  !
  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
  INTEGER(KINT) :: NbType, NbTypeUser

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  LOGICAL :: TChk
  LOGICAL :: Debug, Fini,Flag_Opt_Geom

!  INTEGER(KINT), EXTERNAL :: Iargc
  INTEGER(KINT), EXTERNAL ::  ConvertNumAt

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! For Debugging purposes
  !
  LOGICAL :: FCalcHess
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


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

     SUBROUTINE Read_Geom(Input)
       CHARACTER(32), INTENT(IN) :: input
     END SUBROUTINE READ_GEOM

     subroutine hupdate_all (n, dx, dgrad, HessOld)

       IMPLICIT NONE

       INTEGER, PARAMETER :: KINT=KIND(1)
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)

       INTEGER(KINT), INTENT(IN) :: n
       real(KREAL) :: dx(n), dgrad(n), HessOld(n*n)
     END subroutine hupdate_all

     SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess)

       IMPLICIT NONE
       
       INTEGER, PARAMETER :: KINT=KIND(1)
       INTEGER, PARAMETER :: KREAL=KIND(1.D0)
       
       INTEGER(KINT), INTENT(IN) :: NFREE
       REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree)
       REAL(KREAL), INTENT(IN) :: ds(Nfree)
       REAL(KREAL), INTENT(IN) :: dGrad(Nfree)
       
     END subroutine Hinvup_BFGS_new


  END INTERFACE

  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
       MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, &
       SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, &
       input,prog, ProgExe,RunMode, AtTypes,poscar,PathName,  &
       OptGeom,IniHup, HupNeighbour, hesUpd, &
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
       NGintMax, Align

  Namelist/cartlist/list
  Namelist/frozenlist/list


  Flag_Opt_Geom = .FALSE. ! Initialized.

  ! We read the files names
!  SELECT CASE (iargc())
  SELECT CASE (command_argument_count())
  CASE (2)
     call getarg(1,FileIn)
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
     call getarg(2,FileOut)
     OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN')
  CASE (1)
     call getarg(1,FileIn)
     IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
        WRITE(*,*) "Path mini-help"
        WRITE(*,*) "--------------"
        WRITE(*,*) ""
        WRITE(*,*) "Use: Path Input_file Output_file"
        WRITE(*,*) "Input_file starts with a Namelist called path"
        WRITE(*,*) ""
        WRITE(*,*) "Compulsory variables are:"
        WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
        WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
        WRITE(*,*) "Nat   : Number of atoms"
        WRITE(*,*) ""
        WRITE(*,*) "Other options"
        WRITE(*,*) "Input: string that indicates the type of the input geometries."
        WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac or Ext"
        WRITE(*,*) "      In case of a Gaussian or Mopac calculations, input must be set to Cart."
        WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
        WRITE(*,*) "      input file.See example file Test_G03.path or Test_Mopac.path"
        WRITE(*,*) "      In the case of a VASP calculation, if input is set to Cart, then"
        WRITE(*,*) "    the preamble of a VASP calculation should be added at the end of"
        WRITE(*,*) "    the input file. See example file Test_VASP_cart.path"
        WRITE(*,*) "    In the case of a VASP calculation, one should also give value of the "
        WRITE(*,*) "    Runmode variable"
        WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy"
        WRITE(*,*) "       and gradient of the whole path or not. If one wants to use VASP,"
        WRITE(*,*) "       Runmode must be set to PARA."
        WRITE(*,*) "PathName: Prefix used to save the path. Default is Path"
        WRITE(*,*) "Poscar: string that will be used as the prefix for the different "
        WRITE(*,*) "        POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.,"
        WRITE(*,*) "        not used for internal calculations."
        WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations"
        WRITE(*,*) " ISuffix: Suffix for the input file."
        WRITE(*,*) " OSuffix: suffix for the output file."
        WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates."
        WRITE(*,*) "      Valid only for Coord=Zmat, Hybrid or Mixed"
        WRITE(*,*) "Fact: REAL used to define if two atoms are linked."
        WRITE(*,*) "      If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked."
        WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
        WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
        WRITE(*,*) "       - CART (or Xyz): works in cartesian"
        WRITE(*,*) "       - Zmat: works in internal coordinates (Zmat)"
        WRITE(*,*) "       - Mixed: frozen atoms, as well as atoms defined by the "
        WRITE(*,*) "       'cart' array(see below) are describe in CARTESIAN, whereas the"
        WRITE(*,*) "       others are described in Zmat" 
        WRITE(*,*) "       - Baker: use of Baker coordinates, also called delocalized internal coordinates"
        WRITE(*,*) "       - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
        WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
        WRITE(*,*) "       - RFO: Rational function optimization"
        WRITE(*,*) "       - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
        WRITE(*,*) " HesUpd: method to update the hessian. By default, it is Murtagh-Sargent"
        WRITE(*,*) "      Except for geometry optimization where it is BFGS."
        WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
        WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
        WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
        WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path."
        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
        WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
        WRITE(*,*) ""
        WRITE(*,*) "Arrays:"
        WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements."
        WRITE(*,*) "      You can modify it using, rcov(6)=0.8."
        WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements."
        WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations."
        WRITe(*,*) "If not given, Path will read the POTCAR file."
        WRITE(*,*) ""
        WRITE(*,*) "Flags:"
        WRITE(*,*) "MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
        WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time."
        WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
        WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
        WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
        WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
        WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
        WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path."
        WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
        WRITE(*,*) "         If True, a &frozenlist namelist containing the list of frozen atoms must be given."
        WRITE(*,*) "          If VASP is used, and frozen is not given"
        WRITE(*,*) " here, the program will use the F flags of the input geometry"
        WRITE(*,*) "FCart:  True if one wants to describe some atoms using cartesian coordinates. "
        WRITE(*,*) "        *** Only used in 'mixed' calculations. ***"
        WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
        WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
        WRITE(*,*) ""
        WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
        WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
        WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
        WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
        WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now"
        WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files."
        WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP."
        WRITE(*,*) ""
        STOP
     END IF
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
     IOOUT=6
  CASE (0)
     IOIN=5
     IOOUT=6
  END SELECT

!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 
  ! We print the Version info in file
  !
  WRITE(*,'(1X,A)') "Path v4.1793 (c) PFL/PD 2007-2010"


  ! We initiliaze variables
  Pi=dacos(-1.0d0)
  Nat=1
  Fact=1.1
  IGeomRef=-1
  NGeomI=1
  NGeomF=8
  IReparam=5
  IReparamT=-1
  ISpline=5
  debug=.False.
  MW=.TRUE.
  bohr=.False.
  Coord='ZMAT'
  Step_method='RFO'
  DebugFile='Path.valid'
  PathName="Path"
  MaxCyc=50
  SMax=0.1D0
  SThresh=-1.
  GThresh=-1.
  FTan=0.
  Hinv=.TRUE.
  IniHUp=.FALSE.
  HupNeighbour=.FALSE.
  HesUpd="MS"
  FFrozen=.FALSE.
  FCart=.FALSE.
  List=0
  renum=.TRUE.
  OptReac=.FALSE.
  OptProd=.FALSE.
  EReac=0.d0
  EProd=0.d0
  OptGeom=-1
  PathOnly=.FALSE.
  Prog='EMPTY'
  ProgExe='EMPTY'
  Runmode='SERIAL'
  Input='EMPTY'
  Poscar="POSCAR"
  AutoCart=.FALSE.
  VMD=.TRUE.
  WriteVASP=.FALSE.
  DynMaxStep=.FALSE.
  CalcName="EMPTY"
  ISuffix="EMPTY"
  OSuffix="EMPTY"
  NGintMax=10
  Align=.TRUE.

  ! Number of point used to compute the length of the path,
  ! and to reparameterize the path
  NMaxPtPath=-1 
  FrozTol=1e-4

  ! Read the namelist
  read (IOIN,path)

  debug=valid("Main")
  FCalcHess=valid("CalcHess")
  PathName=AdjustL(PathName)
  Coord=AdjustL(Coord)
  CALL UpCase(coord)

  Step_method=AdjustL(Step_method)
  CALL UpCase(Step_method)

  Prog=AdjustL(Prog)
  CALL UpCase(prog)

  Input=AdjustL(Input)
  CALL UpCase(input)

  Runmode=AdjustL(Runmode)
  Call UpCase(Runmode)
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"

  If (IReparamT<=0) IReparamT=IReparam

  IF (COORD.EQ.'CART') THEN
     Renum=.FALSE.
     IGeomRef=1
  END IF

  if (Prog.EQ.'EMPTY') THEN
     Prog='GAUSSIAN'
     If (ProgExe=="EMPTY") ProgExe="g09"
  END IF
  if (Prog.EQ.'G03') THEN
     Prog='GAUSSIAN'
     If (ProgExe=="EMPTY") ProgExe="g03"
  END IF
  if (Prog.EQ.'G09') THEN
     Prog='GAUSSIAN'
     If (ProgExe=="EMPTY") ProgExe="g09"
  END IF

  if (Prog.EQ.'TM') Prog="TURBOMOLE"
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"


  Select Case (Prog)
    CASE ("EXT")
     if (CalcName=="EMPTY") CalcName="Ext"
     if (ISuffix=="EMPTY")  ISuffix="in"
     if (OSuffix=="EMPTY")  OSuffix="out"
    CASE ('MOPAC')
     if (CalcName=="EMPTY") CalcName="Path"
     if (ISuffix=="EMPTY")  ISuffix="mop"
     if (OSuffix=="EMPTY")  OSuffix="out"
    CASE ("GAUSSIAN")
     if (CalcName=="EMPTY") CalcName="Path"
     if (ISuffix=="EMPTY")  ISuffix="com"
     if (OSuffix=="EMPTY")  OSuffix="log"
    CASE DEFAULT
     if (CalcName=="EMPTY") CalcName="Path"
     if (ISuffix=="EMPTY")  ISuffix="in"
     if (OSuffix=="EMPTY")  OSuffix="out"
  END SELECT 


  if (Input.EQ.'EMPTY') THEN
     Select Case (Prog)
       CASE ("VASP")
          Input=Prog
       CASE ("CHAMFRE")
          Input="VASP"
       CASE DEFAULT
          Input='XYZ'
     END SELECT
    WRITE(*,*) "Input has been set to the default: ",INPUT
 END IF

  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))

  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
     STOP
  END IF

  if (OptGeom.GE.1) THEN
     HesUpd="BFGS"
     IF (IGeomRef.LE.0) THEN
        IGeomRef=OptGeom
     ELSE
        IF (IGeomRef.NE.OptGeom) THEN
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom

           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
           STOP
        END IF
     END IF
  END IF

  IF (ProgExe.EQ.'EMPTY') THEN
     SELECT CASE (PROG)
     CASE ("GAUSSIAN")
        ProgExe="g03"
     CASE ("MOPAC")
        ProgExe="mopac"
     CASE ("EXT")
        ProgExe="./Ext.exe"
     CASE ("VASP")
        ProgExe='VASP'
     CASE ("TURBOMOLE")
        ProgExe='jobex -c 1 -keep'
     CASE ("TEST")
        ProgExe=""
     CASE ("CHAMFRE")
        ProgExe=""
     CASE DEFAULT
        WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
        STOP
     END SELECT
  END IF

  if (FCart.AND.FFrozen) THEN
     ! We have to be carreful because user might use any order
     ! to give the two lists CartList and FrozenList
     READ(IOIN,'(A)') Line
     Call Upcase(Line)
     if (Index(Line,"CART").NE.0) THEN
        ALLOCATE(Cart(Nat))
        BackSpace(IOIN)
        List=0
        READ(IOIN,cartlist)
        IF (COORD.NE.'CART') THEN
           Cart=List(1:Nat)
           if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
        ELSE
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
           WRITE(*,*) "I will discard the cartlist"
           Cart=0
        END IF

        ALLOCATE(Frozen(Nat))
        BACKSPACE(IOIN)
        List=0
        READ(IOIN,Frozenlist)
        Frozen=List(1:Nat)
        if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
     ELSE
        ALLOCATE(Frozen(Nat))
        BACKSPACE(IOIN)
        List=0
        READ(IOIN,Frozenlist)
        Frozen=List(1:Nat)
        if (debug) WRITE(*,*) "DBG Path L523 Frozen",Frozen
        ALLOCATE(Cart(Nat))
        BackSpace(IOIN)
        List=0
        READ(IOIN,cartlist)
        IF (COORD.NE.'CART') THEN
           Cart=List(1:Nat)
           if (debug) WRITE(*,*) "DBG Path L548 Cart",Cart
        ELSE
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
           WRITE(*,*) "I will discard the cartlist"
           Cart=0
        END IF
     END IF
  ELSE
     if (FCart) THEN
        ALLOCATE(Cart(Nat))
        List=0
        READ(IOIN,cartlist)
        IF (COORD.NE.'CART') THEN
           Cart=List(1:Nat)
           if (debug) WRITE(*,*) "DBG Path L561 Cart",Cart
        ELSE
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
           WRITE(*,*) "I will discard the cartlist"
           Cart=0
        END IF
     ELSE
        IF ((PROG=="VASP").OR.(AutoCart)) THEN
           ALLOCATE(Cart(Nat))
        ELSE
           ALLOCATE(Cart(1))
        END IF
        Cart=0
     END IF

     if (FFrozen) THEN
        ALLOCATE(Frozen(Nat))
        List=0
        READ(IOIN,Frozenlist)
        Frozen=List(1:Nat)
        WRITE(*,*) Frozen
        if (debug) WRITE(*,*) "DBG Path L549 Frozen",Frozen
     ELSE
        IF (PROG=="VASP") THEN
           ALLOCATE(Frozen(Nat))
        ELSE
           ALLOCATE(Frozen(1))
        END IF
        Frozen=0
     END IF
  END IF


  If (FCart) Call Expand(Cart,NCart,Nat)
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)

  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen

  if (NMaxPtPath.EQ.-1) then
     NMaxPtPath=min(50*NGeomF,1000)
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
  end if

  !  NMaxPtPath=10000

  ! We ensure that FTan is in [0.:1.]
  FTan=min(abs(FTan),1.d0)

  if (debug) WRITE(IOOUT,path)

  FTan=1.-FTan

  !Thresholds...
  if (SThresh.LE.0) SThresh=0.01d0
  If (GThresh.LE.0) GThresh=SThresh/4.
  SThreshRms=Sthresh*2./3.
  GThreshRms=Gthresh*2./3.

  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
  ALLOCATE(MassAt(NAt),OrderInv(Nat))

  ! We read the initial geometries
  Call Read_Geom(input)

  ! We convert atome names into atom mass number
  DO I=1,NAt
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
     Atome(I)=ConvertNumAt(AtName(I))
  END DO


  ! PFL 23 Sep 2008 ->
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
  ! between the first geometry and the others.
  IF (NFroz.GT.0) THEN
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
     ListAtFroz=.FALSE.

     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)

     SELECT CASE (NFroz)
        ! We have Frozen atoms
        ! PFL 28 Oct 2008 ->
        ! It might happen that the geometries are translated/rotated.
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
        ! three frozen atoms.

     CASE (3:)
        DO I=1,NFroz
           ListAtFroz(Frozen(I))=.TRUE.
        END DO
     CASE (2)
        IAt1=Frozen(1)
        IAt2=Frozen(2)
        ListAtFroz(Iat1)=.TRUE.
        ListAtFroz(Iat2)=.TRUE.
        x2(1)=x1(Iat2)-x1(Iat1)
        y2(1)=y1(Iat2)-y1(Iat1)
        z2(1)=z1(Iat2)-z1(Iat1)
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
        Dist=-1.
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
        DO I=1,Nat
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
           x2(2)=x1(Iat2)-x1(I)
           y2(2)=y1(Iat2)-y1(I)
           z2(2)=z1(Iat2)-z1(I)
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
           if (abs(ca)<0.9) THEN
              IF (Norm2>Dist) THEN
                 Iat3=I
                 Dist=Norm2
              END IF
           END IF
        END DO
        ListAtFroz(Iat3)=.TRUE.
        if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3
     CASE (1)
        IAt1=Frozen(1)
        ListAtFroz(Iat1)=.TRUE.
        Dist=-1.
        ! We scan all atoms to find one that is further from At1
        DO I=1,Nat
           if (I.EQ.Iat1) Cycle
           x2(1)=x1(Iat1)-x1(I)
           y2(1)=y1(Iat1)-y1(I)
           z2(1)=z1(Iat1)-z1(I)
           Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
           IF (Norm1>Dist) THEN
              Iat2=I
              Dist=Norm1
           END IF
        END DO
        IAt2=Frozen(2)
        ListAtFroz(Iat2)=.TRUE.
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
        x2(1)=x1(Iat2)-x1(Iat1)
        y2(1)=y1(Iat2)-y1(Iat1)
        z2(1)=z1(Iat2)-z1(Iat1)
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
        Dist=-1.
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
        DO I=1,Nat
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
           x2(2)=x1(Iat2)-x1(I)
           y2(2)=y1(Iat2)-y1(I)
           z2(2)=z1(Iat2)-z1(I)
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
           if (abs(ca)<0.9) THEN
              IF (Norm2>Dist) THEN
                 Iat3=I
                 Dist=Norm2
              END IF
           END IF
        END DO
        ListAtFroz(Iat3)=.TRUE.
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3
     END SELECT

     DO I=2,NGeomI
        ! First, we align the Ith geometry on the first one, using only
        ! the positions of the "frozen" atoms. (see above: it might be that
        ! ListAtFroz contains non frozen atoms usd to align the geometries
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)


        FChkFrozen=.FALSE.
        DO J=1,NFroz
           IAt=Frozen(J)
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
              FChkFrozen=.TRUE.
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
              WRITE(*,*) "from geometry 1 to geometry ",I
           END IF
        END DO
     END DO
     IF (FChkFrozen) THEN
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
        STOP
     END IF
  END IF


  N3at=3*Nat
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.

  ! if prog=gaussian, we have to read a full input
  ! to mimic it later
  if (prog=='GAUSSIAN') THEN
     ! We read the Gaussian input file
     ! First, the root
     IF (DEBUG) WRITE(*,*) "Reading Gauss Root"
     ALLOCATE(Gauss_Root)
     NULLIFY(Gauss_Root%next)
     Current => Gauss_root
     LineL=1
     DO WHILE (LineL.NE.0)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        Idx=INDEX(Line,"chk")
        IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN
           current%Line=TRIM(Line)
           ALLOCATE(current%next)
           Current => Current%next
           Nullify(Current%next)
        END IF
     END DO

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

     ! Now the comment... 
     IF (DEBUG) WRITE(*,*) "Reading Gauss Comment"
     ALLOCATE(Gauss_Comment)
     NuLLIFY(Gauss_Comment%Next)
     Current => Gauss_comment
     LineL=1
     DO WHILE (LineL.NE.0)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        IF (LineL.NE.0) THEN
           current%Line=TRIM(Line)
           ALLOCATE(current%next)
           Current => Current%next
           Nullify(Current%next)
        END IF
     END DO

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

     ! Now the charge
     IF (DEBUG) WRITE(*,*) "Reading Gauss Charge"
     READ(IOIN,'(A)') Gauss_Charge
!     WRITE(*,*) Gauss_charge
     ! We now read the Paste part...
     ALLOCATE(Gauss_Paste(NAt))
     LenLine=1
     Iat=0
     Gauss_paste=" "
     DO While (LenLine.GT.0)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LenLine=Len_TRIM(Line)
        IF (LenLine.GT.0) Iat=Iat+1
        Idx=Index(Line,'.',BACK=.TRUE.)
        !          WRITE(*,*) Idx,TRIM(Line)
        !          WRITE(*,*) Idx,TRIM(Line(Idx+1:))
        Line=ADJUSTL(Line(Idx+1:))
        Idx=Index(Line," ")
        LineL=LEN_TRIM(Line(Idx:))
        !          WRITE(*,*) LineL,TRIM(Line(Idx:))
        IF (LineL.GT.0) THEN 
           Gauss_paste(Iat)=Line(Idx:)
        END IF
     END DO
     IF (Iat.NE.Nat) THEN
        WRITE(*,*) "I found ", Iat," lines for the geometry instead of ",Nat
        WRITE(*,*) "ERROR. STOP"
        WRITE(IOOUT,*) "I found ", Iat," lines for the geometry instead of ",Nat
        WRITE(IOOUT,*) "ERROR. STOP"
        STOP
     END IF
     ! We now read the last part
     IF (DEBUG) WRITE(*,*) "Reading Gauss End"
     !     READ(IOIN,'(A)') Line
     ALLOCATE(Gauss_End)
     NuLLIFY(Gauss_End%Next)
     Current => Gauss_End
     LineL=1
     DO WHILE (1.EQ.1)
        READ(IOIN,'(A)',END=999) Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        current%Line=TRIM(Line)
        ALLOCATE(current%next)
        Current => Current%next
        Nullify(Current%next)
     END DO
999  CONTINUE

     !     ! Write the gaussian input file for testing purposes
     !     Current => Gauss_root
     !     DO WHILE (ASSOCIATED(Current%next))
     !        WRITE(*,'(1X,A)') Trim(current%line)
     !        Current => current%next
     !     END DO

     !        WRITE(*,*) 
     !        WRITE(*,*) 'Comment original'

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

     !        WRITE(*,*) 
     !        WRITE(*,*) Trim(Gauss_charge)

     !        DO I=1,Nat
     !           WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XyzGeomI(1:3,I,1), TRIM(Gauss_Paste(I))
     !        END DO

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

     !        WRITE(*,*) 

  END IF

! if prog=mopac, we have to read a full input to reproduce it later
! The structure is:
! A MOPAC data set normally consists of one line of keywords, two lines of user-defined text, then the coordinates
! Then  a blank line or a line of 0.
! then the symmetry description.
! comment lines start with * and can be anywhere !!!

  if (prog=='MOPAC') THEN
     ! We read the Gaussian input file
     ! First, the root
     IF (DEBUG) WRITE(*,*) "Reading Mopac input"
     ALLOCATE(Mopac_Root)
     NULLIFY(Mopac_Root%next)
     ALLOCATE(Mopac_Comment)
     NULLIFY(Mopac_Comment%next)
     ALLOCATE(Mopac_End)
     NuLLIFY(Mopac_End%Next)
     Current => Mopac_root
     CurCom => Mopac_Comment
     LineL=1
     NTmp=0
     DO WHILE (NTmp.LT.3)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        IF (Line(1:1)/="*") THEN
           IF (NTmp==0) THEN
              Line2=Line
              Call UpCase(Line2)
              Idx=Index(Line2,'GRADIENTS')
              If (Idx==0) Line=TRIM(Line) // " GRADIENTS"
              Idx=Index(Line2,'1SCF')
              If (Idx==0) Line=TRIM(Line) // " 1SCF"
           END IF
           current%Line=TRIM(Line)
           ALLOCATE(current%next)
           Current => Current%next
           Nullify(Current%next)
           NTmp=NTmp+1
        ELSE
           CurCom%Line=TRIM(LINE)
           ALLOCATE(CurCom%Next)
           CurCom => CurCom%Next
           NULLIFY(CurCom%Next)
        END IF
     END DO

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

! Now the geometry... that we just skip :)
     IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry"
     Mopac_EndGeom=""
     LineL=1
     DO WHILE (LineL.NE.0)
        READ(IOIN,'(A)',END=989) Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
! The last line might be either blank or filled with 0
        IF (Line(1:1)=="0") THEN
           LineL=0
           Mopac_EndGeom=Trim(Line)
        END IF
        IF (Line(1:1)=="*") THEN
           CurCom%Line=TRIM(LINE)
           ALLOCATE(CurCom%Next)
           CurCom => CurCom%Next
           NULLIFY(CurCom%Next)
        END IF
     END DO

! If we are here, there might be something else to read: Mopac_end

     ! We now read the last part
     IF (DEBUG) WRITE(*,*) "Reading Mopac End"
     !     READ(IOIN,'(A)') Line
     Current => Mopac_End
     LineL=1
     DO WHILE (1.EQ.1)
        READ(IOIN,'(A)',END=989) Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        IF (Line(1:1)/="*") THEN
           current%Line=TRIM(Line)
           ALLOCATE(current%next)
           Current => Current%next
           Nullify(Current%next)
           NTmp=NTmp+1
        ELSE
           CurCom%Line=TRIM(LINE)
           ALLOCATE(CurCom%Next)
           CurCom => CurCom%Next
           NULLIFY(CurCom%Next)
        END IF
     END DO
989  CONTINUE

  END IF ! IF (PROG="MOPAC")

  if ((Prog=="VASP").AND.(input/="VASP")) THEN
     ! Input was not Vasp, so many parameters are missing like lattice 
     ! constants...
     ! we read them now !
     ALLOCATE(FFF(3,nat))
     ! First geometry is a bit special for VASP as we have to set
     ! many things
     IF (DEBUG) WRITE(*,*) "Reading Vasp Parameters"
     READ(IOIN,'(A)') Vasp_Title
     READ(IOIN,*) Vasp_param

     READ(IOIN,*) Lat_a
     READ(IOIN,*) Lat_b
     READ(IOIN,*) Lat_c

     Lat_a=Lat_a*Vasp_param
     Lat_b=Lat_b*Vasp_param
     Lat_c=Lat_c*Vasp_param

     ALLOCATE(NbAtType(nat))
     READ(IOIN,'(A)') Vasp_types
     ! First, how many different types ?
     NbAtType=0
     READ(Vasp_types,*,END=998) NbAtType
998  CONTINUE
     NbType=0
     DO WHILE (NbAtType(NbType+1).NE.0)
        NbType=NbType+1
     END DO

     ! Do we know the atom types ?
     IF (AtTypes(1).EQ.'  ') THEN
        ! user has not provided atom types... we have to find them ourselves
        ! by looking into the POTCAR file...
        INQUIRE(File="POTCAR", EXIST=Tchk)
        IF (.NOT.Tchk) THEN
           WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found"
           STOP
        END IF
        OPEN(IOTMP,File="POTCAR")
        DO I=1,NbType
           Line='Empty'
           DO WHILE (Line(1:2).NE.'US')
              READ(IOTMP,'(A)') Line
              Line=AdjustL(Line)
           END DO
           Line=adjustl(Line(3:))
           AtTypes(I)=Line(1:2)
        END DO
        if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType)
        CLOSE(IOTMP)

     ELSE  !AtTypes(1).EQ.'  '
        ! user has provided atom types
        NbTypeUser=0
        DO WHILE (AtTypes(NbTypeUser+1).NE.'  ')
           NbTypeUser=NbTypeUser+1
        END DO
        IF (NbType.NE.NbTypeUser) THEN
           WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes"
           STOP
        END IF
     END IF

     IAt=1
     DO I=1,NbType
        DO J=1,NbAtType(I)
           AtName(Iat)=AtTypes(I)
           Iat=Iat+1
        END DO
     END DO
     DEALLOCATE(NbAtType)

     NbTypes=NbType

     READ(IOIN,'(A)') Vasp_comment
     READ(IOIN,'(A)') Vasp_direct
     V_direct=Adjustl(Vasp_direct)
     Call UpCase(V_direct)

  END IF

  ! In the case of VASP there is always the problem  of moving from one side
  ! of the box to the other...
! PFL 2010 Dec 2 ->
! Here we deal with input and not prog in fact 
!  if (prog.EQ.'VASP') THEN
  if (Input.EQ.'VASP') THEN
! <- PFL 2010 Dec 2
     Renum=.TRUE.

     ! V_direct has been set in Read_geom
     IF (V_direct(1:6).EQ.'DIRECT') THEN
        Latr(1:3,1)=Lat_a
        Latr(1:3,2)=Lat_b
        Latr(1:3,3)=Lat_c
        B=1.
        CALL Gaussj(Latr,3,3,B,1,1)
     ELSE
        Latr=0.
        Latr(1,1)=1.d0
        Latr(2,2)=1.d0
        Latr(3,3)=1.d0
     END IF

     ! Actualization of Frozen using the FFFF... 
     ! Frozen has the advantage ie if given, it imposes _ALL_ the FFF flags.
     IF (Frozen(1).NE.0) THEN
        WRITE(IOOUT,*) "Frozen is given. Flags of the given POSCAR are overriden"
        FFF='T'

        NFroz=0
        DO WHILE (Frozen(NFroz+1).NE.0)
           NFroz=NFroz+1
           FFF(1:3,Frozen(NFroz))='F'
        END DO
     ELSE
        WRITE(IOOUT,*) "Frozen not given : using  Flags of the given POSCAR"
        NFroz=0
        Frozen=0
        DO I=1,Nat
           IF ((FFF(1,I).EQ.'F').OR.(FFF(2,I).EQ.'F').OR.(FFF(3,I).EQ.'F')) THEN
              FFF(1:3,I)='F'
              NFroz=NFroz+1
              Frozen(NFroz)=I
           END IF
        END DO
        WRITE(IOOUT,*) "Frozen atoms:",Frozen(1:NFroz)
     END IF


  END IF

  IF (Prog=="VASP") THEN
     IF (Vmd) THEN
        if (debug) WRITE(*,*) "DBG MAIN, L803, VMD=T,NbTypes,AtTypes",NbTypes,AtTypes(1:NbTypes)
        Line=""
        DO I=1,NbTypes
           Line=TRIM(Line) // " " // TRIM(AdjustL(AtTypes(I))) 
        END DO
        Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title))
        if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title)
     END IF
     Call CheckPeriodicBound
  END IF

  IF (COORD=="MIXED") Call TestCart(AutoCart)

  SELECT CASE(NFroz)
  CASE (2)
     IntFroz=1
  CASE (3)
     IntFroz=3
  CASE (4:)
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
  END SELECT

  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart

  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
     COORD="ZMAT"
  END IF

  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz

  SELECT CASE(COORD)
  CASE ('ZMAT')
     NCoord=Nfree
     ALLOCATE(dzdc(3,nat,3,nat))
  CASE ('MIXED')
     SELECT CASE (NCart+NFroz)
     CASE (1) 
        NCoord=N3At-3
     CASE (2)
        NCoord=N3At-1
     CASE (3:)
        NCoord=N3At
     END SELECT
     ALLOCATE(dzdc(3,nat,3,nat))
  CASE ('BAKER')
     Symmetry_elimination=0
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
     ALLOCATE(BTransInv(NCoord,3*Nat))
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
  CASE ('HYBRID')
     NCoord=N3at
  CASE ('CART')
     NCoord=N3at
  END SELECT

  if (debug) WRITE(*,*) "DBG Path, L826, Coord, NCart,NCoord, N3At",TRIM(COORD),Ncart, NCoord,N3at

  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))

  ZeroVec=0.d0
  DO I =1, NGeomF
     IntTangent(I,:)=0.d0
  END DO

  if (debug) THEN
     Print *, 'L886, IntTangent(1,:)='
     Print *, IntTangent(1,:)
  END IF

  if (.NOT.OptReac) Energies(1)=Ereac
  if (.NOT.OptProd) Energies(NGeomF)=EProd
  MaxStep=SMax

  ! We initialize the mass array
  if (MW) THEN
     WRITE(*,*) 'Working in MW coordinates'
     DO I=1,Nat
        massAt(I)=Mass(Atome(I))
     END DO
  ELSE
     DO I=1,Nat
        MassAt(I)=1.d0
     END DO
  END IF

  WRITE(*,*) "Prog=",TRIM(PROG)

  ALLOCATE(FrozAtoms(Nat))

  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
  IF (NFroz.EQ.0) THEN
     FrozAtoms=.TRUE.
  ELSE
     I=1
     NFr=0
     FrozAtoms=.False.
     DO WHILE (Frozen(I).NE.0)
        IF (Frozen(I).LT.0) THEN
           DO J=Frozen(I-1),abs(Frozen(I))
              IF (.NOT.FrozAtoms(J)) THEN
                 NFr=NFr+1
                 FrozAtoms(J)=.TRUE.
              END IF
           END DO
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
           FrozAtoms(Frozen(I))=.TRUE.
           NFr=NFr+1
        END IF
        I=I+1
     END DO
     IF (NFr.NE.NFroz) THEN
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
        STOP
     END IF
  END IF

  if (debug) THEN
     !Print *, 'L968, IntTangent(1,:)='
     !Print *, IntTangent(1,:)
  END IF

  ! We have read everything we needed... time to work :-)
  IOpt=0
  FirstTimePathCreate = .TRUE. ! Initialized.
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
  FirstTimePathCreate = .FALSE.

  if (debug) THEN
     Print *, 'L980, IntTangent(1,:)='
     Print *, IntTangent(1,:)
  END IF

  ! PFL 30 Mar 2008 ->
  ! The following is now allocated in Calc_Baker. This is more logical to me
  !   IF (COORD .EQ. "BAKER") Then
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
  !   END IF
  ! <- PFL 30 mar 2008

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! For Debugging purposes
  !
  ! OptGeom is the index of the geometry which is to be optimized.
  IF (Optgeom.GT.0) THEN
     Flag_Opt_Geom=.TRUE.
     SELECT CASE(Coord)
     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(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
!!! <- CAUTION : PFL 29.AUG.2008
     CASE ('ZMAT','MIXED')
        !GeomTmp=IntCoordF(OptGeom,:)
        GeomTmp=IntCoordI(OptGeom,:)
     CASE ('BAKER')
        !GeomTmp=IntCoordF(OptGeom,:)
        GeomTmp=IntCoordI(OptGeom,:)
     CASE DEFAULT
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
        STOP
     END SELECT

     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
     Flag_Opt_Geom = .TRUE.
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)

     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
     WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
     STOP
  END IF ! IF (Optgeom.GT.0) THEN

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! End of GeomOpt
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  IF (PathOnly) THEN
     WRITE(*,*) "PathOnly=.T. , Stop here"
     Call Write_path(-1)
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
     STOP
  END IF

  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)

  DEALLOCATE(XyzGeomI,IntCoordI)
  NGeomI=NGeomF
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))

  IF (DEBUG) WRITE(*,*) ' Back from PathCreate, L965, Path.f90'
  ! we print this path, which is the first one :-) called Iteration0
  ! we now have to calculate the gradient for each point (and the energy
  ! if we work at 0K)

  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
     DO IGeom=1,NGeomF
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
        DO I=1,3*Nat-6
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
        END DO
     END DO
  END IF

  ! To calculate B^T^-1 for all extrapolated final geometries:
  ! PFL 18 Jan 2008 ->
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
  ! <-- PFL 18 Jan 2008

  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
     DO IGeom=1,NGeomF
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
        DO I=1,3*Nat-6
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
        END DO
     END DO
  END IF

  IOpt=0
  Call EgradPath(IOpt,Flag_Opt_Geom)

  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
  Call Write_path(IOpt)
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)

  ! We have the geometries, the first gradients... we will generate the first hessian matrices :-)
  ! ... v3.71 Only if IniHUp=.TRUE. which is not the default. In this version, we take the hessian
  ! of end points as Unity matrices, and we update them along the path.

  ! By default, Hess=Id
  Hess=0.
  IF(Coord .EQ. "ZMAT") Then
     ! We use the fact that we know that approximate good hessian values
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
     DO IGeom=1,NGeomF
        Hess(IGeom,1,1)=0.5d0
        Hess(IGeom,2,2)=0.5d0
        Hess(IGeom,3,3)=0.1d0
        DO J=1,Nat-3
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
        END DO
        IF (HInv) THEN
           DO I=1,NCoord
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
           END DO
        END IF
     END DO
  ELSE
     DO I=1,NCoord
        DO J=1,NGeomF
           Hess(J,I,I)=1.d0
        END DO
     END DO
  END IF

  IF (COORD .EQ. "BAKER") THEN
     ! UMat(NPrim,NCoord)
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
     Hprim=0.d0
     ScanCoord=>Coordinate
     I=0
     DO WHILE (Associated(ScanCoord%next))
        I=I+1
        SELECT CASE (ScanCoord%Type)
        CASE ('BOND')
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
           Hprim(I,I) = 0.5d0
           !END DO
        CASE ('ANGLE')
           Hprim(I,I) = 0.2d0
        CASE ('DIHEDRAL')
           Hprim(I,I) = 0.1d0
        END SELECT
        ScanCoord => ScanCoord%next
     END DO

     ! HprimU = H^prim U:
     HprimU=0.d0
     DO I=1, 3*Nat-6
        DO J=1,NPrim
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
        END DO
     END DO

     ! Hess = U^T HprimU = U^T H^prim U:
     Hess=0.d0
     DO I=1, 3*Nat-6
        DO J=1,NPrim
           ! UMat^T is needed below.
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
        END DO
     END DO
     DO K=2,NGeomF
        Hess(K,:,:)=Hess(1,:,:)
     END DO
     DEALLOCATE(Hprim,HprimU)
  end if !  ends IF COORD=BAKER
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
  ALLOCATE(HessTmp(NCoord,NCoord))

  IF (IniHUp) THEN
     IF (FCalcHess) THEN
        ! The numerical calculation of the Hessian has been switched on
        DO IGeom=1,NGeomF
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              GeomTmp=IntCoordF(IGeom,:)
           ELSE
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
           END IF
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
        END DO
     ELSE
        ! We generate a forward hessian and a backward one and we mix them.
        ! First, the forward way:
        DO IGeom=2,NGeomF-1
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
           ELSE
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
           END IF
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
           HessTmp=Hess(IGeom-1,:,:)
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
           Hess(IGeom,:,:)=HessTmp
        END DO

        ! Now, the backward way:
        HessTmp=Hess(NGeomF,:,:)
        DO IGeom=NGeomF-1,2,-1
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
           ELSE
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
           END IF
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
           ca=cos(alpha)
           sa=sin(alpha)
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
        END DO
     END IF ! matches IF (FCalcHess)
  END IF ! matches IF (IniHUp) THEN

  ! For debug purposes, we diagonalize the Hessian matrices
  IF (debug) THEN
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
     DO I=1,NGeomF
        WRITE(*,*) "DBG Main - Point ",I
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
        DO J=1,NCoord
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
        END DO
     END DO
  END IF ! matches IF (debug) THEN

  ! we have the hessian matrices and the gradients, we can start the optimizations !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Main LOOP : Optimization of the Path
  !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (debug)   Print *, 'NGeomF=', NGeomF
  allocation_flag = .TRUE.

  Fini=.FALSE.

  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
     IOpt=IOpt+1

     SELECT CASE (COORD)
     CASE ('ZMAT','MIXED','BAKER')
        GeomOld_all=IntCoordF
     CASE ('CART','HYBRID')
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
     CASE DEFAULT
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
        STOP
     END SELECT

     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
          valid("printtangent")) THEN
        WRITE(*,*) "Visualization of tangents"
        Idx=min(12,NCoord)
        DO I=1,NGeomF
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
           WRITE(*,*) 
        END DO
        WRITE(*,*) "END of tangents"
     END IF

     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
          valid("printgrad")) THEN
        WRITE(*,*) "Visualization of gradients"
        Idx=min(12,NCoord)
        DO I=1,NGeomF
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
           WRITe(*,*) 
        END DO
        WRITE(*,*) "END of gradients"
     END IF

     Fini=.TRUE.
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
        IGeom=1
        if (debug) WRITE(*,*) '**************************************'
        if (debug) WRITE(*,*) 'Optimizing reactant'
        if (debug) WRITE(*,*) '**************************************'
        SELECT CASE (COORD)
        CASE ('ZMAT','MIXED','BAKER')
           SELECT CASE (Step_method)
           CASE ('RFO')
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
              Print *, GradTmp
           CASE ('GDIIS')
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
              Print *, GradTmp
           CASE ('GEDIIS')
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
              ! Energies are updated in EgradPath.f90
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
                   NCoord,allocation_flag,ZeroVec)
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
              !GradTmp = GeomTmp - IntCoordF(1,:)
              !Print *, 'Old Geometry:'
              !Print *, IntCoordF(1,:)
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
              Print *, GradTmp
              !Print *, 'GeomTmp:'
              !Print *, GeomTmp
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
           CASE DEFAULT
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
              STOP
           END SELECT
        CASE ('HYBRID','CART')
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
           STOP
        END SELECT

        IF (debug) THEN
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
        END IF

        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
        Fini=(Fini.AND.(NormStep.LE.SThresh))
        OptReac=(NormStep.GT.SThresh)
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
        !Print *, 'Grad(1,:):'
        !Print *, Grad(1,:)
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac

        GradTmp=GradTmp*FactStep

        ! we take care that frozen atoms do not move.
        IF (NFroz.NE.0) THEN
           SELECT CASE (COORD)
           CASE ('ZMAT','MIXED')
              GradTmp(1:IntFroz)=0.d0
           CASE ('CART','HYBRID')
              DO I=1,Nat
                 IF (FrozAtoms(I)) THEN
                    Iat=Order(I)
                    GradTmp(3*Iat-2:3*Iat)=0.d0
                 END IF
              END DO
           CASE('BAKER')
              GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..)
           CASE DEFAULT
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP'
              STOP
           END SELECT
        END IF ! matches IF (NFroz.NE.0) THEN

        IF (debug) THEN
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
        END IF

        IF (debug) THEN
           !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1'
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
           ELSE
              DO Iat=1,Nat
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
                      XyzGeomF(IGeom,1:3,Iat)
              END DO
           END IF
        END IF ! matches IF (debug) THEN

        SELECT CASE (COORD)
        CASE ('ZMAT','MIXED','BAKER')
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
        CASE ('HYBRID','CART')
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
           STOP
        END SELECT

        IF (debug) THEN
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
           ELSE
              DO Iat=1,Nat

!!!!!!!!!! There was an OrderInv here... should I put it back ?
                 ! It was with indzmat. IndZmat cannot appear here...
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
                      XyzGeomF(IGeom,1:3,Iat)
              END DO
           END IF
        END IF ! matches IF (debug) THEN

        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
     ELSE ! Optreac
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
     END IF ! matches IF (OptReac) THEN

     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
        IGeom=NGeomF
        if (debug) WRITE(*,*) '******************'
        if (debug) WRITE(*,*) 'Optimizing product'
        if (debug) WRITE(*,*) '******************'
        SELECT CASE (COORD)
        CASE ('ZMAT','MIXED','BAKER')
           Print *, 'Optimizing product'
           SELECT CASE (Step_method)
           CASE ('RFO')
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
              Print *, GradTmp
           CASE ('GDIIS')
              ! GradTmp becomes "step" below and has INTENT(OUT):
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
              Print *, GradTmp
           CASE ('GEDIIS')
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
              Print *, GradTmp
           CASE DEFAULT
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
              STOP
           END SELECT
        CASE ('HYBRID','CART')
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
           STOP
        END SELECT

        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
        Fini=(Fini.AND.(NormStep.LE.SThresh))
        OptProd=.NOT.(NormStep.LE.SThresh)
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd

        GradTmp=GradTmp*FactStep

        ! we take care that frozen atoms do not move
        IF (NFroz.NE.0) THEN
           SELECT CASE (COORD)
           CASE ('ZMAT','MIXED','BAKER')
              GradTmp(1:IntFroz)=0.d0
           CASE ('CART','HYBRID')
              DO I=1,Nat
                 IF (FrozAtoms(I)) THEN
                    Iat=Order(I)
                    GradTmp(3*Iat-2:3*Iat)=0.d0
                 END IF
              END DO
           CASE DEFAULT
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
              STOP
           END SELECT
        END IF ! matches IF (NFroz.NE.0) THEN

        IF (debug) THEN
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
           ELSE
              DO Iat=1,Nat
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
                      XyzGeomF(IGeom,1:3,Iat)
              END DO
           END IF
        END IF

        SELECT CASE (COORD)
        CASE ('ZMAT','MIXED','BAKER')
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
        CASE ('HYBRID','CART')
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
           STOP
        END SELECT

        IF (debug) THEN
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
           ELSE
              DO Iat=1,Nat
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
                      XyzGeomF(IGeom,1:3,Iat)
              END DO
           END IF
        END IF
     ELSE ! Optprod
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
     END IF !matches IF (OptProd) THEN 

     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !
     !  Optimizing other geometries
     !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



     DO IGeom=2,NGeomF-1 ! matches at L1556
        if (debug) WRITE(*,*) '****************************'
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
        if (debug) WRITE(*,*) '****************************'

        SELECT CASE (COORD)
        CASE ('ZMAT','BAKER','MIXED')
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
           ! PFL 6 Apr 2008 ->
           ! Special case, if FTan=0. we do not consider Tangent at all.
           ! To DO: add the full treatment of FTan
           if (FTan==0.) GradTmp2=ZeroVec
           ! <- PFL 6 Apr 2008
           if (debug) THEN
              Print *, 'L1500, IntTangent(',IGeom,',:)='
              Print *, IntTangent(IGeom,:)
           END IF
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
           SELECT CASE (Step_method)
           CASE ('RFO')
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
           CASE ('GDIIS')
              !The following has no effect at IOpt==1
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&			  
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
              Print *, GradTmp
           CASE ('GEDIIS')
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
              Print *, GradTmp
           CASE DEFAULT
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
              STOP
           END SELECT
        CASE ('HYBRID','CART')
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
           ! so we change it
           DO I=1,Nat
              DO J=1,3
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
              END DO
           END DO
           !           GradTmp2=XyzTangent(IGeom,:)
           ! PFL 6 Apr 2008 ->
           ! Special case, if FTan=0. we do not consider Tangent at all.
           ! To DO: add the full treatment of FTan
           if (FTan==0.) GradTmp2=ZeroVec
           ! <- PFL 6 Apr 2008

           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
           STOP
        END SELECT

        ! we take care that frozen atoms do not move
        IF (NFroz.NE.0) THEN
           SELECT CASE (COORD)
           CASE ('ZMAT','MIXED','BAKER')
              IF (debug) THEN 
                 WRITE(*,*) 'Step computed. Coord=',Coord
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
              END IF
              GradTmp(1:IntFroz)=0.d0
           CASE ('CART','HYBRID')
              if (debug) THEN
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
                 DO I=1,Nat
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
                 END DO
              END IF
              DO I=1,Nat
                 IF (FrozAtoms(I)) THEN
                    if (debug) THEN
                       write(*,*) 'Step Computed. Coord=',Coord
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
                    END IF
                    Iat=Order(I)
                    GradTmp(3*Iat-2:3*Iat)=0.d0
                 END IF
              END DO

                 if (debug) THEN
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
                    DO I=1,Nat
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
                    END DO
                 END IF

           CASE DEFAULT
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
              STOP
           END SELECT
        END IF !matches IF (NFroz.NE.0) THEN

        IF (debug) THEN
           SELECT CASE (COORD)
           CASE ('ZMAT','MIXED','BAKER')
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
           CASE ('CART','HYBRID')
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
           CASE DEFAULT
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
              STOP
           END SELECT
        END IF ! matches if (debug) THEN

        ! we check that the maximum displacement is not greater than Smax
        If (debug) WRITE(*,*) "Fini before test:",Fini
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
        Fini=(Fini.AND.(NormStep.LE.SThresh))
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini

        GradTmp=GraDTmp*FactStep

        if (debug) WRITE(*,*) "DBG Main, L1475, FactStep=",FactStep
        if (debug.AND.(COORD.EQ.'ZMAT')) WRITE(*,'(A,12(1X,F10.6))') 'Scaled Step:',GradTmp(1:min(12,NFree))

        IF (debug) THEN
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
           ELSE
              DO Iat=1,Nat
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
                      XyzGeomF(IGeom,1:3,Iat)
              END DO
           END IF
        END IF ! matches if (debug) THEN

        SELECT CASE (COORD)
        CASE ('ZMAT')
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
           Idx=4
           DO Iat=4,Nat
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
              Idx=Idx+3
           END DO

        CASE ('BAKER')
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."

        CASE ('MIXED')
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
           DO Iat=1,NCart
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
           END DO
           Idx=3*NCart+1
           IF (NCart.EQ.1) THEN
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi

              Idx=7
           END IF
           IF (NCart.EQ.2) THEN
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
              Idx=9
           END IF

           
           DO Iat=max(NCart+1,4),Nat
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
              Idx=Idx+3
           END DO
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)

        CASE ('HYBRID','CART')
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
           STOP
        END SELECT

        IF (debug) THEN
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
           ELSE
              DO Iat=1,Nat
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
                      XyzGeomF(IGeom,1:3,Iat)
              END DO
           END IF
        END IF ! matches if (debug) THEN

        ! We project out the tangential part of the gradient to know if we are done
        GradTmp=Grad(IGeom,:)
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
        SELECT CASE (COORD)
        CASE ('ZMAT','MIXED','BAKER')
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
           Ca=S/(sqrt(Norm)*NormGrad)
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
        CASE ('CART','HYBRID')
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
           Ca=S/(sqrt(Norm)*NormGrad)
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
        CASE DEFAULT
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
           STOP
        END SELECT
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini

     END DO ! matches DO IGeom=2,NGeomF-1
     ! we save the old gradients
     GradOld=Grad
     EPathOld=Energies

     ! Shall we re-parameterize the path ?
     ! PFL 2007/Apr/10 ->
     ! We call PathCreate in all cases, it will take care of the 
     ! reparameterization, as well as calculating the tangents.
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
     XyzGeomI=XyzGeomF
     IntCoordI=IntCoordF

     Call PathCreate(IOpt)
     !     END IF
     ! <- PFL 2007/Apr/10

     if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)')  'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini
     IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN

        ! We have the new path, we calculate its energy and gradient

        Call EgradPath(IOpt,Flag_Opt_Geom)
        !IF(IOPT .EQ. 10) Then
        !	       Print *, 'Stopping at my checkpoint.'
        !           STOP !This is my temporary checkpoint.
        !ENDIF

        ! PFL 31 Mar 2008 ->
        ! Taken from v3.99 but we added a flag...
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
        ! else it is decreased by 0.8

        if (DynMaxStep) THEN
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
           if (OptReac) THEN
              If (Energies(1)<EPathOld(1)) Then
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
              ELSE
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
              END IF
           END IF

           if (OptProd) THEN
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
              ELSE
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
              END IF
           END IF

           DO IGeom=2,NGeomF-1
              If (Energies(IGeom)<EPathOld(IGeom)) Then
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
              ELSE
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
              END IF
           END DO
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
        END IF ! test on DynMaxStep
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
        ! <- PFL 31 MAr 2008
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
        !  It should get updated for other cases too.

        DO IGeom=1,NGeomF
           SELECT CASE (COORD)
           CASE('ZMAT')
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
              GeomTmp=IntCoordF(IGeom,:)
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
                   OrderInv(indzmat(2,2)),Geomtmp(1)
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
              Idx=4
              DO Iat=4,Nat
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
                 Idx=Idx+3
              END DO
           CASE('BAKER')
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
              !GeomTmp=IntCoordF(IGeom,:)
           CASE('CART','HYBRID','MIXED')
              WRITE(IOOUT,'(1X,I5)') Nat
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
              DO I=1,Nat
                 Iat=I
                 If (renum) Iat=OrderInv(I)
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
              END DO
           CASE DEFAULT
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
              STOP
           END SELECT
        END DO ! matches DO IGeom=1,NGeomF

        Call Write_path(IOpt)

        ! We update the Hessian matrices
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
        ! First using the displacement from the old path
        IGeom0=2
        IGeomF=NGeomF-1
        IF (OptReac) IGeom0=1
        IF (OptProd) IGeomF=NGeomF

        IF (FCalcHess) THEN
           DO IGeom=IGeom0,IGeomF
              SELECT CASE (COORD)
              CASE ('ZMAT','MIXED','BAKER')
                 GeomTmp2=IntCoordF(IGeom,:)
              CASE ('CART','HYBRID')
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
              CASE DEFAULT
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
                 STOP
              END SELECT
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
           END DO
        ELSE
           DO IGeom=IGeom0,IGeomF
              GeomTmp=GeomOld_all(IGeom,:)
              SELECT CASE (COORD)
              CASE ('ZMAT','MIXED','BAKER')
                 GeomTmp2=IntCoordF(IGeom,:)
              CASE ('CART','HYBRID')
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
              CASE DEFAULT
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
                 STOP
              END SELECT
              GeomTmp=GeomTmp2-GeomTmp
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
              HessTmp=Hess(IGeom,:,:)
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
              Hess(IGeom,:,:)=HessTmp
           END DO ! matches DO IGeom=IGeom0,IGeomF

           ! Second using the neighbour points
           IF (HupNeighbour) THEN
              ! Only one point for end points.
              IF (OptReac) THEN
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
                 ELSE
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
                 END IF
                 GradTmp=Grad(1,:)-Grad(2,:)
                 HessTmp=Hess(1,:,:)
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
                 Hess(1,:,:)=HessTmp
              END IF

              IF (OptProd) THEN
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
                 ELSE
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
                 END IF
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
                 HessTmp=Hess(NGeomF,:,:)
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
                 Hess(NGeomF,:,:)=HessTmp
              END IF

              ! Two points for the rest of the path
              DO IGeom=2,NGeomF-1
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
                 ELSE
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
                 END IF
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
                 HessTmp=Hess(IGeom,:,:)
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        

                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
                 ELSE
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
                 END IF
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)

                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
                 Hess(IGeom,:,:)=HessTmp
              END DO ! matches DO IGeom=2,NGeomF-1
           END IF !matches IF (HupNeighbour) THEN
        END IF ! matches IF (FCalcHess)
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))

  IF (PROG=="GAUSSIAN") THEN
     DEALLOCATE(Gauss_paste)
     DEALLOCATE(Gauss_root,Gauss_comment,Gauss_end,current)
  END IF
  DEALLOCATE(XyzGeomF, IntCoordF)
  DEALLOCATE(XyzGeomI, IntCoordI)
  DEALLOCATE(XyzTangent,IntTangent)
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
  DEALLOCATE(GeomOld_all)

  if (PROG=="GAUSSIAN") THEN
     ! We de-allocate the Gauss_XX lists
     ! transverse the list and deallocate each node
     previous => Gauss_root ! make current point to head of list
     DO
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
        current => previous%next ! make list point to next node of head
        DEALLOCATE(previous) ! deallocate current head node
        previous => current  ! make current point to new head
     END DO

     previous => Gauss_comment ! make current point to head of list
     DO
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
        current => previous%next ! make list point to next node of head
        DEALLOCATE(previous) ! deallocate current head node
        previous => current  ! make current point to new head
     END DO

     previous => Gauss_end ! make current point to head of list
     DO
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
        current => previous%next ! make list point to next node of head
        DEALLOCATE(previous) ! deallocate current head node
        previous => current  ! make current point to new head
     END DO

     DEALLOCATE(current)
  END IF

  DEALLOCATE(GeomTmp, Hess, GradTmp)
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
  IF (COORD.EQ.'BAKER') THEN
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
  END IF

  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt

END PROGRAM PathOpt
