Révision 9

doc/Mini_help.tex (revision 9)
1 1
\documentclass[a4paper,11pt]{article}
2
\usepackage{a4wide}
3 2
%\usepackage{times}
4 3
\usepackage{graphicx}
5 4
\usepackage[body={16cm,24.5cm}]{geometry}
......
34 33
\lhead[]{}
35 34
\rhead[]{}
36 35
\cfoot{}
37
\title{OpenPath mini-help}
36
\title{Opt'n Path mini-help}
38 37
\author{P. Fleurat-Lessard, P. Dayal \\
39
Laboratoire de Chimie de l'ENS Lyon, 46 all?e d'Italie, F-69364 Lyon
38
Laboratoire de Chimie de l'ENS Lyon, 46 all\'ee d'Italie, F-69364 Lyon
40 39
Cedex 7}
41
\date{Feb. 2011}
40
\date{March 2013}
42 41

  
43
\def\Path{\texttt{OpenPath}}
42
\def\Path{\texttt{Opt'n Path}}
44 43
\begin{document}
45 44
\maketitle
46 45

  
47
\section{introduction}
46
\section{Introduction}
48 47
\Path{} is a program that can optimize a reaction path between two
49 48
  structures. The algorithm to optimize the path is close to the string
50 49
  method. The originality of this program lies in the coordinate set it
......
52 51

  
53 52
This program is an independant program that calls standard electronic
54 53
  structure codes to get the energies and forces it needs to optimize
55
  the reaction path. For now, it is coupled to Gaussian, MOPAC2009,
56
  Vasp and Turbomole.
54
  the reaction path. For now, it is coupled to Gaussian, MOPAC,
55
  Vasp, Turbomole and Siesta.
57 56

  
58 57
\section{Installation}
59 58
We suppose here that you will install \Path{} into a directory called
60
\texttt{OpenPath}. First, create the directory: 
59
\texttt{optnpath}. First, create the directory: 
61 60
\begin{verbatim}
62
 mkdir OpenPath
63
 cd OpenPath
64
 mv ../OpenPath.tgz .
61
 mkdir optnpath
62
 cd optnpath
63
 mv ../optnpath_x.yy.tgz .
65 64
\end{verbatim}
66 65
Uncompress the archive:
67 66
\begin{verbatim}
68
 gunzip OpenPath.tgz
69
 tar -xvf OpenPath.tar
67
 gunzip optnpath_x.yy.tgz
68
 tar -xvf Optnpath.tar
70 69
\end{verbatim}
71 70
You should now have this \texttt{Mini\_help.pdf} file and 4 directories (doc, src, utils, examples).
72 71

  
......
133 132

  
134 133
\end{verbatim}
135 134

  
136
   Compulsory variables are: 
135
\subsubsection{Compulsory variables are:}
137 136
\begin{description}
138
\item          NGeomi: Number of geometries defining the Initial path 
139
\item          NGeomf: Number of geometries defining the Final path 
140
\item         Nat   : Number of atoms 
137
\item[NGeomi:] Number of geometries defining the Initial path 
138
\item[NGeomf:] Number of geometries defining the Final path 
139
\item[Nat:] Number of atoms 
141 140
\end{description}
142
           
143
          Other options 
141

  
142
          
143
\subsubsection{Other options:}
144 144
\begin{description}
145
\item  Input: string that indicates the type of the input geometries.
146
          Accepted values are: Cart (or Xmol or Xyz) or Vasp 
147
\item    Prog: string that indicates the program that will be used for energy and gradient calculations. 
148
               Accepted values are: Gaussian, Vasp or Ext. \\
149
                In case of a Gaussian calculations, input must be set to Cart. 
145
\item[Input:] String that indicates the type of the input geometries.
146
          Accepted values are: Cart (or Xmol or Xyz), Vasp, Turbomole or Siesta.
147
\item[Prog:] string that indicates the program that will be used for energy and gradient calculations. 
148
               Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext. \\
149
\begin{itemize}
150
\item   In case of a Gaussian calculations, input must be set to Cart. 
150 151
                One example of a gaussian input should be added at the end of the 
151 152
                input file.See example file \texttt{Test\_HCN\_zmat\_g03.path}. \\
152
                In the case of a VASP calculation, if input is set to Cart, then 
153
              the preamble of a VASP calculation should be added at the end of 
153
\item       In the case of a VASP calculation, if input is set to Cart, then 
154
              the preamble of a VASP calculation s\item \texttt{Mopac}: Examples  using sequential call
155
  to MOPAC2009 to calculate the energies and forces along the path.
156
hould be added at the end of 
154 157
              the input file. See example file \texttt{Test\_VASP\_cart.path}. 
155 158
              In the case of a VASP calculation, one should also give value of the  
156 159
              \texttt{RunMode} variable .
160
\item  In the case of a SIESTA calculation, an example of a Siesta input file
161
     should be added at the end of the input file. See \texttt{Test\_Siesta.path}.
162
\end{itemize}
157 163

  
158
\item RunMode: This indicates wether \Path{} should use VASP
159
              routine to calculate the energy and gradient of the whole path
164
\item[RunMode:] When running on a multi-processor machine, this indicates
165
  wether \Path{} should calculate the energy and gradient of the whole path in parallel
160 166
              or not. User has two options. One is to calculate the energy and
161 167
              gradient of each point sequentially. This is usefull when
162 168
              running on one (or two) processors. In this case,
......
164 170
              in parallel with 8 or more processors, one can use VASP to
165 171
              calculate simultaneously the energies and gradients for all
166 172
              points, as in a normal NEB calculation. In this case,
167
              \texttt{RunMode} must be set to \texttt{PARA}.
168
\item ProgExe: Name (with full path) of the executable to be used to get
173
              \texttt{RunMode} must be set to \texttt{PARA}. 
174
 \emph{For now, this is usefull only for VASP.} \\
175

  
176
\item[ProgExe:] Name (with full path) of the executable to be used to get
169 177
  energies and gradients. For example, if VASP is used in parallel, one might
170 178
  have something like: \\ 
171 179
\verb!ProgExe='/usr/local/mpich/bin/mpirun -machinefile machine -np 8 ~/bin/VASP_46'!.
172 180
Another option that I use, is to put \verb!ProgExe='./run_vasp'! and I put every option needed to run VASP
173 181
into the \texttt{run\_vasp} file. 
174
\item   PathName: Prefix used to save the path. Default is Path 
175
\item   Poscar: string that will be used as the prefix for the different  
182
\item[EReac:] (REAL) By default, \Path{} does not compute the energy of the reactants
183
 and products. This thus indicates the reactants energy to \Path{} to have better plots
184
 at the end.
185
\item[EProd:] (REAL) By default, \Path{} does not compute the energy of the reactants
186
 and products. This thus indicates the products energy to \Path{} to have better plots.
187
\item[PathName:] Prefix used to save the path. Default is Path 
188
\item[Poscar:] string that will be used as the prefix for the different  
176 189
           POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE., 
177 190
           not used for internal calculations. 
178
\item   IGeomRef: Index of the geometry used to construct the internal coordinates. 
179
         Valid only for Coord=Zmat, Hybrid or Mixed 
180
\item   Fact: REAL used to define if two atoms are linked. 
191
\item[CalcName:] Prefix for the files used for the energy and gradient calculations
192
\item[ISuffix:] Suffix for the input file used for energy and gradient calculations. 
193
The full inputfile name will thus be \textit{CalcName}.\textit{ISuffix}.
194
\item[OSuffix:] Suffix for the output file used for energy and gradient calculations.
195
 The full outputfile name will thus be \textit{CalcName}.\textit{OSuffix}.
196
\item[IGeomRef:] Index of the geometry used to construct the internal coordinates. 
197
         Valid only for Coord=Zmat, Hybrid or Mixed.
198
\item[Fact:] REAL used to define if two atoms are linked. 
181 199
         If $d(A,B) \leq fact*(rcov(A)+rcov(B))$, then A and B are considered Linked. 
182
\item   debugFile: Name of the file that indicates which subroutine should print debug info. 
183
\item   Coord: System of coordinates to use. Possible choices are: 
184
          - CART (or Xyz): works in cartesian 
185
          - Zmat: works in internal coordinates (Zmat) 
186
          - Mixed: frozen atoms, as well as atoms defined by the  
200
\item[debugFile:] Name of the file that indicates which subroutine should print debug info. 
201
\item[Coord:] System of coordinates to use. Possible choices are: 
202
\begin{itemize}
203
\item CART (or Xyz): works in cartesian 
204
\item Zmat: works in internal coordinates (Zmat) 
205
\item Mixed: frozen atoms, as well as atoms defined by the  
187 206
          'cart' array(see below) are describe in CARTESIAN, whereas the 
188 207
          others are described in Zmat  
189
          - Baker: use of Baker coordinates, also called delocalized internal coordinates 
190
          - Hybrid: geometries are described in zmat, but the gradient are used in cartesian 
191
\item   Step\_method: method to compute the step for optimizing a geometry; choices are:  
192
          - RFO: Rational function optimization 
193
          - GDIIS: Geometry optimization using direct inversion in the iterative supspace 
194
\item    HUpdate: method to update the hessian. By default, it is Murtagh-Sargent 
208
\item  Baker: use of Baker coordinates, also called delocalized internal coordinates 
209
\item Hybrid: geometries are described in zmat, but the gradient are used in cartesian 
210
\end{itemize}
211
\item[Step\_method:] method to compute the step for optimizing a geometry; choices are:  
212
\begin{itemize}
213
\item RFO: Rational function optimization 
214
\item GDIIS: Geometry optimization using direct inversion in the iterative supspace 
215
\end{itemize}
216
\item[HUpdate:] method to update the hessian. By default, it is Murtagh-Sargent 
195 217
         Except for geometry optimization where it is BFGS. 
196
\item   MaxCyc: maximum number of iterations for the path optimization 
197
\item   Smax: Maximum length of a step during path optimization 
198
\item   SThresh: Step Threshold to consider that the path is stationary 
199
\item   GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken 
200
\item   FTan: We moving the path, this gives the proportion of the displacement tangent to the path 
218
\item[MaxCyc:] maximum number of iterations for the path optimization 
219
\item[Smax:] Maximum length of a step during path optimization 
220
\item[SThresh:] Step Threshold to consider that the path is stationary 
221
\item[GThresh:] Gradient Threshold to consider that the path is stationary, only orthogonal part is taken 
222
\item[FTan:] We moving the path, this gives the proportion of the displacement tangent to the path 
201 223
         that is kept. FTan=1. corresponds to the full displacement,  
202 224
         whereas FTan=0. gives a displacement orthogonal to the path. 
203
\item   IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization. 
204
\item   ISpline: By default, a linear interpolation is used to generate the path. 
225
\item[IReparam:] The path is not reparameterised at each iteration. This gives the 
226
frequency of reparameterization.  
227
\item[IReparamT:] When the path is not reparameterised at each iteration, this gives the 
228
frequency of reparameterization of the \emph{tangents}.
229
\item[ISpline:] By default, a linear interpolation is used to generate the path. 
205 230
            This option indicates the first step where spline interpolation is used. 
231
\item[BoxTol:]  Real between 0. and 1. When doing periodic calculations, 
232
        it might happen that an atom moves out of the unit cell.
233
        Path detects this by comparing the displacement to boxtol:
234
        if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5.
235
\item[FrozTol:]  (Real) This indicates the threshold to determine wether an atom moves between two images. Default is 1e-4.
236
\item[OptGeom:] This INTEGER indicates the index of the geometry you want to optimize. 
237
If \texttt{OptGeom} is set, then \Path{} performs a geometry optimization instead of a path interpolation.
238
\item[GeomFile:]  Name of the file to print the geometries and their energy
239
      during a geometry optimization. If this variable is not given
240
      then nothing is printed.
241
\item[AnaFile:]  Name of the file to print the values of the geometrical parameters
242
that are monitored if \texttt{AnaGeom=.TRUE.}. Default is \textit{PathName}.dat
243
\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
244
%% Not described here: NMaxPtPath, NGintMax (too technical ?)
206 245
\end{description}
207 246
           
208
          Arrays: 
247
\subsubsection{Arrays:}
209 248
\begin{description}
210
\item   Rcov: Array containing the covalent radii of the first 80 elements. 
249
\item[Rcov:] Array containing the covalent radii of the first 86 elements. 
211 250
         You can modify it using, \verb!rcov(6)=0.8!. 
212
\item   Mass: Array containing the atomic mass of the first 80 elements. 
213
\item   AtTypes: Name of the different atoms used in a VASP calculations. 
251
\item[Mass:] Array containing the atomic mass of the first 86 elements. 
252
\item[AtTypes:] Name of the different atoms used in a VASP calculations. 
214 253
          If not given, Path will read the POTCAR file. 
215 254
\end{description}
216 255
           
217
          Flags: 
256
\subsubsection{Flags:}
218 257
\begin{description}
219
\item        MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE. 
220
\item       Renum: Flag. True if one wants to reoder the atoms in the
221
  initial order. default is .TRUE. most of the time.  
222
\item       OptProd: True if one wants to optimize the geometry of the products. 
223
\item       OptReac: True if one wants to optimize the geometry of the reactants. 
224
\item       CalcEProd: if TRUE the product energy will be
225
  computed. Default is FALSE. Not Compatible with RunMode=Para". 
226
\item       CalcEReac: if TRUE the reactants energy will be
227
  computed. Default is FALSE. Not Compatible with RunMode=Para". 
228
\item       PathOnly:TRUE if one wants to generate the initial path, and stops.  
229
\item       Hinv: if True, then Hessian inversed is used. 
230
\item       IniHup: if True, then Hessian inverse is extrapolated
258
\item[MW:]  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE. 
259
\item[Renum:] Flag. True if one wants to reoder the atoms in the
260
  initial order. default is .TRUE. unless for \texttt{Coord} equals \texttt{CART}.
261
\item[OptProd:] True if one wants to optimize the geometry of the products. 
262
\item[OptReac:] True if one wants to optimize the geometry of the reactants. 
263
\item[CalcEProd:] if TRUE the product energy will be
264
  computed. Default is FALSE. Not Compatible with \texttt{RunMode=Para}. 
265
\item[CalcEReac:] if TRUE the reactants energy will be
266
  computed. Default is FALSE. Not Compatible with \texttt{RunMode=Para}. 
267
\item[PathOnly:] TRUE if one wants to generate the initial path, and stops.  
268
\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.
269
\item[Hinv:] if True, then Hessian inversed is used. 
270
\item[IniHup:] if True, then Hessian inverse is extrapolated
231 271
  using the initial path calculations.  
232
\item       HupNeighbour: if True, then Hessian inverse is
272
\item[HupNeighbour:]  if True, then Hessian inverse is
233 273
  extrapolated using the neighbouring points of the path.  
234
\item       FFrozen: True if one wants to freeze the positions of some atoms.  
274
\item[FFrozen:] True if one wants to freeze the positions of some atoms.  
235 275
                If True, a \verb!&frozenlist! namelist containing the
236 276
                list of frozen atoms must be given.  
237 277
                 If VASP is used, and frozen is not given 
238 278
        here, the program will use the F flags of the input geometry 
239
\item       FCart:  True if one wants to describe some atoms using cartesian coordinates.  
279
\item[FCart:]  True if one wants to describe some atoms using cartesian coordinates.  
240 280
               *** Only used in 'mixed' calculations. *** 
241 281
             If True, a \verb!&cartlist! namelist containing the list of cart atoms must be given. 
242 282
             By default, only frozen atoms are described in cartesian coordinates. 
243 283
        
244
\item       Autocart: True if you want to let the program choosing the cartesian atoms. 
245
\item       VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now 
284
\item[Autocart:] True if you want to let the program choosing the cartesian atoms. 
285
\item[VMD:] TRUE if you want to use VMD to look at the Path. Used only for VASP for now.
286
\item[WriteVASP:] TRUE if you want to print the images coordinates in POSCAR files.
287
See also the POSCAR option. This can be used only if prog or input=VASP.
288
\item[AnaGeom:] If TRUE, Opt'n Path will create a file .dat for geometries analysis.
289
        If True, Opt'n Path will look for the AnaList namelist after the Path Namelist.
290
\item[DynMaxStep:] if TRUE, the maximum allowed step is updated at each step, for each geometry.
291
        If energy goes up, $Smax=Smax*0.8$, if not $Smax=Smax*1.2$. 
292
       It is ensured that the dynamical Smax is within $[0.5*SMax_0,2*Smax_0]$.
246 293
\end{description}
247 294

  
295
\subsubsection{Additional namelists:}
296
\begin{description}
297
\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: \\
298
\texttt{\&cartlist list = 1 2 6 12 -20 \&end} \\
299
will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian.
300
\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. 
301
\item[\&Analist nb= \ldots{} \&end] list of variables for geometry analysis. If present and if AnaGeom=T then
302
     Opt'n Path will read it and print the values of the variable in a .dat file
303
    If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. 
304
Analist contains the number of variables to monitor: Nb, and is followed by the description
305
of the variables among:
306
b(ond) At1 At2
307
a(ngle) At1 At2 At3
308
d(ihedral) At1 At2 At3 At4
309
c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added
310
at the end of the real atoms of the system.
311
\end{description}
312

  
313
\subsection{Examples}
248 314
More to come...      have a look at the files provided in the
249 315
\texttt{examples} directory. In particular, you will find there three
250 316
directories containing easy-to-use examples. That is, each directory contains
......
289 355
\end{itemize}
290 356
\end{itemize}
291 357
\item \texttt{Mopac}: Examples  using sequential call
292
  to MOPAC2009 to calculate the energies and forces along the path.
358
  to MOPAC to calculate the energies and forces along the path.
359
\item \texttt{Siesta}: Examples  using sequential call
360
  to Siesta to calculate the energies and forces along the path.
293 361
\item \texttt{Test}: Examples  using the analytical HCN potential
294 362
  energy surface to calculate the energies and forces along the path.
295 363
As this is fast, we also provide other Analysis tools.
src/EgradPath.f90 (revision 9)
696 696
        END IF
697 697

  
698 698
        GeomCart=0.d0
699
        IF(IOpt .EQ. 1) Then
700
           Print *, 'GeomTmp='
701
           WRITE(*,'(12(1X,F6.3))') GeomTmp
702
           Print *, 'GeomCart='
703
           WRITE(*,'(12(1X,F6.3))') GeomCart
704
      END IF
699
!        IF(IOpt .EQ. 1) Then
700
!           Print *, 'GeomTmp='
701
!           WRITE(*,'(12(1X,F6.3))') GeomTmp
702
!           Print *, 'GeomCart='
703
!           WRITE(*,'(12(1X,F6.3))') GeomCart
704
!      END IF
705 705
    
706 706
       ! GradTmp is gradient and calculated in Egrad.F, INTENT(OUT).
707 707
       ! GeomCart has INTENT(OUT)
src/ReadInput_siesta.f90 (revision 9)
179 179
           Idx=Index(Line," ")
180 180
           Line2=Trim(AdjustL(Line(Idx+1:)))
181 181
           Read(Line2,*) Siesta_NbSpecies
182
           ALLOCATE(ListSpecies(Siesta_NbSpecies))
182
           ALLOCATE(Siesta_SpeciesMass(Siesta_NbSpecies))
183 183
           ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies))
184 184
        END IF
185 185

  
......
206 206
              Read(Line,*) ITmp
207 207
              Idx=Index(Line,' ')
208 208
              Line=AdjustL(Line(Idx+1:))
209
              Read(Line,*) Ztmp
210
              ListSpecies(ITmp)=ZTmp
209
              Read(Line,*) Idx
210
              Siesta_SpeciesMass(ITmp)=Idx
211 211
              Idx=Index(Line,' ')
212 212
              Line=AdjustL(Line(Idx+1:))
213 213
              Idx=Index(Line,' ')
......
239 239
                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))
240 240
                    Call Die('Readinput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
241 241
                 END IF
242
! We look for the Mass number also
243
                 IF (Atome(I)/=Siesta_SpeciesMass(IdxSpecies(I))) THEN
244
                    Write(Line,'(A,I5,1X,I5,1X,I5)') 'AtMass /= SpeciesMass for atom ',I,Atome(I), &
245
                         Siesta_SpeciesMass(IdxSpecies(I))
246
                    Call Die('ReadInput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)
247
                 END IF
242 248
              ELSE
243 249
                 AtName(I)=AdjustL(Trim(Siesta_SpeciesName(IdxSpecies(I))))
250
                 Atome(I)=Siesta_SpeciesMass(IdxSpecies(I))
244 251
              END IF
245 252
! we look for something else at the end of the line
246 253
! We save everything but the x,y,z, and species description              
src/Header.f90 (revision 9)
1
       subroutine Header(String)
1
      subroutine Header(String)
2 2
! This short subroutine print a header in a nice way
3 3

  
4 4
         IMPLICIT NONE
5 5

  
6

  
6 7
         CHARACTER(*) :: String
8

  
9
         INTERFACE
10
            function valid(string) result (isValid)
11
              CHARACTER(*), intent(in) :: string
12
              logical                  :: isValid
13
            END function VALID
14
         END INTERFACE
15

  
16

  
7 17
         CHARACTER(70) :: Head1,Sep
18
         CHARACTER(70) :: String2
8 19
         INTEGER :: LenS,Len1,Pos1
20
         INTEGER :: I,N,Idx,PosB,PosE
21
         LOGICAL :: Debug
9 22

  
23

  
24
         Debug=Valid("Header")
25

  
10 26
         Head1="====================================================================="
11 27
           Sep="                                                                     "
12 28

  
13
         LenS=len(String)
29
           If (debug) THEN
30
              WRITe(*,'(A)') Head1
31
              WRITE(*,'(A)') "=     Entering Header.... "
32
              WRITe(*,'(A)') Head1
33
           END IF
34

  
35
         WRITe(*,'(A)') Head1
36

  
37
         LenS=len_Trim(String)
38
         DO WHILE (LenS>65)
39
            Idx=Index(String(1:65)," ",BACK=.TRUE.)
40
            If (Idx>0) THEN
41
               PosE=Idx-1
42
            ELSE 
43
               PosE=65
44
            END IF
45
            
46
            String2=AdjustL(Trim(String(1:PosE)))
47
            LenS=Len_Trim(String2)
14 48
 !        WRITE(*,*) "DBG Head:",LenS, String
49
            Pos1=33-Int(LenS/2)
50
            Len1=67-LenS-Pos1
51

  
52
            WRITE(*,'(A)') "=" // Sep(1:Pos1) // TRIM(String2) // Sep(1:Len1) // "="
53
            String2=String(PosE+1:)
54
            String=AdjustL(String2)
55
            LenS=Len_Trim(String)
56
      END DO
57

  
15 58
         Pos1=34-Int(LenS/2)
16 59
         Len1=67-LenS-Pos1
60
         
61
         WRITE(*,'(A)') "=" // Sep(1:Pos1) // TRIM(String) // Sep(1:Len1) // "="
17 62

  
18 63
         WRITe(*,'(A)') Head1
19
         WRITE(*,'(A)') "=" // Sep(1:Pos1) // String(1:LenS) // Sep(1:Len1) // "="
20
         WRITe(*,'(A)') Head1
21 64

  
65
           If (debug) THEN
66
              WRITe(*,'(A)') Head1
67
              WRITE(*,'(A)') "=     Exiting Header.... "
68
              WRITe(*,'(A)') Head1
69
           END IF
70

  
71

  
22 72
       END subroutine Header
src/Path.f90 (revision 9)
409 409
  CASE (1)
410 410
     call getarg(1,FileIn)
411 411
     IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
412
        WRITE(*,*) "Path mini-help"
412
        WRITE(*,*) "Opt'n Path mini-help"
413 413
        WRITE(*,*) "--------------"
414 414
        WRITE(*,*) ""
415 415
        WRITE(*,*) "Use: Path Input_file Output_file"
416 416
        WRITE(*,*) "Input_file starts with a Namelist called path"
417 417
        WRITE(*,*) ""
418
        WRITE(*,*) "Compulsory variables are:"
419
        WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
420
        WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
421
        WRITE(*,*) "Nat   : Number of atoms"
418
! The following has been generated (March 19, 2013) using the Mini_help.tex file
419
! 1) Edit Mini_help.tex
420
! 2) latex Mini_help.tex
421
! 3) catdvi -e 1 -U  Mini_help.dvi | sed -re "s/\[U\+2022\]/-/g" | sed -re "s/([^^[:space:]])\s+/\1 /g" > file.txt
422
! 4) edit file.txt to keep only the following lines, and to reformat a bit
423
! 5) awk '{print "WRITE(*,*) \"" $0 "\""}' file.txt > help
424
! 5) cut&paste help in Path.f90 ...
425
   WRITE(*,*) " &path"
426
   WRITE(*,*) " nat=3, ! Number of atoms"
427
   WRITE(*,*) " ngeomi=3, ! Number of initial geometries"
428
   WRITE(*,*) " ngeomf=12, !Number of geometries along the path"
429
   WRITE(*,*) " OptReac=.T., ! Do you want to optimize the reactants ?"
430
   WRITE(*,*) " OptProd=.T., ! Optimize the products"
431
   WRITE(*,*) " coord='zmat', ! We use Z-matrix coordinates"
432
   WRITE(*,*) " maxcyc=31, ! Max number of iterations"
433
   WRITE(*,*) " IReparam=2,! re-distribution of points along the path every 2 iterations"
434
   WRITE(*,*) " ISpline=50, ! Start using spline interpolation at iteration 50"
435
   WRITE(*,*) " Hinv=.T. , ! Use inverse of the Hessian internally (default: T)"
436
   WRITE(*,*) " MW=T, ! Works in Mass Weighted coordiante (default T)"
437
   WRITE(*,*) " PathName='Path_HCN_zmat_test', ! Name of the file used for path outputs"
438
   WRITE(*,*) " prog='gaussian',! we use G03 to get energy and gradients"
439
   WRITE(*,*) " SMax=0.1 ! Displacement cannot exceed 0.1 atomic units (or mass weighted at. unit)"
440
   WRITE(*,*) " /"
441
   WRITE(*,*) "  3"
442
   WRITE(*,*) " Energy : 0.04937364"
443
   WRITE(*,*) " H 0.0000 0.0000 0.0340"
444
   WRITE(*,*) " C 0.0000 0.0000 1.1030"
445
   WRITE(*,*) " N 0.0000 0.0000 2.2631"
446
   WRITE(*,*) "  3"
447
   WRITE(*,*) " Energy : 0.04937364"
448
   WRITE(*,*) " H 0.0000 1.1000 1.1030"
449
   WRITE(*,*) " C 0.0000 0.0000 1.1030"
450
   WRITE(*,*) " N 0.0000 0.0000 2.2631"
451
   WRITE(*,*) "3"
452
   WRITE(*,*) " CNH"
453
   WRITE(*,*) "H 0.000000 0.000000 3.3"
454
   WRITE(*,*) "C 0.000000 0.000000 1.1"
455
   WRITE(*,*) "N 0.000000 0.000000 2.26"
456
   WRITE(*,*) "%chk=Test"
457
   WRITE(*,*) "#P AM1 FORCE"
458
   WRITE(*,*) " HCN est bien"
459
   WRITE(*,*) ""
460
   WRITE(*,*) "0,1"
461
   WRITE(*,*) "H 0.000000 0.000000 0.000000"
462
   WRITE(*,*) "C 0.000000 0.000000 1.000"
463
   WRITE(*,*) "N 0.000000 0.000000 3.00"
464
   WRITE(*,*) ""
465
   WRITE(*,*) "*******"
466
   WRITE(*,*) "* Compulsory variables are:"
467
   WRITE(*,*) "*******"
468
   WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
469
   WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
470
   WRITE(*,*) "Nat: Number of atoms"
471
   WRITE(*,*) ""
472
   WRITE(*,*) "*******"
473
   WRITE(*,*) "* Other options:"
474
   WRITE(*,*) "*******"
475
   WRITE(*,*) "Input: String that indicates the type of the input geometries. Accepted values are: Cart (or"
476
   WRITE(*,*) "    Xmol or Xyz), Vasp, Turbomole or Siesta."
477
   WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
478
   WRITE(*,*) "    Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext."
479
   WRITE(*,*) ""
480
   WRITE(*,*) ""
481
   WRITE(*,*) "        - In case of a Gaussian calculations, input must be set to Cart. One example of a gaus-"
482
   WRITE(*,*) "          sian input should be added at the end of the input file.See example file Test_HCN_zmat_g03.path."
483
   WRITE(*,*) "        - In the case of a VASP calculation, if input is set to Cart, then the preamble of a"
484
   WRITE(*,*) "          VASP calculation s"
485
   WRITE(*,*) "        - Mopac: Examples using sequential call to MOPAC2009 to calculate the energies and"
486
   WRITE(*,*) "          forces along the path. hould be added at the end of the input file. See example file"
487
   WRITE(*,*) "          Test_VASP_cart.path. In the case of a VASP calculation, one should also give value"
488
   WRITE(*,*) "          of the RunMode variable ."
489
   WRITE(*,*) "        - In the case of a SIESTA calculation, an example of a Siesta input file should be added"
490
   WRITE(*,*) "          at the end of the input file. See Test_Siesta.path."
491
   WRITE(*,*) ""
492
   WRITE(*,*) "RunMode: When running on a multi-processor machine, this indicates wether Opt'n Path"
493
   WRITE(*,*) "    should calculate the energy and gradient of the whole path in parallel or not. User has"
494
   WRITE(*,*) "    two options. One is to calculate the energy and gradient of each point sequentially. This"
495
   WRITE(*,*) "    is usefull when running on one (or two) processors. In this case, RunMode should be put"
496
   WRITE(*,*) "    to SERIAL.When running in parallel with 8 or more processors, one can use VASP to"
497
   WRITE(*,*) "    calculate simultaneously the energies and gradients for all points, as in a normal NEB cal-"
498
   WRITE(*,*) "    culation. In this case, RunMode must be set to PARA. For now, this is usefull only for VASP."
499
   WRITE(*,*) ""
500
   WRITE(*,*) "ProgExe: Name (with full path) of the executable to be used to get energies and gradients."
501
   WRITE(*,*) "    For example, if VASP is used in parallel, one might have something like:"
502
   WRITE(*,*) "    ProgExe='/usr/local/mpich/bin/mpirun -machinefile machine -np 8 ~/bin/VASP_46'."
503
   WRITE(*,*) "    Another option that I use, is to put ProgExe='./run_vasp' and I put every option needed"
504
   WRITE(*,*) "    to run VASP into the run_vasp file."
505
   WRITE(*,*) "EReac: (REAL) By default, Opt'n Path does not compute the energy of the reactants and"
506
   WRITE(*,*) "    products. This thus indicates the reactants energy to Opt'n Path to have better plots at"
507
   WRITE(*,*) "    the end."
508
   WRITE(*,*) "EProd: (REAL) By default, Opt'n Path does not compute the energy of the reactants and"
509
   WRITE(*,*) "    products. This thus indicates the products energy to Opt'n Path to have better plots."
510
   WRITE(*,*) ""
511
   WRITE(*,*) "PathName: Prefix used to save the path. Default is Path."
512
   WRITE(*,*) "Poscar: string that will be used as the prefix for the different POSCAR files in a VASP calcu-"
513
   WRITE(*,*) "    lations. Usefull only if PathOnly=.TRUE., not used for internal calculations."
514
   WRITE(*,*) "CalcName: Prefix for the files used for the energy and gradient calculations"
515
   WRITE(*,*) ""
516
   WRITE(*,*) "ISuffix: Suffix for the input file used for energy and gradient calculations. The full inputfile"
517
   WRITE(*,*) "    name will thus be CalcName.ISuffix."
518
   WRITE(*,*) "OSuffix: Suffix for the output file used for energy and gradient calculations. The full outputfile"
519
   WRITE(*,*) "    name will thus be CalcName.OSuffix."
520
   WRITE(*,*) ""
521
   WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates. Valid only for"
522
   WRITE(*,*) "     Coord=Zmat, Hybrid or Mixed."
523
   WRITE(*,*) "Fact: REAL used to define if two atoms are linked. If d(A,B) =< fact* (rcov(A) + rcov(B)),"
524
   WRITE(*,*) "     then A and B are considered Linked."
525
   WRITE(*,*) ""
526
   WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
527
   WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
528
   WRITE(*,*) "        - CART (or Xyz): works in cartesian"
529
   WRITE(*,*) "        - Zmat: works in internal coordinates (Zmat)"
530
   WRITE(*,*) "        - Mixed: frozen atoms, as well as atoms defined by the 'cart' array(see below) are"
531
   WRITE(*,*) "          describe in CARTESIAN, whereas the others are described in Zmat"
532
   WRITE(*,*) "        - Baker: use of Baker coordinates, also called delocalized internal coordinates"
533
   WRITE(*,*) "        - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
534
   WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:"
535
   WRITE(*,*) ""
536
   WRITE(*,*) "        - RFO: Rational function optimization"
537
   WRITE(*,*) "        - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
538
   WRITE(*,*) "HUpdate: method to update the hessian. By default, it is Murtagh-Sargent Except for geom-"
539
   WRITE(*,*) "     etry optimization where it is BFGS."
540
   WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
541
   WRITE(*,*) ""
542
   WRITE(*,*) "Smax: Maximum length of a step during path optimization"
543
   WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
544
   WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is"
545
   WRITE(*,*) "     taken"
546
   WRITE(*,*) ""
547
   WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
548
   WRITE(*,*) "     that is kept. FTan=1. corresponds to the full displacement, whereas FTan=0. gives a"
549
   WRITE(*,*) "     displacement orthogonal to the path."
550
   WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of"
551
   WRITE(*,*) "     reparameterization."
552
   WRITE(*,*) "IReparamT: When the path is not reparameterised at each iteration, this gives the frequency"
553
   WRITE(*,*) "     of reparameterization of the tangents."
554
   WRITE(*,*) ""
555
   WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path. This option indicates"
556
   WRITE(*,*) "     the first step where spline interpolation is used."
557
   WRITE(*,*) "BoxTol: Real between 0. and 1. When doing periodic calculations, it might happen that an"
558
   WRITE(*,*) "     atom moves out of the unit cell. Path detects this by comparing the displacement to"
559
   WRITE(*,*) "     boxtol: if an atom moves by more than Boxtol, then it is moved back into the unit cell."
560
   WRITE(*,*) "     Default value: 0.5."
561
   WRITE(*,*) "FrozTol: (Real) This indicates the threshold to determine wether an atom moves between two"
562
   WRITE(*,*) "     images. Default is 1e-4."
563
   WRITE(*,*) ""
564
   WRITE(*,*) "OptGeom: This INTEGER indicates the index of the geometry you want to optimize. If"
565
   WRITE(*,*) "    OptGeom is set, then Opt'n Path performs a geometry optimization instead of a path"
566
   WRITE(*,*) "    interpolation."
567
   WRITE(*,*) "GeomFile: Name of the file to print the geometries and their energy during a geometry opti-"
568
   WRITE(*,*) "    mization. If this variable is not given then nothing is printed."
569
   WRITE(*,*) "AnaFile: Name of the file to print the values of the geometrical parameters that are monitored"
570
   WRITE(*,*) "    if AnaGeom=.TRUE.. Default is PathName.dat"
571
   WRITE(*,*) ""
572
   WRITE(*,*) "GplotFile: Name of the gnuplot file to plot the evolution of the geometrical parameters that"
573
   WRITE(*,*) "    are monitored if AnaGeom=.TRUE.. Default is PathName.gplot"
574
   WRITE(*,*) ""
575
   WRITE(*,*) "*******"
576
   WRITE(*,*) "* Arrays:"
577
   WRITE(*,*) "*******"
578
   WRITE(*,*) ""
579
   WRITE(*,*) "Rcov: Array containing the covalent radii of the first 86 elements. You can modify it using,"
580
   WRITE(*,*) "    rcov(6)=0.8."
581
   WRITE(*,*) "Mass: Array containing the atomic mass of the first 86 elements."
582
   WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations. If not given, Opt'n Path will"
583
   WRITE(*,*) "    read the POTCAR file."
584
   WRITE(*,*) ""
585
   WRITE(*,*) "*******"
586
   WRITE(*,*) "* Flags:"
587
   WRITE(*,*) "*******"
588
   WRITE(*,*) ""
589
   WRITE(*,*) "MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
590
   WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE."
591
   WRITE(*,*) "    unless for Coord equals CART."
592
   WRITE(*,*) ""
593
   WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
594
   WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
595
   WRITE(*,*) "CalcEProd: if TRUE the product energy will be computed. Default is FALSE. Not Compatible"
596
   WRITE(*,*) "    with RunMode=Para."
597
   WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compat-"
598
   WRITE(*,*) "    ible with RunMode=Para."
599
   WRITE(*,*) ""
600
   WRITE(*,*) "PathOnly: TRUE if one wants to generate the initial path, and stops."
601
   WRITE(*,*) "Align: If .FALSE., successive geometries along the path are not aligned on each other before"
602
   WRITE(*,*) "    path interpolation. Default is .TRUE. if there are 4 atoms or more."
603
   WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
604
   WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
605
   WRITE(*,*) ""
606
   WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of"
607
   WRITE(*,*) "    the path."
608
   WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms. If True, a &frozenlist"
609
   WRITE(*,*) "    namelist containing the list of frozen atoms must be given. If VASP is used, and frozen is"
610
   WRITE(*,*) "    not given here, the program will use the F flags of the input geometry"
611
   WRITE(*,*) "FCart: True if one wants to describe some atoms using cartesian coordinates. *** Only used in"
612
   WRITE(*,*) "    'mixed' calculations. *** If True, a &cartlist namelist containing the list of cart atoms"
613
   WRITE(*,*) "    must be given. By default, only frozen atoms are described in cartesian coordinates."
614
   WRITE(*,*) ""
615
   WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
616
   WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now."
617
   WRITE(*,*) ""
618
   WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files. See also"
619
   WRITE(*,*) "      the POSCAR option. This can be used only if prog or input=VASP."
620
   WRITE(*,*) "AnaGeom: If TRUE, Opt'n Path will create a file .dat for geometries analysis. If True, "
621
   WRITE(*,*) "Opt'n Path will look for the AnaList namelist after the Path Namelist."
622
   WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
623
   WRITE(*,*) "      If energy goes up, Smax = Smax*0.8, if not Smax = Smax*1.2. It is ensured that the"
624
   WRITE(*,*) "      dynamical Smax is within [0.5*SMax_0 ,2*Smax_0 ]."
625
   WRITE(*,*) ""
626
   WRITE(*,*) "*******"
627
   WRITE(*,*) "* Additional namelists:"
628
   WRITE(*,*) "*******"
629
   WRITE(*,*) ""
630
   WRITE(*,*) "&cartlist list= ... &end This gives the list of atoms that are described using cartesian coor-"
631
   WRITE(*,*) "      dinates. Read only if FCart=.TRUE.. To indicate an atom range, from 1 to 5 for example,"
632
   WRITE(*,*) "      one can put 1 -5 in the list. For example:"
633
   WRITE(*,*) "      &cartlist list = 1 2 6 12 -20 &end"
634
   WRITE(*,*) "      will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian."
635
   WRITE(*,*) "&Frozenlist list= ... &end This gives the list of atoms that are frozen during optimization."
636
   WRITE(*,*) "      Read only if FFrozen=.TRUE.. To indicate an atom range, from 1 to 5 for example, one"
637
   WRITE(*,*) "      can put 1 -5 in the list."
638
   WRITE(*,*) ""
639
   WRITE(*,*) "&Analist nb= ... &end list of variables for geometry analysis. If present and if Ana-"
640
   WRITE(*,*) "      Geom=T then Opt'n Path will read it and print the values of the variable in a .dat file If"
641
   WRITE(*,*) "      AnaGeom is T but Analist is not given, then Path.dat just contains the energy. Analist"
642
   WRITE(*,*) "      contains the number of variables to monitor: Nb, and is followed by the description of the"
643
   WRITE(*,*) "      variables among: b(ond) At1 At2 a(ngle) At1 At2 At3 d(ihedral) At1 At2 At3 At4 c NbAt"
644
   WRITE(*,*) "      At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added at the"
645
   WRITE(*,*) "      end of the real atoms of the system."
646
!!!!!!!!!! end of included file
422 647
        WRITE(*,*) ""
423
        WRITE(*,*) ""
424
        WRITE(*,*) "Other options"
425
        WRITE(*,*) "Input: string that indicates the type of the input geometries."
426
        WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
427
        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
428
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext"
429
        WRITE(*,*) "   *   In case of a Gaussian or Mopac calculations, input must be set to Cart."
430
        WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
431
        WRITE(*,*) "     input file. See example file Test_G03.path or Test_Mopac.path"
432
        WRITE(*,*) "   *   In the case of a VASP calculation, if input is set to Cart, then"
433
        WRITE(*,*) "     the preamble of a VASP calculation should be added at the end of"
434
        WRITE(*,*) "     the input file. See example file Test_VASP_cart.path"
435
        WRITE(*,*) "       In the case of a VASP calculation, one should also give value of the "
436
        WRITE(*,*) "    Runmode variable"
437
        WRITE(*,*) "   *   In the case of a SIESTA calculation, an example of a Siesta input file"
438
        WRITE(*,*) "     should be added at the end of the input file. See Test_Siesta.path"
439
        WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy"
440
        WRITE(*,*) "       and gradient of the whole path or not. If one wants to use VASP,"
441
        WRITE(*,*) "       Runmode must be set to PARA."
442
        WRITE(*,*) "PathName: Prefix used to save the path. Default is Path"
443
        WRITE(*,*) "Poscar: string that will be used as the prefix for the different "
444
        WRITE(*,*) "        POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.,"
445
        WRITE(*,*) "        not used for internal calculations."
446
        WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations"
447
        WRITE(*,*) " ISuffix: Suffix for the input file."
448
        WRITE(*,*) " OSuffix: suffix for the output file."
449
        WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates."
450
        WRITE(*,*) "      Valid only for Coord=Zmat, Hybrid or Mixed"
451
        WRITE(*,*) "Fact: REAL used to define if two atoms are linked."
452
        WRITE(*,*) "      If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked."
453
        WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
454
        WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
455
        WRITE(*,*) "       - CART (or Xyz): works in cartesian"
456
        WRITE(*,*) "       - Zmat: works in internal coordinates (Zmat)"
457
        WRITE(*,*) "       - Mixed: frozen atoms, as well as atoms defined by the "
458
        WRITE(*,*) "       'cart' array(see below) are describe in CARTESIAN, whereas the"
459
        WRITE(*,*) "       others are described in Zmat" 
460
        WRITE(*,*) "       - Baker: use of Baker coordinates, also called delocalized internal coordinates"
461
        WRITE(*,*) "       - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
462
        WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
463
        WRITE(*,*) "       - RFO: Rational function optimization"
464
        WRITE(*,*) "       - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
465
!        WRITE(*,*) " HesUpd: Deprecated. Use Hupdate."
466
        WRITE(*,*) " HUpdate: method to update the hessian. By default, it is Murtagh-Sargent"
467
        WRITE(*,*) "      Except for geometry optimization where it is BFGS."
468
        WRITE(*,*) " GeomFile: Name of the file to print the geometries and their energy"
469
        WRITE(*,*) "      during a geometry optimization. If this variable is not given"
470
        WRITE(*,*) "      then nothing is printed."
471
        WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
472
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
473
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
474
        WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
475
        WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path"
476
        WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
477
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0."
478
        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
479
        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
480
        WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
481
        WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, "
482
        WRITE(*,*) "        it might happen that an atom moves out of the unit cell."
483
        WRITE(*,*) "        Path detects this by comparing the displacement to boxtol:"
484
        WRITE(*,*) "        if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5"
485
        WRITE(*,*) ""
486
        WRITE(*,*) "Arrays:"
487
        WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements."
488
        WRITE(*,*) "      You can modify it using, rcov(6)=0.8."
489
        WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements."
490
        WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations."
491
        WRITe(*,*) "If not given, Path will read the POTCAR file."
492
        WRITE(*,*) ""
493
        WRITE(*,*) "Flags:"
494
        WRITE(*,*) "MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
495
        WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time."
496
        WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
497
        WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
498
        WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para"
499
        WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para"
500
        WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
501
        WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
502
        WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
503
        WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path."
504
        WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
505
        WRITE(*,*) "         If True, a &frozenlist namelist containing the list of frozen atoms must be given."
506
        WRITE(*,*) "          If VASP is used, and frozen is not given"
507
        WRITE(*,*) " here, the program will use the F flags of the input geometry"
508
        WRITE(*,*) "FCart:  True if one wants to describe some atoms using cartesian coordinates. "
509
        WRITE(*,*) "        *** Only used in 'mixed' calculations. ***"
510
        WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
511
        WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
512
        WRITE(*,*) ""
513
        WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis."
514
        WRITE(*,*) "        If True, Opt'n Path will look for the AnaList namelist after the Path Namelist."
515
        WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then"
516
        WRITE(*,*) "     Opt'n Path will read it and print the values of the variable in a .dat file"
517
        WRITE(*,*) "    If AnaGeom is T but Analist is not given, then Path.dat just contains the energy." 
518
        WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description"
519
        WRITE(*,*) "of the variables among:"
520
        WRITE(*,*) "b(ond) At1 At2"
521
        WRITE(*,*) "a(ngle) At1 At2 At3"
522
        WRITE(*,*) "d(ihedral) At1 At2 At3 At4"
523
        WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added"
524
        WRITE(*,*) "at the end of the real atoms of the system"
525
        WRITE(*,*) ""
526
        WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
527
        WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
528
        WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
529
        WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
530
        WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now"
531
        WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files."
532
        WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP."
533
        WRITE(*,*) ""
534 648
        STOP
535 649
     END IF
536 650

  
......
1223 1337

  
1224 1338
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1225 1339
  !
1226
  ! For Debugging purposes
1227 1340
  !
1228 1341
  ! OptGeom is the index of the geometry which is to be optimized.
1229 1342
  IF (Optgeom.GT.0) THEN
......
1252 1365
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1253 1366

  
1254 1367
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1255
     WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1368
     if (debug) WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1256 1369
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1257 1370
     STOP
1258 1371
  END IF ! IF (Optgeom.GT.0) THEN
......
1262 1375
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1263 1376

  
1264 1377
  IF (PathOnly) THEN
1265
     WRITE(*,*) "PathOnly=.T. , Stop here"
1378
     Call Header("PathOnly=.T. , Stop here")
1266 1379
     Call Write_path(-1)
1267 1380
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1268 1381
     STOP
src/Egrad.f90 (revision 9)
454 454
     END DO
455 455
     DEALLOCATE(GradTmp)
456 456
  CASE ("CART","HYBRID")
457
! PFL 2013 Feb 
458
! We could think of just doing:
459
!     Grad=GradCart
460
! however, this does not work as implicitly we have :
461
! Grad(Nat,3)  BUT GradCar(3,Nat)
462
! So we have to use reshape :-(
463
     Grad=reshape(reshape(GradCart,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))
457
     Grad=GradCart
464 458
  CASE DEFAULT
465 459
     WRITE(*,*) "Coord=",AdjustL(Trim(COORD))," not yet implemented in Egradients. STOP"
466 460
     STOP
src/Opt_Geom.f90 (revision 9)
73 73

  
74 74

  
75 75
  LOGICAL  :: debug
76
  LOGICAL :: Fini
76
  LOGICAL :: FiniS,FiniG,Fini
77 77
  LOGICAL, SAVE :: FRST=.TRUE.
78 78

  
79 79
  ! Variables
80 80
  INTEGER(KINT) :: IOpt, I,J,Iat, IBEG
81 81
  REAL(KREAL), ALLOCATABLE :: GradTmp(:),ZeroVec(:) ! 3*Nat or 3*Nat-6
82
  REAL(KREAL), ALLOCATABLE :: Step(:) ! 3*Nat or 3*Nat-6
82 83
  REAL(KREAL), ALLOCATABLE :: GeomOld(:),GradOld(:) ! (3*Nat or 3*Nat-6)
83 84
  REAL(KREAL), ALLOCATABLE :: GeomRef(:) ! (3*Nat or 3*Nat-6)
84 85
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! 3*Nat or 3*Nat-6
......
90 91
  REAL(KREAL), ALLOCATABLE :: NewGeom(:) !NCoord=3*Nat-6
91 92
  REAL(KREAL), ALLOCATABLE :: NewGradTmp(:)
92 93
  REAL(KREAL), ALLOCATABLE :: Hess_local_inv(:,:) ! used in diis
93
  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep
94
  REAL(KREAL) :: NormStep, FactStep, HP, MaxStep,NormGrad
94 95
  REAL(KREAL) :: Eold, Eini
95 96
  ! Values contains the values for the geometry analizes
96 97
  REAL(KREAL), ALLOCATABLE :: Values(:) ! NbVar
98
  CHARACTER(LCHARS) :: Line
97 99

  
98 100
  debug=valid('OptGeom')
99 101

  
......
104 106

  
105 107
  if (debug) Call Header("Entering Geom Optimization")
106 108

  
107
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord))
109
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord),GradTmp(NCoord),Step(NCoord))
108 110
  ALLOCATE(GradOld(NCoord),GeomOld(NCoord),ZeroVec(NCoord))
109 111
  ALLOCATE(GeomRef(NCoord))
110 112
  ALLOCATE(Hess_local(NCoord,NCoord))
......
144 146
     IF (Step_method .EQ. "RFO") Then
145 147
        Hess_local=0.
146 148
        DO I=1,NCoord
147
           Hess_local(I,I)=0.5d0
149
           Hess_local(I,I)=1.d0
148 150
        END DO
149 151
     END IF
150 152

  
......
211 213
     END IF
212 214

  
213 215
     !Print *, 'Hess_local after inversion='
214
     DO I=1,NCoord
216
!     DO I=1,NCoord
215 217
        !   WRITE(*,'(3(1X,F6.3))') Hess_local(I,1:3)
216
     END DO
218
!     END DO
217 219

  
218 220
     IF (Step_method .EQ. "RFO") Then
219 221
        Hess_local=0.
......
243 245
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
244 246
     IOpt=IOpt+1
245 247

  
248
     Write(Line,'(1X,A,I5)') "Iteration ",Iopt
249
     Call Header(TRIM(Line))
250
     WRITE(IoOut,*) "Current Geometry"
251
     Line=""
252
     Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
253
     if (COORD/="CART") THEN
254
        WRITE(IoOut,*) "Current Geometry in Cartesian coordinates"
255
        Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)
256
     END IF
257

  
258
     WRITE(IoOut,*) "Computing energy and gradient"
246 259
     ! Calculate the  energy and gradient
247 260
     IF (debug) THEN 
248 261
        WRITE(*,*) 'L134, DBG Opt_Geom, COORD=',TRIM(COORD),' IOpt=',IOpt
......
266 279
     ELSE
267 280
        Call EGrad(E,Geom,GradTmp,NCoord,IGeom,IOpt,GeomCart)
268 281
     END IF
282
!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!
283
! We have a pb with the format of GradTmp(Nat,3) whereas it is 3,Nat in Egrad
284
!
285
! This is a path for CART coordinates !!!
286
     IF (COORD=="CART") THEN
287
        Step=reshape(reshape(GradTmp,(/3,Nat/),ORDER=(/2,1/)),(/3*Nat/))
288
        GradTmp=Step
289
     END IF
290
!!!!!!!!!!!!!!!!!! PFL March 2013 !!!!!!!!!!!!!!!!!!!!
269 291

  
292
     WRITE(IoOut,'(1X," Energy E=",F15.6,A)') E,TRIM(UnitProg)
293
     WRITE(IoOut,*) "Gradient for current geometry"
294
     Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
295

  
270 296
     If (Iopt==1) EIni=E
271 297

  
272 298
     IF (debug) THEN 
......
329 355
     ! GradTmp becomes Step in Step_RFO_all.
330 356
     SELECT CASE (Step_method)
331 357
     CASE ('RFO')
332
        Call Step_RFO_all(NCoord,GradTmp,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
358
        Call Step_RFO_all(NCoord,Step,IGeom,Geom,GradTmp,Hess_local,ZeroVec)
359
        GradTmp=Step
333 360
     CASE ('GDIIS')
334 361
        Call Step_DIIS(NewGeom,Geom,NewGradTmp,GradTmp,HP,E,Hess_local,NCoord,FRST)
335 362
        ! q_{m+1} = q'_{m+1} - H^{-1}g'_{m+1}
......
367 394
     NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
368 395
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L215, NormStep=', NormStep, 'Step Threshold=', SThresh
369 396
     FactStep=min(1.d0,MaxStep/NormStep)
370
     Fini=(NormStep.LE.SThresh)
371
     NormStep=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
372
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormStep (Gradient)=', NormStep, 'Gradient Threshold=', GThresh
373
     Fini=(Fini.AND.(NormStep.LE.GThresh))
397
     FiniS=((NormStep*FactStep).LE.SThresh)
398
     NormGrad=sqrt(DOT_Product(GradOld,GradOld)) ! GradOld is the value of gradients.
399
     if (debug) WRITE(*,'(A,E10.4,1X,A,E10.4)') 'DBG OptGeom L219, NormGrad (Gradient)=', NormGrad, 'Gradient Threshold=', GThresh
400
     FiniG=(NormGrad.LE.GThresh)
401
     Fini=FiniS.AND.FiniG
374 402
     if (debug) WRITE(*,*) "DBG OptGeom FactStep=",FactStep
375
     GradTmp=GradTmp*FactStep ! GradTmp is step size here, from Step_RFO_all.
403
     GradTmp=GradTmp*FactStep ! GradTmp is step  here, from Step_RFO_all.
376 404

  
405
     WRITE(IoOut,*) " Converged ?"
406
     WRITE(IoOut,'(1X,A18,1X,A15,1X,A13,5X,A3)') " Criterion","Current Value","Threshold","Ok?"
407
     IF (FiniS) THEN
408
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"YES"
409
     ELSE
410
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Stepsize :",NormStep,SThresh,"NO"
411
     END IF
412
     IF (FiniG) THEN
413
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"YES"
414
     ELSE
415
        WRITE(IoOUT,'(1X,A20,5X,F8.3,6X,F8.3,6X,A3)') "Grad norm :",NormGrad,GThresh,"NO"
416
     END IF
377 417

  
378 418
     if (DynMaxStep.AND.(IOpt>1)) THEN
379 419
        If (E<EOld) Then
......
387 427

  
388 428
     Call Check_step(Coord,Nat,NCoord,GeomOld,GradTmp)
389 429

  
390
     Geom=GeomOld+GradTmp ! GradTmp is step size here, from Step_RFO_all.
430
! We add the step to the old geom
431
     Geom=GeomOld+GradTmp
391 432

  
392 433
     EOld=E
393 434

  
......
446 487

  
447 488
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
448 489

  
449
  DEALLOCATE(GradTmp2,GeomTmp2,GradTmp)
450
  DEALLOCATE(GradOld,GeomOld)
451
  DEALLOCATE(Hess_local)
452
  DEALLOCATE(GeomCart_old)
453
  DEALLOCATE(NewGeom,NewGradTmp)
454
  DEALLOCATE(Hess_local_inv)
455

  
456 490
  IF (Fini) THEN
457
     WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles, Energy =",E
491
     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization for geom ",IGeom," converged in ",Iopt," cycles"
458 492
  ELSE
459
     WRITE(*,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles, Last Energy = ",E
493
     WRITE(Line,'(1X,A,I5,A,I5,A,F15.8)') "Optimization  for geom ",IGeom," *NOT* converged in ",Iopt," cycles"
460 494
  END IF
495
  Call Header(Line)
496
  WRITE(IoOut,*) "Last Energy E=",E
461 497

  
462
  WRITE(*,*) "Initial Geometry:"
463
  DO I=1, Nat
464
     WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
465
  END DO
466
  WRITE(*,*) "Final Geometry:"
467
  DO I=1, Nat
468
     WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
498
!  WRITE(*,*) "Initial Geometry:"
499
!  DO I=1, Nat
500
!     WRITE(*,'(3(1X,F8.4))')  XyzGeomI(IGeom,:,I)
501
!  END DO
502
  WRITE(IoOut,*) "Final Geometry:"
503

  
504
     Line=""
505
     Call PrintGeom(Line,Nat,NCoord,Geom,Coord,IOOUT,Atome,Order,OrderInv,IndZmat)
506
     if (COORD/="CART") THEN
507
        WRITE(IoOut,*) "Last Geometry in Cartesian coordinates"
508
        Call PrintGeom(Line,Nat,NCoord,GeomCart,"CART",IOOUT,Atome,Order,OrderInv,IndZmat)
509
     END IF
510

  
511

  
512
!  DO I=1, Nat
513
!     WRITE(*,'(3(1X,F8.4))') GeomCart(I,:)
469 514
     !IF (I .GT. 1) Then
470 515
     !       WRITE(*,'(F6.4)') Sqrt(((GeomCart(I,1)-GeomCart(I-1,1))*(GeomCart(I,1)-GeomCart(I-1,1)))&
471 516
     !        + ((GeomCart(I,2)-GeomCart(I-1,2))*(GeomCart(I,2)-GeomCart(I-1,2)))&
472 517
     !        + ((GeomCart(I,3)-GeomCart(I-1,3))*(GeomCart(I,3)-GeomCart(I-1,3))))
473 518
     !END IF
474
  END DO
519
!  END DO
475 520

  
476 521
  IF (COORD .EQ. "BAKER") Then
477 522
     I=0 ! index for the NPrim (NPrim is the number of primitive coordinates).
......
493 538
  END IF
494 539

  
495 540
  DEALLOCATE(GeomCart)
541
  DEALLOCATE(GradTmp2,GeomTmp2,GradTmp,Step)
542
  DEALLOCATE(GradOld,GeomOld)
543
  DEALLOCATE(Hess_local)
544
  DEALLOCATE(GeomCart_old)
545
  DEALLOCATE(NewGeom,NewGradTmp)
546
  DEALLOCATE(Hess_local_inv)
496 547

  
497 548
  if (debug) Call Header("Geom Optimization Over")
498 549

  
src/Io_module.f90 (revision 9)
69 69
  CHARACTER(LCHARS) :: Siesta_Label, Siesta_CoordFile
70 70
! Number of species used in Siesta
71 71
  INTEGER(KINT) :: Siesta_NbSpecies
72
! Mass number for each species
73
  REAL(KREAL), ALLOCATABLE :: ListSpecies(:) ! NbSpecies
72
! Mass number for each species (atomic number)
73
  INTEGER(KINT), ALLOCATABLE :: Siesta_SpeciesMass(:) ! NbSpecies
74 74
! Name of each species
75 75
  CHARACTER(LCHARS), ALLOCATABLE :: Siesta_SpeciesName(:) ! NbSpecies
76 76
! Species for each atom

Formats disponibles : Unified diff