/ - Diff - Opt'n Path - Forge du Centre Blaise Pascal

Révision 9

     \documentclass[a4paper,11pt]{article}
     \usepackage{a4wide}
     %\usepackage{times}
     \usepackage{graphicx}
     \usepackage[body={16cm,24.5cm}]{geometry}
-...
     \lhead[]{}
     \rhead[]{}
     \cfoot{}
     \title{OpenPath mini-help}
     \title{Opt'n Path mini-help}
     \author{P. Fleurat-Lessard, P. Dayal \\
     Laboratoire de Chimie de l'ENS Lyon, 46 all?e d'Italie, F-69364 Lyon
     Laboratoire de Chimie de l'ENS Lyon, 46 all\'ee d'Italie, F-69364 Lyon
     Cedex 7}
     \date{Feb. 2011}
     \date{March 2013}
     \def\Path{\texttt{OpenPath}}
     \def\Path{\texttt{Opt'n Path}}
     \begin{document}
     \maketitle
     \section{introduction}
     \section{Introduction}
     \Path{} is a program that can optimize a reaction path between two
       structures. The algorithm to optimize the path is close to the string
       method. The originality of this program lies in the coordinate set it
-...
     This program is an independant program that calls standard electronic
       structure codes to get the energies and forces it needs to optimize
       the reaction path. For now, it is coupled to Gaussian, MOPAC2009,
       Vasp and Turbomole.
       the reaction path. For now, it is coupled to Gaussian, MOPAC,
       Vasp, Turbomole and Siesta.
     \section{Installation}
     We suppose here that you will install \Path{} into a directory called
     \texttt{OpenPath}. First, create the directory:
     \texttt{optnpath}. First, create the directory:
     \begin{verbatim}
      mkdir OpenPath
      cd OpenPath
      mv ../OpenPath.tgz .
      mkdir optnpath
      cd optnpath
      mv ../optnpath_x.yy.tgz .
     \end{verbatim}
     Uncompress the archive:
     \begin{verbatim}
      gunzip OpenPath.tgz
      tar -xvf OpenPath.tar
      gunzip optnpath_x.yy.tgz
      tar -xvf Optnpath.tar
     \end{verbatim}
     You should now have this \texttt{Mini\_help.pdf} file and 4 directories (doc, src, utils, examples).
-...
     \end{verbatim}
        Compulsory variables are:
     \subsubsection{Compulsory variables are:}
     \begin{description}
     \item          NGeomi: Number of geometries defining the Initial path
     \item          NGeomf: Number of geometries defining the Final path
     \item         Nat   : Number of atoms
     \item[NGeomi:] Number of geometries defining the Initial path
     \item[NGeomf:] Number of geometries defining the Final path
     \item[Nat:] Number of atoms
     \end{description}
               Other options
     \subsubsection{Other options:}
     \begin{description}
     \item  Input: string that indicates the type of the input geometries.
               Accepted values are: Cart (or Xmol or Xyz) or Vasp
     \item    Prog: string that indicates the program that will be used for energy and gradient calculations.
                    Accepted values are: Gaussian, Vasp or Ext. \\
                     In case of a Gaussian calculations, input must be set to Cart.
     \item[Input:] String that indicates the type of the input geometries.
               Accepted values are: Cart (or Xmol or Xyz), Vasp, Turbomole or Siesta.
     \item[Prog:] string that indicates the program that will be used for energy and gradient calculations.
                    Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext. \\
     \begin{itemize}
     \item   In case of a Gaussian calculations, input must be set to Cart.
                     One example of a gaussian input should be added at the end of the
                     input file.See example file \texttt{Test\_HCN\_zmat\_g03.path}. \\
                     In the case of a VASP calculation, if input is set to Cart, then
                   the preamble of a VASP calculation should be added at the end of
     \item       In the case of a VASP calculation, if input is set to Cart, then
                   the preamble of a VASP calculation s\item \texttt{Mopac}: Examples  using sequential call
       to MOPAC2009 to calculate the energies and forces along the path.
     hould be added at the end of
                   the input file. See example file \texttt{Test\_VASP\_cart.path}.
                   In the case of a VASP calculation, one should also give value of the
                   \texttt{RunMode} variable .
     \item  In the case of a SIESTA calculation, an example of a Siesta input file
          should be added at the end of the input file. See \texttt{Test\_Siesta.path}.
     \end{itemize}
     \item RunMode: This indicates wether \Path{} should use VASP
                   routine to calculate the energy and gradient of the whole path
     \item[RunMode:] When running on a multi-processor machine, this indicates
       wether \Path{} should calculate the energy and gradient of the whole path in parallel
                   or not. User has two options. One is to calculate the energy and
                   gradient of each point sequentially. This is usefull when
                   running on one (or two) processors. In this case,
-...
                   in parallel with 8 or more processors, one can use VASP to
                   calculate simultaneously the energies and gradients for all
                   points, as in a normal NEB calculation. In this case,
                   \texttt{RunMode} must be set to \texttt{PARA}.
     \item ProgExe: Name (with full path) of the executable to be used to get
                   \texttt{RunMode} must be set to \texttt{PARA}.
      \emph{For now, this is usefull only for VASP.} \\
     \item[ProgExe:] Name (with full path) of the executable to be used to get
       energies and gradients. For example, if VASP is used in parallel, one might
       have something like: \\
     \verb!ProgExe='/usr/local/mpich/bin/mpirun -machinefile machine -np 8 ~/bin/VASP_46'!.
     Another option that I use, is to put \verb!ProgExe='./run_vasp'! and I put every option needed to run VASP
     into the \texttt{run\_vasp} file.
     \item   PathName: Prefix used to save the path. Default is Path
     \item   Poscar: string that will be used as the prefix for the different
     \item[EReac:] (REAL) By default, \Path{} does not compute the energy of the reactants
      and products. This thus indicates the reactants energy to \Path{} to have better plots
      at the end.
     \item[EProd:] (REAL) By default, \Path{} does not compute the energy of the reactants
      and products. This thus indicates the products energy to \Path{} to have better plots.
     \item[PathName:] Prefix used to save the path. Default is Path
     \item[Poscar:] string that will be used as the prefix for the different
                POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.,
                not used for internal calculations.
     \item   IGeomRef: Index of the geometry used to construct the internal coordinates.
              Valid only for Coord=Zmat, Hybrid or Mixed
     \item   Fact: REAL used to define if two atoms are linked.
     \item[CalcName:] Prefix for the files used for the energy and gradient calculations
     \item[ISuffix:] Suffix for the input file used for energy and gradient calculations.
     The full inputfile name will thus be \textit{CalcName}.\textit{ISuffix}.
     \item[OSuffix:] Suffix for the output file used for energy and gradient calculations.
      The full outputfile name will thus be \textit{CalcName}.\textit{OSuffix}.
     \item[IGeomRef:] Index of the geometry used to construct the internal coordinates.
              Valid only for Coord=Zmat, Hybrid or Mixed.
     \item[Fact:] REAL used to define if two atoms are linked.
              If $d(A,B) \leq fact*(rcov(A)+rcov(B))$, then A and B are considered Linked.
     \item   debugFile: Name of the file that indicates which subroutine should print debug info.
     \item   Coord: System of coordinates to use. Possible choices are:
               - CART (or Xyz): works in cartesian
               - Zmat: works in internal coordinates (Zmat)
               - Mixed: frozen atoms, as well as atoms defined by the
     \item[debugFile:] Name of the file that indicates which subroutine should print debug info.
     \item[Coord:] System of coordinates to use. Possible choices are:
     \begin{itemize}
     \item CART (or Xyz): works in cartesian
     \item Zmat: works in internal coordinates (Zmat)
     \item Mixed: frozen atoms, as well as atoms defined by the
               'cart' array(see below) are describe in CARTESIAN, whereas the
               others are described in Zmat
               - Baker: use of Baker coordinates, also called delocalized internal coordinates
               - Hybrid: geometries are described in zmat, but the gradient are used in cartesian
     \item   Step\_method: method to compute the step for optimizing a geometry; choices are:
               - RFO: Rational function optimization
               - GDIIS: Geometry optimization using direct inversion in the iterative supspace
     \item    HUpdate: method to update the hessian. By default, it is Murtagh-Sargent
     \item  Baker: use of Baker coordinates, also called delocalized internal coordinates
     \item Hybrid: geometries are described in zmat, but the gradient are used in cartesian
     \end{itemize}
     \item[Step\_method:] method to compute the step for optimizing a geometry; choices are:
     \begin{itemize}
     \item RFO: Rational function optimization
     \item GDIIS: Geometry optimization using direct inversion in the iterative supspace
     \end{itemize}
     \item[HUpdate:] method to update the hessian. By default, it is Murtagh-Sargent
              Except for geometry optimization where it is BFGS.
     \item   MaxCyc: maximum number of iterations for the path optimization
     \item   Smax: Maximum length of a step during path optimization
     \item   SThresh: Step Threshold to consider that the path is stationary
     \item   GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken
     \item   FTan: We moving the path, this gives the proportion of the displacement tangent to the path
     \item[MaxCyc:] maximum number of iterations for the path optimization
     \item[Smax:] Maximum length of a step during path optimization
     \item[SThresh:] Step Threshold to consider that the path is stationary
     \item[GThresh:] Gradient Threshold to consider that the path is stationary, only orthogonal part is taken
     \item[FTan:] We moving the path, this gives the proportion of the displacement tangent to the path
              that is kept. FTan=1. corresponds to the full displacement,
              whereas FTan=0. gives a displacement orthogonal to the path.
     \item   IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization.
     \item   ISpline: By default, a linear interpolation is used to generate the path.
     \item[IReparam:] The path is not reparameterised at each iteration. This gives the
     frequency of reparameterization.
     \item[IReparamT:] When the path is not reparameterised at each iteration, this gives the
     frequency of reparameterization of the \emph{tangents}.
     \item[ISpline:] By default, a linear interpolation is used to generate the path.
                 This option indicates the first step where spline interpolation is used.
     \item[BoxTol:]  Real between 0. and 1. When doing periodic calculations,
             it might happen that an atom moves out of the unit cell.
             Path detects this by comparing the displacement to boxtol:
             if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5.
     \item[FrozTol:]  (Real) This indicates the threshold to determine wether an atom moves between two images. Default is 1e-4.
     \item[OptGeom:] This INTEGER indicates the index of the geometry you want to optimize.
     If \texttt{OptGeom} is set, then \Path{} performs a geometry optimization instead of a path interpolation.
     \item[GeomFile:]  Name of the file to print the geometries and their energy
           during a geometry optimization. If this variable is not given
           then nothing is printed.
     \item[AnaFile:]  Name of the file to print the values of the geometrical parameters
     that are monitored if \texttt{AnaGeom=.TRUE.}. Default is \textit{PathName}.dat
     \item[GplotFile:]  Name of the \texttt{gnuplot} file to plot the evolution of the geometrical parameters that are monitored if \texttt{AnaGeom=.TRUE.}. Default is \textit{PathName}.gplot
     %% Not described here: NMaxPtPath, NGintMax (too technical ?)
     \end{description}
               Arrays:
     \subsubsection{Arrays:}
     \begin{description}
     \item   Rcov: Array containing the covalent radii of the first 80 elements.
     \item[Rcov:] Array containing the covalent radii of the first 86 elements.
              You can modify it using, \verb!rcov(6)=0.8!.
     \item   Mass: Array containing the atomic mass of the first 80 elements.
     \item   AtTypes: Name of the different atoms used in a VASP calculations.
     \item[Mass:] Array containing the atomic mass of the first 86 elements.
     \item[AtTypes:] Name of the different atoms used in a VASP calculations.
               If not given, Path will read the POTCAR file.
     \end{description}
               Flags:
     \subsubsection{Flags:}
     \begin{description}
     \item        MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE.
     \item       Renum: Flag. True if one wants to reoder the atoms in the
       initial order. default is .TRUE. most of the time.
     \item       OptProd: True if one wants to optimize the geometry of the products.
     \item       OptReac: True if one wants to optimize the geometry of the reactants.
     \item       CalcEProd: if TRUE the product energy will be
       computed. Default is FALSE. Not Compatible with RunMode=Para".
     \item       CalcEReac: if TRUE the reactants energy will be
       computed. Default is FALSE. Not Compatible with RunMode=Para".
     \item       PathOnly:TRUE if one wants to generate the initial path, and stops.
     \item       Hinv: if True, then Hessian inversed is used.
     \item       IniHup: if True, then Hessian inverse is extrapolated
     \item[MW:]  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE.
     \item[Renum:] Flag. True if one wants to reoder the atoms in the
       initial order. default is .TRUE. unless for \texttt{Coord} equals \texttt{CART}.
     \item[OptProd:] True if one wants to optimize the geometry of the products.
     \item[OptReac:] True if one wants to optimize the geometry of the reactants.
     \item[CalcEProd:] if TRUE the product energy will be
       computed. Default is FALSE. Not Compatible with \texttt{RunMode=Para}.
     \item[CalcEReac:] if TRUE the reactants energy will be
       computed. Default is FALSE. Not Compatible with \texttt{RunMode=Para}.
     \item[PathOnly:] TRUE if one wants to generate the initial path, and stops.
     \item[Align:] If .FALSE., successive geometries along the path are not aligned on each other before path interpolation. Default is .TRUE. if there are 4 atoms or more.
     \item[Hinv:] if True, then Hessian inversed is used.
     \item[IniHup:] if True, then Hessian inverse is extrapolated
       using the initial path calculations.
     \item       HupNeighbour: if True, then Hessian inverse is
     \item[HupNeighbour:]  if True, then Hessian inverse is
       extrapolated using the neighbouring points of the path.
     \item       FFrozen: True if one wants to freeze the positions of some atoms.
     \item[FFrozen:] True if one wants to freeze the positions of some atoms.
                     If True, a \verb!&frozenlist! namelist containing the
                     list of frozen atoms must be given.
                      If VASP is used, and frozen is not given
             here, the program will use the F flags of the input geometry
     \item       FCart:  True if one wants to describe some atoms using cartesian coordinates.
     \item[FCart:]  True if one wants to describe some atoms using cartesian coordinates.
                    *** Only used in 'mixed' calculations. ***
                  If True, a \verb!&cartlist! namelist containing the list of cart atoms must be given.
                  By default, only frozen atoms are described in cartesian coordinates.
     \item       Autocart: True if you want to let the program choosing the cartesian atoms.
     \item       VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now
     \item[Autocart:] True if you want to let the program choosing the cartesian atoms.
     \item[VMD:] TRUE if you want to use VMD to look at the Path. Used only for VASP for now.
     \item[WriteVASP:] TRUE if you want to print the images coordinates in POSCAR files.
     See also the POSCAR option. This can be used only if prog or input=VASP.
     \item[AnaGeom:] If TRUE, Opt'n Path will create a file .dat for geometries analysis.
             If True, Opt'n Path will look for the AnaList namelist after the Path Namelist.
     \item[DynMaxStep:] if TRUE, the maximum allowed step is updated at each step, for each geometry.
             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]$.
     \end{description}
     \subsubsection{Additional namelists:}
     \begin{description}
     \item[\&cartlist list= \ldots{} \&end] This gives the list of atoms that are described using cartesian coordinates. Read only if \texttt{FCart=.TRUE.}. To indicate an atom range, from 1 to 5 for example, one can put 1 -5 in the list. For example: \\
     \texttt{\&cartlist list = 1 2 6 12 -20 \&end} \\
     will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian.
     \item[\&Frozenlist list= \ldots{}  \&end] This gives the list of atoms that are frozen during optimization. Read only if \texttt{FFrozen=.TRUE.}. To indicate an atom range, from 1 to 5 for example, one can put 1 -5 in the list.
     \item[\&Analist nb= \ldots{} \&end] list of variables for geometry analysis. If present and if AnaGeom=T then
          Opt'n Path will read it and print the values of the variable in a .dat file
         If AnaGeom is T but Analist is not given, then Path.dat just contains the energy.
     Analist contains the number of variables to monitor: Nb, and is followed by the description
     of the variables among:
     b(ond) At1 At2
     a(ngle) At1 At2 At3
     d(ihedral) At1 At2 At3 At4
     c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added
     at the end of the real atoms of the system.
     \end{description}
     \subsection{Examples}
     More to come...      have a look at the files provided in the
     \texttt{examples} directory. In particular, you will find there three
     directories containing easy-to-use examples. That is, each directory contains
-...
     \end{itemize}
     \end{itemize}
     \item \texttt{Mopac}: Examples  using sequential call
       to MOPAC2009 to calculate the energies and forces along the path.
       to MOPAC to calculate the energies and forces along the path.
     \item \texttt{Siesta}: Examples  using sequential call
       to Siesta to calculate the energies and forces along the path.
     \item \texttt{Test}: Examples  using the analytical HCN potential
       energy surface to calculate the energies and forces along the path.
     As this is fast, we also provide other Analysis tools.

             END IF
             GeomCart=0.d0
             IF(IOpt .EQ. 1) Then
                Print *, 'GeomTmp='
                WRITE(*,'(12(1X,F6.3))') GeomTmp
                Print *, 'GeomCart='
                WRITE(*,'(12(1X,F6.3))') GeomCart
           END IF
     !        IF(IOpt .EQ. 1) Then
     !           Print *, 'GeomTmp='
     !           WRITE(*,'(12(1X,F6.3))') GeomTmp
     !           Print *, 'GeomCart='
     !           WRITE(*,'(12(1X,F6.3))') GeomCart
     !      END IF
            ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
            ! GeomCart has INTENT(OUT)

                Idx=Index(Line," ")
                Line2=Trim(AdjustL(Line(Idx+1:)))
                Read(Line2,*) Siesta_NbSpecies
                ALLOCATE(ListSpecies(Siesta_NbSpecies))
                ALLOCATE(Siesta_SpeciesMass(Siesta_NbSpecies))
                ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies))
             END IF
-...
                   Read(Line,*) ITmp
                   Idx=Index(Line,' ')
                   Line=AdjustL(Line(Idx+1:))
                   Read(Line,*) Ztmp
                   ListSpecies(ITmp)=ZTmp
                   Read(Line,*) Idx
                   Siesta_SpeciesMass(ITmp)=Idx
                   Idx=Index(Line,' ')
                   Line=AdjustL(Line(Idx+1:))
                   Idx=Index(Line,' ')
-...
                              TRIM(Siesta_SpeciesName(IdxSpecies(I)))
                         Call Die('Readinput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
                      END IF
     ! We look for the Mass number also
                      IF (Atome(I)/=Siesta_SpeciesMass(IdxSpecies(I))) THEN
                         Write(Line,'(A,I5,1X,I5,1X,I5)') 'AtMass /= SpeciesMass for atom ',I,Atome(I), &
                              Siesta_SpeciesMass(IdxSpecies(I))
                         Call Die('ReadInput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
                      END IF
                   ELSE
                      AtName(I)=AdjustL(Trim(Siesta_SpeciesName(IdxSpecies(I))))
                      Atome(I)=Siesta_SpeciesMass(IdxSpecies(I))
                   END IF
     ! we look for something else at the end of the line
     ! We save everything but the x,y,z, and species description

            subroutine Header(String)
           subroutine Header(String)
     ! This short subroutine print a header in a nice way
              IMPLICIT NONE
              CHARACTER(*) :: String
              INTERFACE
                 function valid(string) result (isValid)
                   CHARACTER(*), intent(in) :: string
                   logical                  :: isValid
                 END function VALID
              END INTERFACE
              CHARACTER(70) :: Head1,Sep
              CHARACTER(70) :: String2
              INTEGER :: LenS,Len1,Pos1
              INTEGER :: I,N,Idx,PosB,PosE
              LOGICAL :: Debug
              Debug=Valid("Header")
              Head1="====================================================================="
                Sep="                                                                     "
              LenS=len(String)
                If (debug) THEN
                   WRITe(*,'(A)') Head1
                   WRITE(*,'(A)') "=     Entering Header.... "
                   WRITe(*,'(A)') Head1
                END IF
              WRITe(*,'(A)') Head1
              LenS=len_Trim(String)
              DO WHILE (LenS>65)
                 Idx=Index(String(1:65)," ",BACK=.TRUE.)
                 If (Idx>0) THEN
                    PosE=Idx-1
                 ELSE
                    PosE=65
                 END IF
                 String2=AdjustL(Trim(String(1:PosE)))
                 LenS=Len_Trim(String2)
      !        WRITE(*,*) "DBG Head:",LenS, String
                 Pos1=33-Int(LenS/2)
                 Len1=67-LenS-Pos1
                 WRITE(*,'(A)') "=" // Sep(1:Pos1) // TRIM(String2) // Sep(1:Len1) // "="
                 String2=String(PosE+1:)
                 String=AdjustL(String2)
                 LenS=Len_Trim(String)
           END DO
              Pos1=34-Int(LenS/2)
              Len1=67-LenS-Pos1
              WRITE(*,'(A)') "=" // Sep(1:Pos1) // TRIM(String) // Sep(1:Len1) // "="
              WRITe(*,'(A)') Head1
              WRITE(*,'(A)') "=" // Sep(1:Pos1) // String(1:LenS) // Sep(1:Len1) // "="
              WRITe(*,'(A)') Head1
                If (debug) THEN
                   WRITe(*,'(A)') Head1
                   WRITE(*,'(A)') "=     Exiting Header.... "
                   WRITe(*,'(A)') Head1
                END IF
            END subroutine Header

       CASE (1)
          call getarg(1,FileIn)
          IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
             WRITE(*,*) "Path mini-help"
             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(*,*) ""
             WRITE(*,*) "Compulsory variables are:"
             WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
             WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
             WRITE(*,*) "Nat   : Number of atoms"
     ! 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(*,*) ""
             WRITE(*,*) ""
             WRITE(*,*) "Other options"
             WRITE(*,*) "Input: string that indicates the type of the input geometries."
             WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
             WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
             WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext"
             WRITE(*,*) "   *   In case of a Gaussian or Mopac calculations, input must be set to Cart."
             WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
             WRITE(*,*) "     input file. See example file Test_G03.path or Test_Mopac.path"
             WRITE(*,*) "   *   In the case of a VASP calculation, if input is set to Cart, then"
             WRITE(*,*) "     the preamble of a VASP calculation should be added at the end of"
             WRITE(*,*) "     the input file. See example file Test_VASP_cart.path"
             WRITE(*,*) "       In the case of a VASP calculation, one should also give value of the "
             WRITE(*,*) "    Runmode variable"
             WRITE(*,*) "   *   In the case of a SIESTA calculation, an example of a Siesta input file"
             WRITE(*,*) "     should be added at the end of the input file. See Test_Siesta.path"
             WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy"
             WRITE(*,*) "       and gradient of the whole path or not. If one wants to use VASP,"
             WRITE(*,*) "       Runmode must be set to PARA."
             WRITE(*,*) "PathName: Prefix used to save the path. Default is Path"
             WRITE(*,*) "Poscar: string that will be used as the prefix for the different "
             WRITE(*,*) "        POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.,"
             WRITE(*,*) "        not used for internal calculations."
             WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations"
             WRITE(*,*) " ISuffix: Suffix for the input file."
             WRITE(*,*) " OSuffix: suffix for the output file."
             WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates."
             WRITE(*,*) "      Valid only for Coord=Zmat, Hybrid or Mixed"
             WRITE(*,*) "Fact: REAL used to define if two atoms are linked."
             WRITE(*,*) "      If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked."
             WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
             WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
             WRITE(*,*) "       - CART (or Xyz): works in cartesian"
             WRITE(*,*) "       - Zmat: works in internal coordinates (Zmat)"
             WRITE(*,*) "       - Mixed: frozen atoms, as well as atoms defined by the "
             WRITE(*,*) "       'cart' array(see below) are describe in CARTESIAN, whereas the"
             WRITE(*,*) "       others are described in Zmat"
             WRITE(*,*) "       - Baker: use of Baker coordinates, also called delocalized internal coordinates"
             WRITE(*,*) "       - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
             WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:"
             WRITE(*,*) "       - RFO: Rational function optimization"
             WRITE(*,*) "       - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
     !        WRITE(*,*) " HesUpd: Deprecated. Use Hupdate."
             WRITE(*,*) " HUpdate: method to update the hessian. By default, it is Murtagh-Sargent"
             WRITE(*,*) "      Except for geometry optimization where it is BFGS."
             WRITE(*,*) " GeomFile: Name of the file to print the geometries and their energy"
             WRITE(*,*) "      during a geometry optimization. If this variable is not given"
             WRITE(*,*) "      then nothing is printed."
             WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
             WRITE(*,*) "Smax: Maximum length of a step during path optimization"
             WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
             WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
             WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path"
             WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
             WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0."
             WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
             WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
             WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
             WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, "
             WRITE(*,*) "        it might happen that an atom moves out of the unit cell."
             WRITE(*,*) "        Path detects this by comparing the displacement to boxtol:"
             WRITE(*,*) "        if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5"
             WRITE(*,*) ""
             WRITE(*,*) "Arrays:"
             WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements."
             WRITE(*,*) "      You can modify it using, rcov(6)=0.8."
             WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements."
             WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations."
             WRITe(*,*) "If not given, Path will read the POTCAR file."
             WRITE(*,*) ""
             WRITE(*,*) "Flags:"
             WRITE(*,*) "MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
             WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time."
             WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
             WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
             WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para"
             WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para"
             WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops."
             WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
             WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
             WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path."
             WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms."
             WRITE(*,*) "         If True, a &frozenlist namelist containing the list of frozen atoms must be given."
             WRITE(*,*) "          If VASP is used, and frozen is not given"
             WRITE(*,*) " here, the program will use the F flags of the input geometry"
             WRITE(*,*) "FCart:  True if one wants to describe some atoms using cartesian coordinates. "
             WRITE(*,*) "        *** Only used in 'mixed' calculations. ***"
             WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
             WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
             WRITE(*,*) ""
             WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis."
             WRITE(*,*) "        If True, Opt'n Path will look for the AnaList namelist after the Path Namelist."
             WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then"
             WRITE(*,*) "     Opt'n Path will read it and print the values of the variable in a .dat file"
             WRITE(*,*) "    If AnaGeom is T but Analist is not given, then Path.dat just contains the energy."
             WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description"
             WRITE(*,*) "of the variables among:"
             WRITE(*,*) "b(ond) At1 At2"
             WRITE(*,*) "a(ngle) At1 At2 At3"
             WRITE(*,*) "d(ihedral) At1 At2 At3 At4"
             WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added"
             WRITE(*,*) "at the end of the real atoms of the system"
             WRITE(*,*) ""
             WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
             WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
             WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
             WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
             WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now"
             WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files."
             WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP."
             WRITE(*,*) ""
             STOP
          END IF
-...
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      !
       ! For Debugging purposes
+      !
       ! OptGeom is the index of the geometry which is to be optimized.
       IF (Optgeom.GT.0) THEN
-...
          Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
          WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
          WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
          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
-...
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       IF (PathOnly) THEN
          WRITE(*,*) "PathOnly=.T. , Stop here"
          Call Header("PathOnly=.T. , Stop here")
          Call Write_path(-1)
          if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
          STOP

          END DO
          DEALLOCATE(GradTmp)
       CASE ("CART","HYBRID")
     ! PFL 2013 Feb
     ! We could think of just doing:
     !     Grad=GradCart
     ! however, this does not work as implicitly we have :
     ! Grad(Nat,3)  BUT GradCar(3,Nat)
     ! So we have to use reshape :-(
          Grad=reshape(reshape(GradCart,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))
          Grad=GradCart
       CASE DEFAULT
          WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP"
          STOP

       LOGICAL  :: debug
       LOGICAL :: Fini
       LOGICAL :: FiniS,FiniG,Fini
       LOGICAL, SAVE :: FRST=.TRUE.
       ! Variables
       INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
       REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
       REAL(KREAL), ALLOCATABLE :: Step(:) ! 3*Nat or 3*Nat-6
       REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
       REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6)
       REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6
-...
       REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6
       REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
       REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
       REAL(KREAL) :: NormStep, FactStep, HP, MaxStep
       REAL(KREAL) :: NormStep, FactStep, HP, MaxStep,NormGrad
       REAL(KREAL) :: Eold, Eini
       ! Values contains the values for the geometry analizes
       REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
       CHARACTER(LCHARS) :: Line
       debug=valid('OptGeom')
-...
       if (debug) Call Header("Entering Geom Optimization")
       ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
       ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord),Step(NCoord))
       ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
       ALLOCATE(GeomRef(NCoord))
       ALLOCATE(Hess_local(NCoord,NCoord))
-...
          IF (Step_method .EQ. "RFO") Then
             Hess_local=0.
             DO I=1,NCoord
                Hess_local(I,I)=0.5d0
                Hess_local(I,I)=1.d0
             END DO
          END IF
-...
          END IF
          !Print *, 'Hess_local after inversion='
          DO I=1,NCoord
     !     DO I=1,NCoord
             !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
          END DO
     !     END DO
          IF (Step_method .EQ. "RFO") Then
             Hess_local=0.
-...
       DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
          IOpt=IOpt+1
          Write(Line,'(1X,A,I5)') "Iteration ",Iopt
          Call Header(TRIM(Line))
          WRITE(IoOut,*) "Current Geometry"
          Line=""
          Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
          if (COORD/="CART") THEN
             WRITE(IoOut,*) "Current Geometry in Cartesian coordinates"
             Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)
          END IF
          WRITE(IoOut,*) "Computing energy and gradient"
          ! Calculate the  energy and gradient
          IF (debug) THEN
             WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
-...
          ELSE
             Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
          END IF
     !!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!
     ! We have a pb with the format of GradTmp(Nat,3) whereas it is 3,Nat in Egrad
+    !
     ! This is a path for CART coordinates !!!
          IF (COORD=="CART") THEN
             Step=reshape(reshape(GradTmp,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))
             GradTmp=Step
          END IF
     !!!!!!!!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!!!!!!!!
          WRITE(IoOut,'(1X," Energy E=",F15.6,A)') E,TRIM(UnitProg)
          WRITE(IoOut,*) "Gradient for current geometry"
          Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
          If (Iopt==1) EIni=E
          IF (debug) THEN
-...
          ! GradTmp becomes Step in Step_RFO_all.
          SELECT CASE (Step_method)
          CASE ('RFO')
             Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
             Call Step_RFO_all(NCoord,Step,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
             GradTmp=Step
          CASE ('GDIIS')
             Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
             ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
-...
          NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
          if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
          FactStep=min(1.d0,MaxStep/NormStep)
          Fini=(NormStep.LE.SThresh)
          NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
          if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
          Fini=(Fini.AND.(NormStep.LE.GThresh))
          FiniS=((NormStep*FactStep).LE.SThresh)
          NormGrad=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
          if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormGrad (Gradient)=', NormGrad, 'Gradient Threshold=', GThresh
          FiniG=(NormGrad.LE.GThresh)
          Fini=FiniS.AND.FiniG
          if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
          GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
          GradTmp=GradTmp*FactStep ! GradTmp is step  here, from Step_RFO_all.
          WRITE(IoOut,*) " Converged ?"
          WRITE(IoOut,'(1X,A18,1X,A15,1X,A13,5X,A3)') " Criterion","Current Value","Threshold","Ok?"
          IF (FiniS) THEN
             WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"YES"
          ELSE
             WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"NO"
          END IF
          IF (FiniG) THEN
             WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"YES"
          ELSE
             WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"NO"
          END IF
          if (DynMaxStep.AND.(IOpt>1)) THEN
             If (E<EOld) Then
-...
          Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)
          Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
     ! We add the step to the old geom
          Geom=GeomOld+GradTmp
          EOld=E
-...
       END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
       DEALLOCATE(GradTmp2,GeomTmp2,GradTmp)
       DEALLOCATE(GradOld,GeomOld)
       DEALLOCATE(Hess_local)
       DEALLOCATE(GeomCart_old)
       DEALLOCATE(NewGeom,NewGradTmp)
       DEALLOCATE(Hess_local_inv)
       IF (Fini) THEN
          WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
          WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles"
       ELSE
          WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E
          WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles"
       END IF
       Call Header(Line)
       WRITE(IoOut,*) "Last Energy E=",E
       WRITE(*,*) "Initial Geometry:"
       DO I=1, Nat
          WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
       END DO
       WRITE(*,*) "Final Geometry:"
       DO I=1, Nat
          WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
     !  WRITE(*,*) "Initial Geometry:"
     !  DO I=1, Nat
     !     WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
     !  END DO
       WRITE(IoOut,*) "Final Geometry:"
          Line=""
          Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
          if (COORD/="CART") THEN
             WRITE(IoOut,*) "Last Geometry in Cartesian coordinates"
             Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)
          END IF
     !  DO I=1, Nat
     !     WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
          !IF (I .GT. 1) Then
          !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
          !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
          !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
          !END IF
       END DO
     !  END DO
       IF (COORD .EQ. "BAKER") Then
          I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
-...
       END IF
       DEALLOCATE(GeomCart)
       DEALLOCATE(GradTmp2,GeomTmp2,GradTmp,Step)
       DEALLOCATE(GradOld,GeomOld)
       DEALLOCATE(Hess_local)
       DEALLOCATE(GeomCart_old)
       DEALLOCATE(NewGeom,NewGradTmp)
       DEALLOCATE(Hess_local_inv)
       if (debug) Call Header("Geom Optimization Over")

       CHARACTER(LCHARS) :: Siesta_Label, Siesta_CoordFile
     ! Number of species used in Siesta
       INTEGER(KINT) :: Siesta_NbSpecies
     ! Mass number for each species
       REAL(KREAL), ALLOCATABLE :: ListSpecies(:) ! NbSpecies
     ! Mass number for each species (atomic number)
       INTEGER(KINT), ALLOCATABLE :: Siesta_SpeciesMass(:) ! NbSpecies
     ! Name of each species
       CHARACTER(LCHARS), ALLOCATABLE :: Siesta_SpeciesName(:) ! NbSpecies
     ! Species for each atom

Formats disponibles : Unified diff

Chimie Théorique » Opt'n Path

Révision 9