Chimie Théorique » Opt'n Path

! This programs generates and optimizes a reaction path

! between at least two endpoints.

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs 

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

!

! 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 namelist.

! 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

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

! OpenPath v1.4  2011 May-June

! - Also added two flags: CalcEReac and CalcEProd. If .TRUE. Path will

! compute the energy of the reactant

! (and/or Product) at the begining of the calculations.

! This is useful if you start from non optimized points.

!

! - Also added more methods for Hessian updates. The following methods are

! now availables:

!   For inverse Hessian: BFGS, DFP, MS. 

!   For Hessian: BFGS, DFP, MS, Bofill

! 

! -Small change: HesUpd is replaced by Hupdate

!

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

! TO DO

! - Trying to improve the output of Path. In particular, we now create

! a .datl file containing the image index, curvilinear coordinate and energy

! in the spirit of the AnaPath tool.

! AnaPath is still usefull as it can extract geometrical parameters on top

! of the energy

! - When using Hessian update, it is more efficient to use Cholesky 

! decomposition, and to update this decomposition rather than the Hessian.

! See Nocedal&Wright p141.

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

!

! Project Name has changed to Opt'n Path

! OpenPaht 2011 May-June -> Opt'n Path v1.4 

!

! 2013 Jan 18 Opt'n Path v1.47 (going to 1.5 slowly)

! We add Siesta (and start thinking of add CP2K)

!  as a new "energy engine" code

! TO DO to go to 1.5: finish cleaning...

!

! 2013 Feb Opt'n Path v1.48

! We add some keyword for more output for Geometry Optimizations

! GeomFile: name of the file to print out the geometries and their energy

! as Xyz file. (only format for now)

!

! 2013 Feb opt'n Path v1.49

! We add the possibility to analyze the geometries in terms

! of Bonds, Angles, dihedrals

! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path)

! This is basicaly the merging of AnaPath into Opt'n Path... at last !

! This is done by adding a Flag into the Path namelist:

! AnaGeom: logical, if true Path looks for the AnaList namelist

! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path

! will read it and print the values of the variable in a Path.dat file

! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. 

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

! 2014 March Opt'n Path v1.50 Released under AGPL

! We used this release to clean the source files (and delete some unused files!)

! The Chamfre reactive Force field is not public and has been removed:

! if you want it, contact optnpath@gmail.com to see how to get it.

! it was anyway an old version of CH4+Ni(111) with old parameters.

!

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

PROGRAM PathOpt

  use VarTypes

  use Path_module

  use Io_module

  IMPLICIT NONE

  CHARACTER(132) :: FileIn, FileOut, Line

  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF

  INTEGER(KINT) :: IOpt

  INTEGER(KINT) :: Idx

  INTEGER(KINT) :: N3at, NFree, Nfr

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

  REAL(KREAL) :: E,FactStep !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

  ! Maximum step allowed for this geometry

  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF

! these are used to read temporary coordinates

  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

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

  LOGICAL :: Debug, Fini,Flag_Opt_Geom

!  INTEGER(KINT), EXTERNAL :: Iargc

  INTEGER(KINT), EXTERNAL ::  ConvertNumAt

  INTEGER(KINT) :: IoS

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

  !

  ! 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

    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)

      Use VarTypes

      implicit none

! Input

      CHARACTER(*), INTENT(IN) :: Line

      CHARACTER(*), INTENT(IN) :: String

      LOGICAL, OPTIONAL, INTENT(IN) :: Case

      LOGICAL, OPTIONAL, INTENT(IN) :: Back

      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean

! Output

! the position of String in Line (the first one) unless Back is present

      INTEGER(KINT) :: Pos

    END FUNCTION InString

    SUBROUTINE die(routine, msg, file, line, unit)

      Use VarTypes

      Use io_module

      implicit none

      character(len=*), intent(in)           :: routine, msg

      character(len=*), intent(in), optional :: file

      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die

    SUBROUTINE Warning(routine, msg, file, line, unit)

      Use VarTypes

      Use io_module

      implicit none

      character(len=*), intent(in)           :: routine, msg

      character(len=*), intent(in), optional :: file

      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE Warning

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

       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &

       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &

       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile

  Namelist/cartlist/list

  Namelist/frozenlist/list

  Namelist/analist/nb

  Flag_Opt_Geom = .FALSE. ! Initialized.

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

  ! 

  ! We print the Version info in file

  !

  WRITE(*,'(1X,A)') "Opt'n Path v1.50 (c) PFL/PD 2003-2014"

  WRITE(*,'(1X,A)') "Contact: optnpath@gmail.com"

  ! 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(*,*) "Opt'n Path mini-help"

        WRITE(*,*) "--------------"

        WRITE(*,*) ""

        WRITE(*,*) "Use: Path Input_file Output_file"

        WRITE(*,*) "Input_file starts with a Namelist called path"

        WRITE(*,*) ""

! The following has been generated (March 19, 2013) using the Mini_help.tex file

! 1) Edit Mini_help.tex

! 2) latex Mini_help.tex

! 3) catdvi -e 1 -U  Mini_help.dvi | sed -re "s/\[U\+2022\]/-/g" | sed -re "s/([^^[:space:]])\s+/\1 /g" > file.txt

! 4) edit file.txt to keep only the following lines, and to reformat a bit

! 5) awk '{print "WRITE(*,*) \"" $0 "\""}' file.txt > help

! 5) cut&paste help in Path.f90 ...

        WRITE(*,*) " &path"

        WRITE(*,*) " nat=3, ! Number of atoms"

        WRITE(*,*) " ngeomi=3, ! Number of initial geometries"

        WRITE(*,*) " ngeomf=12, !Number of geometries along the path"

        WRITE(*,*) " OptReac=.T., ! Do you want to optimize the reactants ?"

        WRITE(*,*) " OptProd=.T., ! Optimize the products"

        WRITE(*,*) " coord='zmat', ! We use Z-matrix coordinates"

        WRITE(*,*) " maxcyc=31, ! Max number of iterations"

        WRITE(*,*) " IReparam=2,! re-distribution of points along the path every 2 iterations"

        WRITE(*,*) " ISpline=50, ! Start using spline interpolation at iteration 50"

        WRITE(*,*) " Hinv=.T. , ! Use inverse of the Hessian internally (default: T)"

        WRITE(*,*) " MW=T, ! Works in Mass Weighted coordiante (default T)"

        WRITE(*,*) " PathName='Path_HCN_zmat_test', ! Name of the file used for path outputs"

        WRITE(*,*) " prog='gaussian',! we use G03 to get energy and gradients"

        WRITE(*,*) " SMax=0.1 ! Displacement cannot exceed 0.1 atomic units (or mass weighted at. unit)"

        WRITE(*,*) " /"

        WRITE(*,*) "  3"

        WRITE(*,*) " Energy : 0.04937364"

        WRITE(*,*) " H 0.0000 0.0000 0.0340"

        WRITE(*,*) " C 0.0000 0.0000 1.1030"

        WRITE(*,*) " N 0.0000 0.0000 2.2631"

        WRITE(*,*) "  3"

        WRITE(*,*) " Energy : 0.04937364"

        WRITE(*,*) " H 0.0000 1.1000 1.1030"

        WRITE(*,*) " C 0.0000 0.0000 1.1030"

        WRITE(*,*) " N 0.0000 0.0000 2.2631"

        WRITE(*,*) "3"

        WRITE(*,*) " CNH"

        WRITE(*,*) "H 0.000000 0.000000 3.3"

        WRITE(*,*) "C 0.000000 0.000000 1.1"

        WRITE(*,*) "N 0.000000 0.000000 2.26"

        WRITE(*,*) "%chk=Test"

        WRITE(*,*) "#P AM1 FORCE"

        WRITE(*,*) " HCN est bien"

        WRITE(*,*) ""

        WRITE(*,*) "0,1"

        WRITE(*,*) "H 0.000000 0.000000 0.000000"

        WRITE(*,*) "C 0.000000 0.000000 1.000"

        WRITE(*,*) "N 0.000000 0.000000 3.00"

        WRITE(*,*) ""

        WRITE(*,*) "*******"

        WRITE(*,*) "* Compulsory variables are:"

        WRITE(*,*) "*******"

        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(*,*) "*******"

        WRITE(*,*) "* Other options:"

        WRITE(*,*) "*******"

        WRITE(*,*) "Input: String that indicates the type of the input geometries. Accepted values are: Cart (or"

        WRITE(*,*) "    Xmol or Xyz), Vasp, Turbomole or Siesta."

        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."

        WRITE(*,*) "    Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext."

        WRITE(*,*) ""

        WRITE(*,*) ""

        WRITE(*,*) "        - In case of a Gaussian calculations, input must be set to Cart. One example of a gaus-"

        WRITE(*,*) "          sian input should be added at the end of the input file.See example file Test_HCN_zmat_g03.path."

        WRITE(*,*) "        - In the case of a VASP calculation, if input is set to Cart, then the preamble of a"

        WRITE(*,*) "          VASP calculation s"

        WRITE(*,*) "        - Mopac: Examples using sequential call to MOPAC2009 to calculate the energies and"

        WRITE(*,*) "          forces along the path. hould be added at the end of the input file. See example file"

        WRITE(*,*) "          Test_VASP_cart.path. In the case of a VASP calculation, one should also give value"

        WRITE(*,*) "          of the RunMode variable ."

        WRITE(*,*) "        - In the case of a SIESTA calculation, an example of a Siesta input file should be added"

        WRITE(*,*) "          at the end of the input file. See Test_Siesta.path."

        WRITE(*,*) ""

        WRITE(*,*) "RunMode: When running on a multi-processor machine, this indicates wether Opt'n Path"

        WRITE(*,*) "    should calculate the energy and gradient of the whole path in parallel or not. User has"

        WRITE(*,*) "    two options. One is to calculate the energy and gradient of each point sequentially. This"

        WRITE(*,*) "    is usefull when running on one (or two) processors. In this case, RunMode should be put"

        WRITE(*,*) "    to SERIAL.When running in parallel with 8 or more processors, one can use VASP to"

        WRITE(*,*) "    calculate simultaneously the energies and gradients for all points, as in a normal NEB cal-"

        WRITE(*,*) "    culation. In this case, RunMode must be set to PARA. For now, this is usefull only for VASP."

        WRITE(*,*) ""

        WRITE(*,*) "ProgExe: Name (with full path) of the executable to be used to get energies and gradients."

        WRITE(*,*) "    For example, if VASP is used in parallel, one might have something like:"

        WRITE(*,*) "    ProgExe='/usr/local/mpich/bin/mpirun -machinefile machine -np 8 ~/bin/VASP_46'."

        WRITE(*,*) "    Another option that I use, is to put ProgExe='./run_vasp' and I put every option needed"

        WRITE(*,*) "    to run VASP into the run_vasp file."

        WRITE(*,*) "EReac: (REAL) By default, Opt'n Path does not compute the energy of the reactants and"

        WRITE(*,*) "    products. This thus indicates the reactants energy to Opt'n Path to have better plots at"

        WRITE(*,*) "    the end."

        WRITE(*,*) "EProd: (REAL) By default, Opt'n Path does not compute the energy of the reactants and"

        WRITE(*,*) "    products. This thus indicates the products energy to Opt'n Path to have better plots."

        WRITE(*,*) ""

        WRITE(*,*) "PathName: Prefix used to save the path. Default is Path."

        WRITE(*,*) "Poscar: string that will be used as the prefix for the different POSCAR files in a VASP calcu-"

        WRITE(*,*) "    lations. Usefull only if PathOnly=.TRUE., not used for internal calculations."

        WRITE(*,*) "CalcName: Prefix for the files used for the energy and gradient calculations"

        WRITE(*,*) ""

        WRITE(*,*) "ISuffix: Suffix for the input file used for energy and gradient calculations. The full inputfile"

        WRITE(*,*) "    name will thus be CalcName.ISuffix."

        WRITE(*,*) "OSuffix: Suffix for the output file used for energy and gradient calculations. The full outputfile"

        WRITE(*,*) "    name will thus be CalcName.OSuffix."

        WRITE(*,*) ""

        WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates. Valid only for"

        WRITE(*,*) "     Coord=Zmat, Hybrid or Mixed."

        WRITE(*,*) "Fact: REAL used to define if two atoms are linked. If d(A,B) =< fact* (rcov(A) + rcov(B)),"

        WRITE(*,*) "     then A and B are considered Linked."

        WRITE(*,*) ""

        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 'cart' array(see below) are"

        WRITE(*,*) "          describe in CARTESIAN, whereas the 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(*,*) ""

        WRITE(*,*) "        - RFO: Rational function optimization"

        WRITE(*,*) "        - GDIIS: Geometry optimization using direct inversion in the iterative supspace"

        WRITE(*,*) "HUpdate: method to update the hessian. By default, it is Murtagh-Sargent Except for geom-"

        WRITE(*,*) "     etry optimization where it is BFGS."

        WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"

        WRITE(*,*) ""

        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"

        WRITE(*,*) "     taken"

        WRITE(*,*) ""

        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, whereas FTan=0. gives a"

        WRITE(*,*) "     displacement orthogonal to the path."

        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of"

        WRITE(*,*) "     reparameterization."

        WRITE(*,*) "IReparamT: When the path is not reparameterised at each iteration, this gives the frequency"

        WRITE(*,*) "     of reparameterization of the tangents."

        WRITE(*,*) ""

        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path. This option indicates"

        WRITE(*,*) "     the first step where spline interpolation is used."

        WRITE(*,*) "BoxTol: Real between 0. and 1. When doing periodic calculations, it might happen that an"

        WRITE(*,*) "     atom moves out of the unit cell. Path detects this by comparing the displacement to"

        WRITE(*,*) "     boxtol: if an atom moves by more than Boxtol, then it is moved back into the unit cell."

        WRITE(*,*) "     Default value: 0.5."

        WRITE(*,*) "FrozTol: (Real) This indicates the threshold to determine wether an atom moves between two"

        WRITE(*,*) "     images. Default is 1e-4."

        WRITE(*,*) ""

        WRITE(*,*) "OptGeom: This INTEGER indicates the index of the geometry you want to optimize. If"

        WRITE(*,*) "    OptGeom is set, then Opt'n Path performs a geometry optimization instead of a path"

        WRITE(*,*) "    interpolation."

        WRITE(*,*) "GeomFile: Name of the file to print the geometries and their energy during a geometry opti-"

        WRITE(*,*) "    mization. If this variable is not given then nothing is printed."

        WRITE(*,*) "AnaFile: Name of the file to print the values of the geometrical parameters that are monitored"

        WRITE(*,*) "    if AnaGeom=.TRUE.. Default is PathName.dat"

        WRITE(*,*) ""

        WRITE(*,*) "GplotFile: Name of the gnuplot file to plot the evolution of the geometrical parameters that"

        WRITE(*,*) "    are monitored if AnaGeom=.TRUE.. Default is PathName.gplot"

        WRITE(*,*) ""

        WRITE(*,*) "*******"

        WRITE(*,*) "* Arrays:"

        WRITE(*,*) "*******"

        WRITE(*,*) ""

        WRITE(*,*) "Rcov: Array containing the covalent radii of the first 86 elements. You can modify it using,"

        WRITE(*,*) "    rcov(6)=0.8."

        WRITE(*,*) "Mass: Array containing the atomic mass of the first 86 elements."

        WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations. If not given, Opt'n Path will"

        WRITE(*,*) "    read the POTCAR file."

        WRITE(*,*) ""

        WRITE(*,*) "*******"

        WRITE(*,*) "* Flags:"

        WRITE(*,*) "*******"

        WRITE(*,*) ""

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

        WRITE(*,*) "    unless for Coord equals CART."

        WRITE(*,*) ""

        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(*,*) "CalcEProd: if TRUE the product energy will be computed. Default is FALSE. Not Compatible"

        WRITE(*,*) "    with RunMode=Para."

        WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compat-"

        WRITE(*,*) "    ible with RunMode=Para."

        WRITE(*,*) ""

        WRITE(*,*) "PathOnly: TRUE if one wants to generate the initial path, and stops."

        WRITE(*,*) "Align: If .FALSE., successive geometries along the path are not aligned on each other before"

        WRITE(*,*) "    path interpolation. Default is .TRUE. if there are 4 atoms or more."

        WRITE(*,*) "Hinv: if True, then Hessian inversed is used."

        WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."

        WRITE(*,*) ""

        WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of"

        WRITE(*,*) "    the path."

        WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms. If True, a &frozenlist"

        WRITE(*,*) "    namelist containing the list of frozen atoms must be given. If VASP is used, and frozen is"

        WRITE(*,*) "    not given 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. *** Only used in"

        WRITE(*,*) "    'mixed' calculations. *** If True, a &cartlist namelist containing the list of cart atoms"

        WRITE(*,*) "    must be given. By default, only frozen atoms are described in cartesian coordinates."

        WRITE(*,*) ""

        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(*,*) ""

        WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files. See also"

        WRITE(*,*) "      the POSCAR option. This can be used only if prog or input=VASP."

        WRITE(*,*) "AnaGeom: If TRUE, Opt'n Path will create a file .dat for geometries analysis. If True, "

        WRITE(*,*) "Opt'n Path will look for the AnaList namelist after the Path Namelist."

        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. It is ensured that the"

        WRITE(*,*) "      dynamical Smax is within [0.5*SMax_0 ,2*Smax_0 ]."

        WRITE(*,*) ""

        WRITE(*,*) "*******"

        WRITE(*,*) "* Additional namelists:"

        WRITE(*,*) "*******"

        WRITE(*,*) ""

        WRITE(*,*) "&cartlist list= ... &end This gives the list of atoms that are described using cartesian coor-"

        WRITE(*,*) "      dinates. Read only if FCart=.TRUE.. To indicate an atom range, from 1 to 5 for example,"

        WRITE(*,*) "      one can put 1 -5 in the list. For example:"

        WRITE(*,*) "      &cartlist list = 1 2 6 12 -20 &end"

        WRITE(*,*) "      will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian."

        WRITE(*,*) "&Frozenlist list= ... &end This gives the list of atoms that are frozen during optimization."

        WRITE(*,*) "      Read only if FFrozen=.TRUE.. To indicate an atom range, from 1 to 5 for example, one"

        WRITE(*,*) "      can put 1 -5 in the list."

        WRITE(*,*) ""

        WRITE(*,*) "&Analist nb= ... &end list of variables for geometry analysis. If present and if Ana-"

        WRITE(*,*) "      Geom=T then Opt'n Path will read it and print the values of the variable in a .dat file If"

        WRITE(*,*) "      AnaGeom is T but Analist is not given, then Path.dat just contains the energy. Analist"

        WRITE(*,*) "      contains the number of variables to monitor: Nb, and is followed by the description of the"

        WRITE(*,*) "      variables among: b(ond) At1 At2 a(ngle) At1 At2 At3 d(ihedral) At1 At2 At3 At4 c NbAt"

        WRITE(*,*) "      At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added at the"

        WRITE(*,*) "      end of the real atoms of the system."

!!!!!!!!!! end of included file

        WRITE(*,*) ""

        STOP

     END IF

     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')

     IOOUT=6

  CASE (0)

     IOIN=5

     IOOUT=6

  END SELECT

! We set the file delimiter (system dependent)

  Call Set_FileDelim

  ! 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="EMPTY"

  HUpdate="EMPTY"

  FFrozen=.FALSE.

  FCart=.FALSE.

  List=0

  renum=.TRUE.

  OptReac=.FALSE.

  OptProd=.FALSE.

  CalcEReac=.FALSE.

  CalcEProd=.FALSE.

  EReac=0.d0

  EProd=0.d0 

  OptGeom=-1

  GeomFile="EMPTY"

  AnaGeom=.FALSE.

  AnaFile="EMPTY"

  GplotFile="EMPTY"

  Nb=0

  NbCom=0

  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").or.valid("Path")

  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 ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN

     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"

     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"

     CalcEReac=.FALSE.

     CalcEProd=.FALSE.

  END IF

  If (IReparamT<=0) IReparamT=IReparam

! We take care of HesUpd flag

! this is needed because some users might still use it

  if (HesUpd.NE."EMPTY") THEN

! We are pityless:

     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"

     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"

     STOP

  END IF

  IF (HUpdate=='EMPTY') THEN

     IF (OptGeom>=1) THEN

        HupDate="BFGS"

     ELSE

        HUpdate="MS"

     END IF

  END IF

  If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 

     COORD='CART'

  END IF

  IF (COORD.EQ.'CART') THEN

     Renum=.FALSE.

     IF (OptGeom>0) THEN

        IGeomRef=OptGeom

     ELSE

        IGeomRef=1

     END IF

  END IF

  if (Prog.EQ.'TEST2'.OR.Prog.EQ.'CHAMFRE') Then

     WRITE(IOOUT,*) 'The CHAMFRE Reactive Force Field is not public.'

     WRITE(IOOUT,*) 'Contact optnpath@gmail.com to get it.'

     STOP

  END IF

  if (Prog.EQ.'2D') Prog="TEST2D"

  if (Prog.EQ.'TEST3') Prog="LEPS"

  Select Case (Prog)

    CASE ("EMPTY","G09","GAUSSIAN")

     Prog='GAUSSIAN'

     If (ProgExe=="EMPTY") ProgExe="g09"

     if (CalcName=="EMPTY") CalcName="Path"

     if (ISuffix=="EMPTY")  ISuffix="com"

     if (OSuffix=="EMPTY")  OSuffix="log"

     UnitProg="au"

     ConvE=au2kcal

    CASE ("G03")

     Prog='GAUSSIAN'

     If (ProgExe=="EMPTY") ProgExe="g03"

     if (CalcName=="EMPTY") CalcName="Path"

     if (ISuffix=="EMPTY")  ISuffix="com"

     if (OSuffix=="EMPTY")  OSuffix="log"

     UnitProg="au"

     ConvE=au2kcal

    CASE ("EXT")

     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"

     if (CalcName=="EMPTY") CalcName="Ext"

     if (ISuffix=="EMPTY")  ISuffix="in"

     if (OSuffix=="EMPTY")  OSuffix="out"

     UnitProg="au"

     ConvE=au2kcal

    CASE ('MOPAC')

     If (ProgExe=="EMPTY") ProgExe="mopac"

     if (CalcName=="EMPTY") CalcName="Path"

     if (ISuffix=="EMPTY")  ISuffix="mop"

     if (OSuffix=="EMPTY")  OSuffix="out"

     UnitProg="au"

     ConvE=au2kcal

    CASE ("SIESTA")

     If (ProgExe=="EMPTY") ProgExe="siesta"

     if (CalcName=="EMPTY") CalcName="Path"

     if (ISuffix=="EMPTY")  ISuffix="fdf"

     if (OSuffix=="EMPTY")  OSuffix="out"

     UnitProg="eV"

     ConvE=eV2kcal

    CASE ("VASP")

     If (ProgExe=="EMPTY") ProgExe="VASP"

     UnitProg="eV"

     ConvE=eV2kcal

     FPBC=.TRUE.

     IPer=3

     Kabeg=-1

     kaEnd=1

     kbBeg=-1

     kbEnd=1

     kcBeg=-1

     kcEnd=1

    CASE ("TM","TURBOMOLE")

     Prog="TURBOMOLE"

     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"

     UnitProg="au"

     ConvE=au2kcal

     FPBC=.FALSE.

    CASE ("TEST","CHAMFRE","TEST2D","LEPS")

     ProgExe=""

     UnitProg="au"

     ConvE=au2kcal

    CASE DEFAULT

     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"

     STOP

  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

     FPrintGeom=.FALSE.

     If (GeomFile/='EMPTY') THEN

        FPrintGeom=.TRUE.

     END IF

     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

  ALLOCATE(Cart(Nat))

  Cart=0

  ALLOCATE(Frozen(Nat))

  Frozen=0

  IF (FCart) THEN

! We rewind the file to be sure to browse all of it to read the Cart namelist

     REWIND(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

  END IF

  if (FFrozen) THEN

     REWIND(IOIN)

     List=0

     READ(IOIN,Frozenlist)

     Frozen=List(1:Nat)

     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen

  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 (AnaGeom) THEN

     REWIND(IOIN)

     READ(IOIN,AnaList,IOSTAT=IoS)

     IF (IOS==0) THEN ! we have a AnaList namelist 

        Call ReadAnaList

     ELSE

! No AnaList namelist, we will print only the energy

        FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)'

     END IF

     IF (AnaFile=="EMPTY") THEN

        AnaFile=TRIM(PathName) // '.dat'

     END IF

     OPEN(IODat,File=AnaFile)

     IF (GplotFile=="EMPTY") THEN

        GplotFile=TRIM(PathName) // '.gplot'

     END IF

     OPEN(IODat,File=AnaFile)

     If (OptGeom<=0) THEN

        OPEN(IOGplot,File=GPlotFile)

        WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist"

        WRITE(IoGplot,'(A)') " set pointsize 2"

        WRITE(IoGplot,'(A)') " xcol=1"

        WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3

     END IF

  END IF

  if (NMaxPtPath.EQ.-1) then

     NMaxPtPath=min(50*NGeomF,2000)

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

  AtName=""

  MassAt=1.

! As we have used rewind many times, the position in the input file

! is blurry (at least !) we thus have to go back to the Geometries description

! (TO DO: use a better Input parser !)

  REWIND(IOIN)

! We search the '&path' namelist

  READ(IOIN,'(A)',Iostat=Ios) Line

  Fini=.FALSE.

  DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0)

     if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line)

     READ(IOIN,'(A)',Iostat=Ios) Line

  END DO

  if (debug) WRITE(*,*) "Path found:",TRIM(LINE)

  ! We are on the &path line, we move to the end

  Call NoComment(Line)

  DO WHILE ((Ios==0).AND.InString(Line,"/")==0)

     if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)

     READ(IOIN,'(A)',Iostat=Ios) Line

     Call NoComment(Line)

     Call noString(Line)

  END DO

! We are at then end of &Path / namelist

! We scan for other Namelists

  READ(IOIN,'(A)',Iostat=Ios) Line

  if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)

  DO WHILE  ((IoS==0).AND.(Index(Line,'&')/=0))

     Call NoComment(Line)

     if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)

     Idx=InString(Line,'Analist')

     DO WHILE ((Ios==0).AND.InString(Line,"/")==0)

        if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)

        READ(IOIN,'(A)',Iostat=Ios) Line

        Call NoComment(Line)

     END DO

     If (Idx/=0) THEN

! we have just read &Analist namelist, we have to skip the variables description

        DO I=1,Nb

           READ(IOIN,'(A)',Iostat=Ios) Line

           if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE)

        END DO

     END IF

     READ(IOIN,'(A)',Iostat=Ios) Line

     if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE)

  END DO

  BACKSPACE(IOIN)

  ! 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

 If (AnaGeom) THEN

    If (NbVar>0) THEN

! We print the description of the varialbes in AnaFile

       Call PrintAnaList(IoDat)

       If (OptGeom<=0) THEN

          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"

       ELSE

          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"

       END IF

    ELSE

       If (OptGeom<=0) THEN

          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"

       ELSE

          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"       

       END IF

    END IF

 END IF

  ! 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

        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

        Iat3=1

        DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat))

           Iat3=Iat3+1

        END DO

        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 used 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.

  Call ReadInput

  if (FPBC)  THEN

     Allocate(XGeomRefPBC(Nat),YGeomRefPBC(Nat),ZGeomRefPBC(Nat))

     XGeomRefPBC(1:Nat)=XyzGeomI(1,1,1:Nat)

     YGeomRefPBC(1:Nat)=XyzGeomI(1,2,1:Nat)

     ZGeomRefPBC(1:Nat)=XyzGeomI(1,3,1:Nat)

     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

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

  !

  !

  ! 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(XyzGeomI(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

     if (debug) 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

     Call Header("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)

  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt

  Call Write_path(IOpt)

! Shall we analyze the geometries ?

     IF (AnaGeom) Call AnaPath(Iopt,IoDat)

  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