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