Chimie Théorique » OpenPath

! 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

! these are used to read temporary coordinates

  REAL(KREAL) :: xtmp,ytmp,ztmp

  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)') "OpenPath v1.0 (c) PFL/PD 2007-2011"

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

! PFL 2011 Mar 14 ->

! Added some printing for analyses with Anapath

  if (debug) THEN

     WRITE(IOOUT,path)

  ELSE

! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC

! and PATHNAME

     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat

     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI

     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF

     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC

     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)

  END IF

  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