Statistics
| Revision:

root / src / Path.f90 @ 9

History | View | Annotate | Download (89.5 kB)

1 1 pfleura2
! This programs generates and optimizes a reaction path
2 1 pfleura2
! between at least two endpoints.
3 1 pfleura2
!
4 1 pfleura2
! It uses NAMELIST for input, see below !
5 1 pfleura2
!
6 1 pfleura2
! v3.0
7 1 pfleura2
! First version entirely in Fortran.
8 1 pfleura2
! It uses routines from Cart (path creation and parameterisation)
9 1 pfleura2
! and from IRC06, especially routines for point optimization.
10 1 pfleura2
!
11 1 pfleura2
! TO DO:
12 1 pfleura2
! 1) Possibility to read and/or calculate some hessian structures
13 1 pfleura2
! 2) Use of internal coordinate to estimate the hessian for stable
14 1 pfleura2
! species in a much better way...
15 1 pfleura2
!
16 1 pfleura2
!!!!!!!!!!!!!!!
17 1 pfleura2
! v3.1
18 1 pfleura2
! The Hessian update includes neighbour points.
19 1 pfleura2
!
20 1 pfleura2
! v3.2
21 1 pfleura2
! We add the option Hinv that uses the BFGS update of the Hessian inverse
22 1 pfleura2
!
23 1 pfleura2
! v3.3
24 1 pfleura2
! We add the full option for ZMAT coordinates, where all is done using
25 1 pfleura2
! internal coordinates.
26 1 pfleura2
! v3.31
27 1 pfleura2
! The step size is controlled using the step norm and not the maximum step
28 1 pfleura2
! component.
29 1 pfleura2
! v3.5
30 1 pfleura2
! Big modifications of Calc_zmat_const_frag in order to be able
31 1 pfleura2
! to treat VASP geometries.
32 1 pfleura2
! For now, only XYZ geometries are read and written. VASP interface
33 1 pfleura2
! has to be written.
34 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
35 1 pfleura2
! v3.6
36 1 pfleura2
! Some minor modif in Calc_zmat_const_frag to have nicer programming
37 1 pfleura2
! and a new mode added: Coord=MIXED where part of the system
38 1 pfleura2
! is extrapolated in cartesian coordinates and the rest of it in zmat.
39 1 pfleura2
! it might be Useful for VASP, or when lots of atoms are frozen
40 1 pfleura2
! by default, when Coord=MIXED, all frozen atoms are described in cartesian
41 1 pfleura2
! coordinates.
42 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43 1 pfleura2
! v3.61
44 1 pfleura2
! We move further and further to VASP program. New variables added to the 'path' namelist:
45 1 pfleura2
! input: what is the program used for the geometries ?
46 1 pfleura2
!     Possible values: xyz (or cart) for Xmol format; VASP. Xyz should be used for Gaussian
47 1 pfleura2
! prog: what is the program used for the calculations ?
48 1 pfleura2
!     Possible values for now: gaussian, vasp.
49 1 pfleura2
! For now (02/2007), when prog='vasp' is used, only the initial path is created.
50 1 pfleura2
! That means that we have only added input/output subroutines :-)
51 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
52 1 pfleura2
! v3.7
53 1 pfleura2
! We try to improve the optimization steps. For that, we first modify many routines
54 1 pfleura2
! so that the new step is calculated in the 'free' degrees of freedom, as it is done
55 1 pfleura2
! in IRC07/ADF for example. That will allow us to impose constraints in an easier way.
56 1 pfleura2
!
57 1 pfleura2
! v3.701
58 1 pfleura2
! Add a new option for the initial extrapolation of the Hessian
59 1 pfleura2
!
60 1 pfleura2
!v3.71
61 1 pfleura2
! The optimization step has been improved... but the vfree part has not yet been included.
62 1 pfleura2
! This has still to be done :-)
63 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
64 1 pfleura2
! v3.8
65 1 pfleura2
! I have included part of the vfree procedure: for now freemv returns the Id
66 1 pfleura2
! matrix, that is used to project out the tangent vector.
67 1 pfleura2
! This improves the optimization a lot, but I will have to implement
68 1 pfleura2
! a real freemw subroutine to have constraints there !
69 1 pfleura2
! v3.81 Technical version
70 1 pfleura2
! I have added a new program 'Ext'. When this one is used, Path creates a Ext.in file,
71 1 pfleura2
! calls Ext.exe file to get a Ext.out file that contains the Energy and the gradient.
72 1 pfleura2
! Format:
73 1 pfleura2
! Ext.in : Xmol format. The second line (comment in Xmol) contains the COORD name.
74 1 pfleura2
! Ext.out: Energy on firt line (1X,D25.15), then gradient in CARTESIAN coordinates
75 1 pfleura2
!  (3(1X,D25.15)) format. Natoms lines.
76 1 pfleura2
! TO DO:  Gradient could be given in COORD format, ie cartesian for COORD=CART, XYZ or HYBRID
77 1 pfleura2
! ZMAT for CORD=ZMAT (in that case, 6 zeros are there at the begining !)
78 1 pfleura2
! For MIXED, it is equivalent to CART :-)
79 1 pfleura2
! v3.811
80 1 pfleura2
! I have added a prog="TEST" option, where egradient calls egrad_test subroutine.
81 1 pfleura2
! It then calculates the energy and gradient as it wants and it returns the cartesian
82 1 pfleura2
! gradient. The purpose is to have an analytical approximation of the PES of the system
83 1 pfleura2
! under study to gain some time in testing the program. This is not documented :-)
84 1 pfleura2
!
85 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
86 1 pfleura2
!  v3.9
87 1 pfleura2
! In order to be fully interfaced with VASP, I have to change the architecture
88 1 pfleura2
! of the program. In particular, with VASP, one can use a Para8+NEB routine
89 1 pfleura2
! of VASP to calculate energy+forces for all points of the path.
90 1 pfleura2
! v3.91
91 1 pfleura2
! One small modification: in the cart list, one can freeze a whole sequence by
92 1 pfleura2
! entering a negative number for the last atom of the sequence
93 1 pfleura2
! for example : cart=1 -160 165 -170 will define atoms from 1 to 160 and from
94 1 pfleura2
! 165 to 170 (included) as cartesian atoms.
95 1 pfleura2
! Of course, Same applies to Frozen
96 1 pfleura2
! v3.92 / v3.93
97 1 pfleura2
! Slight modifications. One noticeable: when using Vasp, the program analyses the
98 1 pfleura2
! displacements to suggest which atoms could be described in CART.
99 1 pfleura2
! v3.94
100 1 pfleura2
! VaspMode=Para now implemented. One needs to be careful in defining 'ProgExe'.
101 1 pfleura2
! In particular, one must put the mpirun command in ProgExe !!! There is NO test.
102 1 pfleura2
! v3.95
103 1 pfleura2
! When using internal coordinates (zmat, hybrid or mixes), Path calculates the number
104 1 pfleura2
! of fragments for each geometry. The reference one (ie the one used to compute the internal
105 1 pfleura2
! coordinates) is now the one with the least fragments.
106 5 pfleura2
! user can also choose one geometry by specifying IGeomRef in the Path namelist.
107 1 pfleura2
! v3.96
108 1 pfleura2
! I have added a test for the geometry step: it does not allow valence angle to
109 1 pfleura2
! be negative or bigger than 180. I hope that this might help geometry optimization
110 1 pfleura2
! when dealing with nearly linear molecules.
111 1 pfleura2
! v3.961
112 1 pfleura2
! The options cart and frozen are now LOGICAL:
113 1 pfleura2
! Fcart=T indicates that a &cartlist namelist should be read !
114 1 pfleura2
! Ffrozen=T indicates that a &frozenlist namelist should be read !
115 1 pfleura2
! TO DO: Autocart=T indicates that the atoms described in cartesian coordiantes should
116 1 pfleura2
! be found automatically
117 1 pfleura2
! v3.97
118 1 pfleura2
! Baker has been implemented by P. Dayal: we are interpolating Xint for now.
119 1 pfleura2
! TO do: interpolate Umat ?
120 1 pfleura2
! -----
121 1 pfleura2
! v3.971
122 1 pfleura2
! AutoCart=T is being implemented. Goal: having a 'black box' version for Vasp users.
123 1 pfleura2
! Vmd: True if user wants to watch the path using VMD. Only used for VASP for now.
124 1 pfleura2
!------------
125 1 pfleura2
! v3.98
126 1 pfleura2
! Calc_zmat has been replaced by Calc_zmat_frag which is more robust.
127 1 pfleura2
! Still missing: linear molecules and some tests.
128 1 pfleura2
!
129 1 pfleura2
! v3.981
130 1 pfleura2
! * Added a non document flag in .valid file that switch on the numerical
131 1 pfleura2
! calculation of the hessian instead of the update.
132 1 pfleura2
! * Added HesUpd variable in the namelist. Default value is MS for path
133 1 pfleura2
!   optimization because it does not imposes a positive hessian, and BFGS
134 1 pfleura2
!  for geometry optimization
135 1 pfleura2
! v 4.0
136 1 pfleura2
! * This version has been made public within the CARTE project.
137 1 pfleura2
! * Added:
138 1 pfleura2
! - Step_method to choose the optimization method: RFO or Diis
139 1 pfleura2
! - DynMaxStep: if TRUE, the maximum size of a step is updated at each step.
140 1 pfleura2
!             If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2.
141 1 pfleura2
!             It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]
142 1 pfleura2
!
143 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!
144 1 pfleura2
! v4.1
145 1 pfleura2
! * Para mode has been partly implemented for Gaussian.
146 1 pfleura2
!  VaspMode is thus now 'RunMode' and can be Serial or Para
147 1 pfleura2
! Only Gaussian and VASP are supported for Para mode now.
148 1 pfleura2
! v4.12
149 1 pfleura2
! Some bugs in unconstrained zmat construction are corrected.
150 1 pfleura2
! v4.13
151 1 pfleura2
! Prakash has implemented the GEDIIS optimization procedure.
152 1 pfleura2
! v4.14
153 1 pfleura2
! There is a small problem for some inputs in VASP, where the
154 1 pfleura2
! last geometry is not the one given by the user.
155 1 pfleura2
! It seems to come from the fact that some people try to freeze
156 1 pfleura2
! some unfrozen atoms
157 1 pfleura2
! So some tests have been added to check that frozen atoms do not move.
158 1 pfleura2
!!!
159 1 pfleura2
! v4.15
160 1 pfleura2
! There were some bugs when reading and expanding cartlist and frozenlist
161 1 pfleura2
! I have changed the way frozen atoms are compared by adding partial alignment.
162 1 pfleura2
! v4.16
163 1 pfleura2
! Some bugs were corrected.
164 1 pfleura2
! v4.17
165 1 pfleura2
! Added MOPAC as a program.
166 1 pfleura2
! v4.175
167 1 pfleura2
! Added Three more parameters for defining the input and output names for
168 1 pfleura2
! the calculations.
169 1 pfleura2
! CalcName: Prefix for the files used for the energy and gradient calculations
170 1 pfleura2
! ISuffix: Suffix for the input file
171 1 pfleura2
! OSuffix: suffix for the output file.
172 1 pfleura2
! This maters for Gaussian for example.
173 1 pfleura2
!
174 1 pfleura2
! v4.177 (P) 2010 Nov 22
175 1 pfleura2
! - We add TurboMole as a new external code.
176 1 pfleura2
! - Subtle change for Input: the default is XYZ whatever the code, except
177 1 pfleura2
! for VASP.
178 1 pfleura2
!
179 1 pfleura2
! v4.178 (P) 2010 Nov 28
180 1 pfleura2
! - We add the possibility to call CHAMFRE reactive force field using
181 1 pfleura2
! prog=TEST2 or CHAMFRE
182 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183 1 pfleura2
! OpenPath v1.4  2011 May-June
184 1 pfleura2
! - Also added two flags: CalcEReac and CalcEProd. If .TRUE. Path will
185 1 pfleura2
! compute the energy of the reactant
186 1 pfleura2
! (and/or Product) at the begining of the calculations.
187 1 pfleura2
! This is useful if you start from non optimized points.
188 1 pfleura2
!
189 1 pfleura2
! - Also added more methods for Hessian updates. The following methods are
190 1 pfleura2
! now availables:
191 1 pfleura2
!   For inverse Hessian: BFGS, DFP, MS.
192 1 pfleura2
!   For Hessian: BFGS, DFP, MS, Bofill
193 1 pfleura2
!
194 1 pfleura2
! -Small change: HesUpd is replaced by Hupdate
195 1 pfleura2
!
196 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197 1 pfleura2
! TO DO
198 1 pfleura2
! - Trying to improve the output of Path. In particular, we now create
199 1 pfleura2
! a .datl file containing the image index, curvilinear coordinate and energy
200 1 pfleura2
! in the spirit of the AnaPath tool.
201 1 pfleura2
! AnaPath is still usefull as it can extract geometrical parameters on top
202 1 pfleura2
! of the energy
203 1 pfleura2
! - When using Hessian update, it is more efficient to use Cholesky
204 1 pfleura2
! decomposition, and to update this decomposition rather than the Hessian.
205 1 pfleura2
! See Nocedal&Wright p141.
206 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
207 5 pfleura2
!
208 1 pfleura2
! Project Name has changed to Opt'n Path
209 1 pfleura2
! OpenPaht 2011 May-June -> Opt'n Path v1.4
210 5 pfleura2
!
211 5 pfleura2
! 2013 Jan 18 Opt'n Path v1.47 (going to 1.5 slowly)
212 5 pfleura2
! We add Siesta (and start thinking of add CP2K)
213 5 pfleura2
!  as a new "energy engine" code
214 5 pfleura2
! TO DO to go to 1.5: finish cleaning...
215 7 pfleura2
!
216 5 pfleura2
! 2013 Feb Opt'n Path v1.48
217 5 pfleura2
! We add some keyword for more output for Geometry Optimizations
218 5 pfleura2
! GeomFile: name of the file to print out the geometries and their energy
219 5 pfleura2
! as Xyz file. (only format for now)
220 7 pfleura2
!
221 7 pfleura2
! 2013 Feb opt'n Path v1.49
222 7 pfleura2
! We add the possibility to analyze the geometries in terms
223 7 pfleura2
! of Bonds, Angles, dihedrals
224 7 pfleura2
! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path)
225 7 pfleura2
! This is basicaly the merging of AnaPath into Opt'n Path... at last !
226 7 pfleura2
! This is done by adding a Flag into the Path namelist:
227 7 pfleura2
! AnaGeom: logical, if true Path looks for the AnaList namelist
228 7 pfleura2
! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path
229 7 pfleura2
! will read it and print the values of the variable in a Path.dat file
230 7 pfleura2
! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy.
231 1 pfleura2
232 1 pfleura2
PROGRAM PathOpt
233 1 pfleura2
234 1 pfleura2
  use VarTypes
235 1 pfleura2
  use Path_module
236 1 pfleura2
  use Io_module
237 1 pfleura2
238 1 pfleura2
  IMPLICIT NONE
239 1 pfleura2
240 5 pfleura2
  CHARACTER(132) :: FileIn, FileOut, Line
241 1 pfleura2
242 1 pfleura2
  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF
243 1 pfleura2
  INTEGER(KINT) :: IOpt
244 5 pfleura2
  INTEGER(KINT) :: Idx
245 1 pfleura2
  INTEGER(KINT) :: N3at, NFree, Nfr
246 1 pfleura2
247 1 pfleura2
  INTEGER(KINT) :: List(2*MaxFroz)
248 1 pfleura2
249 5 pfleura2
  REAL(KREAL) :: E,FactStep !E=Energy
250 1 pfleura2
  REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad
251 1 pfleura2
  REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E
252 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord
253 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord)
254 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord
255 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: HessTmp(:,:)
256 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
257 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
258 1 pfleura2
  LOGICAL, SAVE :: allocation_flag
259 1 pfleura2
260 1 pfleura2
  ! Flag to see if frozen atoms are frozen (see below)
261 1 pfleura2
  LOGICAL :: FChkFrozen
262 1 pfleura2
263 1 pfleura2
  ! Energies for all points along the path at the previous iteration
264 7 pfleura2
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF
265 1 pfleura2
  ! Maximum step allowed for this geometry
266 7 pfleura2
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF
267 1 pfleura2
268 1 pfleura2
! these are used to read temporary coordinates
269 1 pfleura2
  LOGICAL :: FFrozen,FCart
270 1 pfleura2
271 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272 1 pfleura2
  !
273 1 pfleura2
  ! To analyse the frozen atoms
274 1 pfleura2
  !
275 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat
276 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat
277 1 pfleura2
  REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist
278 1 pfleura2
  LOGICAL, ALLOCATABLE :: ListAtFroz(:)
279 1 pfleura2
  INTEGER(KINT) :: Iat1, Iat2, Iat3
280 1 pfleura2
281 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282 1 pfleura2
283 1 pfleura2
  LOGICAL :: Debug, Fini,Flag_Opt_Geom
284 1 pfleura2
285 1 pfleura2
!  INTEGER(KINT), EXTERNAL :: Iargc
286 1 pfleura2
  INTEGER(KINT), EXTERNAL ::  ConvertNumAt
287 7 pfleura2
  INTEGER(KINT) :: IoS
288 1 pfleura2
289 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
290 1 pfleura2
  !
291 1 pfleura2
  ! For Debugging purposes
292 1 pfleura2
  !
293 1 pfleura2
  LOGICAL :: FCalcHess
294 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
295 1 pfleura2
296 1 pfleura2
297 1 pfleura2
  INTERFACE
298 1 pfleura2
     function valid(string) result (isValid)
299 1 pfleura2
       CHARACTER(*), intent(in) :: string
300 1 pfleura2
       logical                  :: isValid
301 1 pfleura2
     END function VALID
302 1 pfleura2
303 1 pfleura2
     SUBROUTINE Read_Geom(Input)
304 1 pfleura2
       CHARACTER(32), INTENT(IN) :: input
305 1 pfleura2
     END SUBROUTINE READ_GEOM
306 1 pfleura2
307 5 pfleura2
308 1 pfleura2
     subroutine hupdate_all (n, dx, dgrad, HessOld)
309 1 pfleura2
310 1 pfleura2
       IMPLICIT NONE
311 1 pfleura2
312 1 pfleura2
       INTEGER, PARAMETER :: KINT=KIND(1)
313 1 pfleura2
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
314 1 pfleura2
315 1 pfleura2
       INTEGER(KINT), INTENT(IN) :: n
316 1 pfleura2
       real(KREAL) :: dx(n), dgrad(n), HessOld(n*n)
317 1 pfleura2
     END subroutine hupdate_all
318 1 pfleura2
319 1 pfleura2
     SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess)
320 1 pfleura2
321 1 pfleura2
       IMPLICIT NONE
322 1 pfleura2
323 1 pfleura2
       INTEGER, PARAMETER :: KINT=KIND(1)
324 1 pfleura2
       INTEGER, PARAMETER :: KREAL=KIND(1.D0)
325 1 pfleura2
326 1 pfleura2
       INTEGER(KINT), INTENT(IN) :: NFREE
327 1 pfleura2
       REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree)
328 1 pfleura2
       REAL(KREAL), INTENT(IN) :: ds(Nfree)
329 1 pfleura2
       REAL(KREAL), INTENT(IN) :: dGrad(Nfree)
330 1 pfleura2
331 1 pfleura2
     END subroutine Hinvup_BFGS_new
332 1 pfleura2
333 1 pfleura2
334 8 pfleura2
    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)
335 8 pfleura2
336 8 pfleura2
      Use VarTypes
337 8 pfleura2
338 8 pfleura2
      implicit none
339 8 pfleura2
! Input
340 8 pfleura2
      CHARACTER(*), INTENT(IN) :: Line
341 8 pfleura2
      CHARACTER(*), INTENT(IN) :: String
342 8 pfleura2
      LOGICAL, OPTIONAL, INTENT(IN) :: Case
343 8 pfleura2
      LOGICAL, OPTIONAL, INTENT(IN) :: Back
344 8 pfleura2
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
345 8 pfleura2
346 8 pfleura2
! Output
347 8 pfleura2
! the position of String in Line (the first one) unless Back is present
348 8 pfleura2
      INTEGER(KINT) :: Pos
349 8 pfleura2
    END FUNCTION InString
350 8 pfleura2
351 8 pfleura2
    SUBROUTINE die(routine, msg, file, line, unit)
352 8 pfleura2
353 8 pfleura2
      Use VarTypes
354 8 pfleura2
      Use io_module
355 8 pfleura2
356 8 pfleura2
      implicit none
357 8 pfleura2
      character(len=*), intent(in)           :: routine, msg
358 8 pfleura2
      character(len=*), intent(in), optional :: file
359 8 pfleura2
      integer(KINT), intent(in), optional      :: line, unit
360 8 pfleura2
361 8 pfleura2
    END SUBROUTINE die
362 8 pfleura2
363 8 pfleura2
    SUBROUTINE Warning(routine, msg, file, line, unit)
364 8 pfleura2
365 8 pfleura2
      Use VarTypes
366 8 pfleura2
      Use io_module
367 8 pfleura2
368 8 pfleura2
      implicit none
369 8 pfleura2
370 8 pfleura2
      character(len=*), intent(in)           :: routine, msg
371 8 pfleura2
      character(len=*), intent(in), optional :: file
372 8 pfleura2
      integer(KINT), intent(in), optional      :: line, unit
373 8 pfleura2
374 8 pfleura2
    END SUBROUTINE Warning
375 8 pfleura2
376 8 pfleura2
377 1 pfleura2
  END INTERFACE
378 1 pfleura2
379 1 pfleura2
  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
380 1 pfleura2
       MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, &
381 1 pfleura2
       SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, &
382 1 pfleura2
       input,prog, ProgExe,RunMode, AtTypes,poscar,PathName,  &
383 1 pfleura2
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
384 1 pfleura2
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
385 1 pfleura2
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
386 8 pfleura2
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile
387 1 pfleura2
388 1 pfleura2
  Namelist/cartlist/list
389 1 pfleura2
  Namelist/frozenlist/list
390 7 pfleura2
  Namelist/analist/nb
391 1 pfleura2
392 1 pfleura2
393 1 pfleura2
  Flag_Opt_Geom = .FALSE. ! Initialized.
394 1 pfleura2
395 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!
396 1 pfleura2
  !
397 1 pfleura2
  ! We print the Version info in file
398 1 pfleura2
  !
399 5 pfleura2
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
400 1 pfleura2
401 8 pfleura2
  ! We read the files  names
402 1 pfleura2
!  SELECT CASE (iargc())
403 1 pfleura2
  SELECT CASE (command_argument_count())
404 1 pfleura2
  CASE (2)
405 1 pfleura2
     call getarg(1,FileIn)
406 1 pfleura2
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
407 1 pfleura2
     call getarg(2,FileOut)
408 1 pfleura2
     OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN')
409 1 pfleura2
  CASE (1)
410 1 pfleura2
     call getarg(1,FileIn)
411 1 pfleura2
     IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
412 9 pfleura2
        WRITE(*,*) "Opt'n Path mini-help"
413 1 pfleura2
        WRITE(*,*) "--------------"
414 1 pfleura2
        WRITE(*,*) ""
415 1 pfleura2
        WRITE(*,*) "Use: Path Input_file Output_file"
416 1 pfleura2
        WRITE(*,*) "Input_file starts with a Namelist called path"
417 1 pfleura2
        WRITE(*,*) ""
418 9 pfleura2
! The following has been generated (March 19, 2013) using the Mini_help.tex file
419 9 pfleura2
! 1) Edit Mini_help.tex
420 9 pfleura2
! 2) latex Mini_help.tex
421 9 pfleura2
! 3) catdvi -e 1 -U  Mini_help.dvi | sed -re "s/\[U\+2022\]/-/g" | sed -re "s/([^^[:space:]])\s+/\1 /g" > file.txt
422 9 pfleura2
! 4) edit file.txt to keep only the following lines, and to reformat a bit
423 9 pfleura2
! 5) awk '{print "WRITE(*,*) \"" $0 "\""}' file.txt > help
424 9 pfleura2
! 5) cut&paste help in Path.f90 ...
425 9 pfleura2
   WRITE(*,*) " &path"
426 9 pfleura2
   WRITE(*,*) " nat=3, ! Number of atoms"
427 9 pfleura2
   WRITE(*,*) " ngeomi=3, ! Number of initial geometries"
428 9 pfleura2
   WRITE(*,*) " ngeomf=12, !Number of geometries along the path"
429 9 pfleura2
   WRITE(*,*) " OptReac=.T., ! Do you want to optimize the reactants ?"
430 9 pfleura2
   WRITE(*,*) " OptProd=.T., ! Optimize the products"
431 9 pfleura2
   WRITE(*,*) " coord='zmat', ! We use Z-matrix coordinates"
432 9 pfleura2
   WRITE(*,*) " maxcyc=31, ! Max number of iterations"
433 9 pfleura2
   WRITE(*,*) " IReparam=2,! re-distribution of points along the path every 2 iterations"
434 9 pfleura2
   WRITE(*,*) " ISpline=50, ! Start using spline interpolation at iteration 50"
435 9 pfleura2
   WRITE(*,*) " Hinv=.T. , ! Use inverse of the Hessian internally (default: T)"
436 9 pfleura2
   WRITE(*,*) " MW=T, ! Works in Mass Weighted coordiante (default T)"
437 9 pfleura2
   WRITE(*,*) " PathName='Path_HCN_zmat_test', ! Name of the file used for path outputs"
438 9 pfleura2
   WRITE(*,*) " prog='gaussian',! we use G03 to get energy and gradients"
439 9 pfleura2
   WRITE(*,*) " SMax=0.1 ! Displacement cannot exceed 0.1 atomic units (or mass weighted at. unit)"
440 9 pfleura2
   WRITE(*,*) " /"
441 9 pfleura2
   WRITE(*,*) "  3"
442 9 pfleura2
   WRITE(*,*) " Energy : 0.04937364"
443 9 pfleura2
   WRITE(*,*) " H 0.0000 0.0000 0.0340"
444 9 pfleura2
   WRITE(*,*) " C 0.0000 0.0000 1.1030"
445 9 pfleura2
   WRITE(*,*) " N 0.0000 0.0000 2.2631"
446 9 pfleura2
   WRITE(*,*) "  3"
447 9 pfleura2
   WRITE(*,*) " Energy : 0.04937364"
448 9 pfleura2
   WRITE(*,*) " H 0.0000 1.1000 1.1030"
449 9 pfleura2
   WRITE(*,*) " C 0.0000 0.0000 1.1030"
450 9 pfleura2
   WRITE(*,*) " N 0.0000 0.0000 2.2631"
451 9 pfleura2
   WRITE(*,*) "3"
452 9 pfleura2
   WRITE(*,*) " CNH"
453 9 pfleura2
   WRITE(*,*) "H 0.000000 0.000000 3.3"
454 9 pfleura2
   WRITE(*,*) "C 0.000000 0.000000 1.1"
455 9 pfleura2
   WRITE(*,*) "N 0.000000 0.000000 2.26"
456 9 pfleura2
   WRITE(*,*) "%chk=Test"
457 9 pfleura2
   WRITE(*,*) "#P AM1 FORCE"
458 9 pfleura2
   WRITE(*,*) " HCN est bien"
459 9 pfleura2
   WRITE(*,*) ""
460 9 pfleura2
   WRITE(*,*) "0,1"
461 9 pfleura2
   WRITE(*,*) "H 0.000000 0.000000 0.000000"
462 9 pfleura2
   WRITE(*,*) "C 0.000000 0.000000 1.000"
463 9 pfleura2
   WRITE(*,*) "N 0.000000 0.000000 3.00"
464 9 pfleura2
   WRITE(*,*) ""
465 9 pfleura2
   WRITE(*,*) "*******"
466 9 pfleura2
   WRITE(*,*) "* Compulsory variables are:"
467 9 pfleura2
   WRITE(*,*) "*******"
468 9 pfleura2
   WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
469 9 pfleura2
   WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
470 9 pfleura2
   WRITE(*,*) "Nat: Number of atoms"
471 9 pfleura2
   WRITE(*,*) ""
472 9 pfleura2
   WRITE(*,*) "*******"
473 9 pfleura2
   WRITE(*,*) "* Other options:"
474 9 pfleura2
   WRITE(*,*) "*******"
475 9 pfleura2
   WRITE(*,*) "Input: String that indicates the type of the input geometries. Accepted values are: Cart (or"
476 9 pfleura2
   WRITE(*,*) "    Xmol or Xyz), Vasp, Turbomole or Siesta."
477 9 pfleura2
   WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
478 9 pfleura2
   WRITE(*,*) "    Accepted values are: Gaussian, Mopac, Vasp, Turbomole, Siesta or Ext."
479 9 pfleura2
   WRITE(*,*) ""
480 9 pfleura2
   WRITE(*,*) ""
481 9 pfleura2
   WRITE(*,*) "        - In case of a Gaussian calculations, input must be set to Cart. One example of a gaus-"
482 9 pfleura2
   WRITE(*,*) "          sian input should be added at the end of the input file.See example file Test_HCN_zmat_g03.path."
483 9 pfleura2
   WRITE(*,*) "        - In the case of a VASP calculation, if input is set to Cart, then the preamble of a"
484 9 pfleura2
   WRITE(*,*) "          VASP calculation s"
485 9 pfleura2
   WRITE(*,*) "        - Mopac: Examples using sequential call to MOPAC2009 to calculate the energies and"
486 9 pfleura2
   WRITE(*,*) "          forces along the path. hould be added at the end of the input file. See example file"
487 9 pfleura2
   WRITE(*,*) "          Test_VASP_cart.path. In the case of a VASP calculation, one should also give value"
488 9 pfleura2
   WRITE(*,*) "          of the RunMode variable ."
489 9 pfleura2
   WRITE(*,*) "        - In the case of a SIESTA calculation, an example of a Siesta input file should be added"
490 9 pfleura2
   WRITE(*,*) "          at the end of the input file. See Test_Siesta.path."
491 9 pfleura2
   WRITE(*,*) ""
492 9 pfleura2
   WRITE(*,*) "RunMode: When running on a multi-processor machine, this indicates wether Opt'n Path"
493 9 pfleura2
   WRITE(*,*) "    should calculate the energy and gradient of the whole path in parallel or not. User has"
494 9 pfleura2
   WRITE(*,*) "    two options. One is to calculate the energy and gradient of each point sequentially. This"
495 9 pfleura2
   WRITE(*,*) "    is usefull when running on one (or two) processors. In this case, RunMode should be put"
496 9 pfleura2
   WRITE(*,*) "    to SERIAL.When running in parallel with 8 or more processors, one can use VASP to"
497 9 pfleura2
   WRITE(*,*) "    calculate simultaneously the energies and gradients for all points, as in a normal NEB cal-"
498 9 pfleura2
   WRITE(*,*) "    culation. In this case, RunMode must be set to PARA. For now, this is usefull only for VASP."
499 9 pfleura2
   WRITE(*,*) ""
500 9 pfleura2
   WRITE(*,*) "ProgExe: Name (with full path) of the executable to be used to get energies and gradients."
501 9 pfleura2
   WRITE(*,*) "    For example, if VASP is used in parallel, one might have something like:"
502 9 pfleura2
   WRITE(*,*) "    ProgExe='/usr/local/mpich/bin/mpirun -machinefile machine -np 8 ~/bin/VASP_46'."
503 9 pfleura2
   WRITE(*,*) "    Another option that I use, is to put ProgExe='./run_vasp' and I put every option needed"
504 9 pfleura2
   WRITE(*,*) "    to run VASP into the run_vasp file."
505 9 pfleura2
   WRITE(*,*) "EReac: (REAL) By default, Opt'n Path does not compute the energy of the reactants and"
506 9 pfleura2
   WRITE(*,*) "    products. This thus indicates the reactants energy to Opt'n Path to have better plots at"
507 9 pfleura2
   WRITE(*,*) "    the end."
508 9 pfleura2
   WRITE(*,*) "EProd: (REAL) By default, Opt'n Path does not compute the energy of the reactants and"
509 9 pfleura2
   WRITE(*,*) "    products. This thus indicates the products energy to Opt'n Path to have better plots."
510 9 pfleura2
   WRITE(*,*) ""
511 9 pfleura2
   WRITE(*,*) "PathName: Prefix used to save the path. Default is Path."
512 9 pfleura2
   WRITE(*,*) "Poscar: string that will be used as the prefix for the different POSCAR files in a VASP calcu-"
513 9 pfleura2
   WRITE(*,*) "    lations. Usefull only if PathOnly=.TRUE., not used for internal calculations."
514 9 pfleura2
   WRITE(*,*) "CalcName: Prefix for the files used for the energy and gradient calculations"
515 9 pfleura2
   WRITE(*,*) ""
516 9 pfleura2
   WRITE(*,*) "ISuffix: Suffix for the input file used for energy and gradient calculations. The full inputfile"
517 9 pfleura2
   WRITE(*,*) "    name will thus be CalcName.ISuffix."
518 9 pfleura2
   WRITE(*,*) "OSuffix: Suffix for the output file used for energy and gradient calculations. The full outputfile"
519 9 pfleura2
   WRITE(*,*) "    name will thus be CalcName.OSuffix."
520 9 pfleura2
   WRITE(*,*) ""
521 9 pfleura2
   WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates. Valid only for"
522 9 pfleura2
   WRITE(*,*) "     Coord=Zmat, Hybrid or Mixed."
523 9 pfleura2
   WRITE(*,*) "Fact: REAL used to define if two atoms are linked. If d(A,B) =< fact* (rcov(A) + rcov(B)),"
524 9 pfleura2
   WRITE(*,*) "     then A and B are considered Linked."
525 9 pfleura2
   WRITE(*,*) ""
526 9 pfleura2
   WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
527 9 pfleura2
   WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
528 9 pfleura2
   WRITE(*,*) "        - CART (or Xyz): works in cartesian"
529 9 pfleura2
   WRITE(*,*) "        - Zmat: works in internal coordinates (Zmat)"
530 9 pfleura2
   WRITE(*,*) "        - Mixed: frozen atoms, as well as atoms defined by the 'cart' array(see below) are"
531 9 pfleura2
   WRITE(*,*) "          describe in CARTESIAN, whereas the others are described in Zmat"
532 9 pfleura2
   WRITE(*,*) "        - Baker: use of Baker coordinates, also called delocalized internal coordinates"
533 9 pfleura2
   WRITE(*,*) "        - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
534 9 pfleura2
   WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:"
535 9 pfleura2
   WRITE(*,*) ""
536 9 pfleura2
   WRITE(*,*) "        - RFO: Rational function optimization"
537 9 pfleura2
   WRITE(*,*) "        - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
538 9 pfleura2
   WRITE(*,*) "HUpdate: method to update the hessian. By default, it is Murtagh-Sargent Except for geom-"
539 9 pfleura2
   WRITE(*,*) "     etry optimization where it is BFGS."
540 9 pfleura2
   WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
541 9 pfleura2
   WRITE(*,*) ""
542 9 pfleura2
   WRITE(*,*) "Smax: Maximum length of a step during path optimization"
543 9 pfleura2
   WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
544 9 pfleura2
   WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is"
545 9 pfleura2
   WRITE(*,*) "     taken"
546 9 pfleura2
   WRITE(*,*) ""
547 9 pfleura2
   WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
548 9 pfleura2
   WRITE(*,*) "     that is kept. FTan=1. corresponds to the full displacement, whereas FTan=0. gives a"
549 9 pfleura2
   WRITE(*,*) "     displacement orthogonal to the path."
550 9 pfleura2
   WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of"
551 9 pfleura2
   WRITE(*,*) "     reparameterization."
552 9 pfleura2
   WRITE(*,*) "IReparamT: When the path is not reparameterised at each iteration, this gives the frequency"
553 9 pfleura2
   WRITE(*,*) "     of reparameterization of the tangents."
554 9 pfleura2
   WRITE(*,*) ""
555 9 pfleura2
   WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path. This option indicates"
556 9 pfleura2
   WRITE(*,*) "     the first step where spline interpolation is used."
557 9 pfleura2
   WRITE(*,*) "BoxTol: Real between 0. and 1. When doing periodic calculations, it might happen that an"
558 9 pfleura2
   WRITE(*,*) "     atom moves out of the unit cell. Path detects this by comparing the displacement to"
559 9 pfleura2
   WRITE(*,*) "     boxtol: if an atom moves by more than Boxtol, then it is moved back into the unit cell."
560 9 pfleura2
   WRITE(*,*) "     Default value: 0.5."
561 9 pfleura2
   WRITE(*,*) "FrozTol: (Real) This indicates the threshold to determine wether an atom moves between two"
562 9 pfleura2
   WRITE(*,*) "     images. Default is 1e-4."
563 9 pfleura2
   WRITE(*,*) ""
564 9 pfleura2
   WRITE(*,*) "OptGeom: This INTEGER indicates the index of the geometry you want to optimize. If"
565 9 pfleura2
   WRITE(*,*) "    OptGeom is set, then Opt'n Path performs a geometry optimization instead of a path"
566 9 pfleura2
   WRITE(*,*) "    interpolation."
567 9 pfleura2
   WRITE(*,*) "GeomFile: Name of the file to print the geometries and their energy during a geometry opti-"
568 9 pfleura2
   WRITE(*,*) "    mization. If this variable is not given then nothing is printed."
569 9 pfleura2
   WRITE(*,*) "AnaFile: Name of the file to print the values of the geometrical parameters that are monitored"
570 9 pfleura2
   WRITE(*,*) "    if AnaGeom=.TRUE.. Default is PathName.dat"
571 9 pfleura2
   WRITE(*,*) ""
572 9 pfleura2
   WRITE(*,*) "GplotFile: Name of the gnuplot file to plot the evolution of the geometrical parameters that"
573 9 pfleura2
   WRITE(*,*) "    are monitored if AnaGeom=.TRUE.. Default is PathName.gplot"
574 9 pfleura2
   WRITE(*,*) ""
575 9 pfleura2
   WRITE(*,*) "*******"
576 9 pfleura2
   WRITE(*,*) "* Arrays:"
577 9 pfleura2
   WRITE(*,*) "*******"
578 9 pfleura2
   WRITE(*,*) ""
579 9 pfleura2
   WRITE(*,*) "Rcov: Array containing the covalent radii of the first 86 elements. You can modify it using,"
580 9 pfleura2
   WRITE(*,*) "    rcov(6)=0.8."
581 9 pfleura2
   WRITE(*,*) "Mass: Array containing the atomic mass of the first 86 elements."
582 9 pfleura2
   WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations. If not given, Opt'n Path will"
583 9 pfleura2
   WRITE(*,*) "    read the POTCAR file."
584 9 pfleura2
   WRITE(*,*) ""
585 9 pfleura2
   WRITE(*,*) "*******"
586 9 pfleura2
   WRITE(*,*) "* Flags:"
587 9 pfleura2
   WRITE(*,*) "*******"
588 9 pfleura2
   WRITE(*,*) ""
589 9 pfleura2
   WRITE(*,*) "MW: Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
590 9 pfleura2
   WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE."
591 9 pfleura2
   WRITE(*,*) "    unless for Coord equals CART."
592 9 pfleura2
   WRITE(*,*) ""
593 9 pfleura2
   WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
594 9 pfleura2
   WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
595 9 pfleura2
   WRITE(*,*) "CalcEProd: if TRUE the product energy will be computed. Default is FALSE. Not Compatible"
596 9 pfleura2
   WRITE(*,*) "    with RunMode=Para."
597 9 pfleura2
   WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compat-"
598 9 pfleura2
   WRITE(*,*) "    ible with RunMode=Para."
599 9 pfleura2
   WRITE(*,*) ""
600 9 pfleura2
   WRITE(*,*) "PathOnly: TRUE if one wants to generate the initial path, and stops."
601 9 pfleura2
   WRITE(*,*) "Align: If .FALSE., successive geometries along the path are not aligned on each other before"
602 9 pfleura2
   WRITE(*,*) "    path interpolation. Default is .TRUE. if there are 4 atoms or more."
603 9 pfleura2
   WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
604 9 pfleura2
   WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
605 9 pfleura2
   WRITE(*,*) ""
606 9 pfleura2
   WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of"
607 9 pfleura2
   WRITE(*,*) "    the path."
608 9 pfleura2
   WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms. If True, a &frozenlist"
609 9 pfleura2
   WRITE(*,*) "    namelist containing the list of frozen atoms must be given. If VASP is used, and frozen is"
610 9 pfleura2
   WRITE(*,*) "    not given here, the program will use the F flags of the input geometry"
611 9 pfleura2
   WRITE(*,*) "FCart: True if one wants to describe some atoms using cartesian coordinates. *** Only used in"
612 9 pfleura2
   WRITE(*,*) "    'mixed' calculations. *** If True, a &cartlist namelist containing the list of cart atoms"
613 9 pfleura2
   WRITE(*,*) "    must be given. By default, only frozen atoms are described in cartesian coordinates."
614 9 pfleura2
   WRITE(*,*) ""
615 9 pfleura2
   WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
616 9 pfleura2
   WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now."
617 9 pfleura2
   WRITE(*,*) ""
618 9 pfleura2
   WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files. See also"
619 9 pfleura2
   WRITE(*,*) "      the POSCAR option. This can be used only if prog or input=VASP."
620 9 pfleura2
   WRITE(*,*) "AnaGeom: If TRUE, Opt'n Path will create a file .dat for geometries analysis. If True, "
621 9 pfleura2
   WRITE(*,*) "Opt'n Path will look for the AnaList namelist after the Path Namelist."
622 9 pfleura2
   WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
623 9 pfleura2
   WRITE(*,*) "      If energy goes up, Smax = Smax*0.8, if not Smax = Smax*1.2. It is ensured that the"
624 9 pfleura2
   WRITE(*,*) "      dynamical Smax is within [0.5*SMax_0 ,2*Smax_0 ]."
625 9 pfleura2
   WRITE(*,*) ""
626 9 pfleura2
   WRITE(*,*) "*******"
627 9 pfleura2
   WRITE(*,*) "* Additional namelists:"
628 9 pfleura2
   WRITE(*,*) "*******"
629 9 pfleura2
   WRITE(*,*) ""
630 9 pfleura2
   WRITE(*,*) "&cartlist list= ... &end This gives the list of atoms that are described using cartesian coor-"
631 9 pfleura2
   WRITE(*,*) "      dinates. Read only if FCart=.TRUE.. To indicate an atom range, from 1 to 5 for example,"
632 9 pfleura2
   WRITE(*,*) "      one can put 1 -5 in the list. For example:"
633 9 pfleura2
   WRITE(*,*) "      &cartlist list = 1 2 6 12 -20 &end"
634 9 pfleura2
   WRITE(*,*) "      will described atoms 1, 2, 6, 12, 13, 14, 15, 16, 17, 18, 19 and 20 in cartesian."
635 9 pfleura2
   WRITE(*,*) "&Frozenlist list= ... &end This gives the list of atoms that are frozen during optimization."
636 9 pfleura2
   WRITE(*,*) "      Read only if FFrozen=.TRUE.. To indicate an atom range, from 1 to 5 for example, one"
637 9 pfleura2
   WRITE(*,*) "      can put 1 -5 in the list."
638 9 pfleura2
   WRITE(*,*) ""
639 9 pfleura2
   WRITE(*,*) "&Analist nb= ... &end list of variables for geometry analysis. If present and if Ana-"
640 9 pfleura2
   WRITE(*,*) "      Geom=T then Opt'n Path will read it and print the values of the variable in a .dat file If"
641 9 pfleura2
   WRITE(*,*) "      AnaGeom is T but Analist is not given, then Path.dat just contains the energy. Analist"
642 9 pfleura2
   WRITE(*,*) "      contains the number of variables to monitor: Nb, and is followed by the description of the"
643 9 pfleura2
   WRITE(*,*) "      variables among: b(ond) At1 At2 a(ngle) At1 At2 At3 d(ihedral) At1 At2 At3 At4 c NbAt"
644 9 pfleura2
   WRITE(*,*) "      At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added at the"
645 9 pfleura2
   WRITE(*,*) "      end of the real atoms of the system."
646 9 pfleura2
!!!!!!!!!! end of included file
647 1 pfleura2
        WRITE(*,*) ""
648 1 pfleura2
        STOP
649 1 pfleura2
     END IF
650 5 pfleura2
651 1 pfleura2
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
652 1 pfleura2
     IOOUT=6
653 1 pfleura2
  CASE (0)
654 1 pfleura2
     IOIN=5
655 1 pfleura2
     IOOUT=6
656 1 pfleura2
  END SELECT
657 1 pfleura2
658 1 pfleura2
659 1 pfleura2
  ! We initiliaze variables
660 1 pfleura2
  Pi=dacos(-1.0d0)
661 1 pfleura2
  Nat=1
662 1 pfleura2
  Fact=1.1
663 1 pfleura2
  IGeomRef=-1
664 1 pfleura2
  NGeomI=1
665 1 pfleura2
  NGeomF=8
666 1 pfleura2
  IReparam=5
667 1 pfleura2
  IReparamT=-1
668 1 pfleura2
  ISpline=5
669 1 pfleura2
  debug=.False.
670 1 pfleura2
  MW=.TRUE.
671 1 pfleura2
  bohr=.False.
672 1 pfleura2
  Coord='ZMAT'
673 1 pfleura2
  Step_method='RFO'
674 1 pfleura2
  DebugFile='Path.valid'
675 1 pfleura2
  PathName="Path"
676 1 pfleura2
  MaxCyc=50
677 1 pfleura2
  SMax=0.1D0
678 1 pfleura2
  SThresh=-1.
679 1 pfleura2
  GThresh=-1.
680 1 pfleura2
  FTan=0.
681 1 pfleura2
  Hinv=.TRUE.
682 1 pfleura2
  IniHUp=.FALSE.
683 1 pfleura2
  HupNeighbour=.FALSE.
684 1 pfleura2
  HesUpd="EMPTY"
685 5 pfleura2
  HUpdate="EMPTY"
686 1 pfleura2
  FFrozen=.FALSE.
687 1 pfleura2
  FCart=.FALSE.
688 1 pfleura2
  List=0
689 1 pfleura2
  renum=.TRUE.
690 1 pfleura2
  OptReac=.FALSE.
691 1 pfleura2
  OptProd=.FALSE.
692 1 pfleura2
  CalcEReac=.FALSE.
693 1 pfleura2
  CalcEProd=.FALSE.
694 1 pfleura2
  EReac=0.d0
695 7 pfleura2
  EProd=0.d0
696 1 pfleura2
  OptGeom=-1
697 5 pfleura2
  GeomFile="EMPTY"
698 7 pfleura2
  AnaGeom=.FALSE.
699 7 pfleura2
  AnaFile="EMPTY"
700 8 pfleura2
  GplotFile="EMPTY"
701 7 pfleura2
  Nb=0
702 7 pfleura2
  NbCom=0
703 1 pfleura2
  PathOnly=.FALSE.
704 1 pfleura2
  Prog='EMPTY'
705 1 pfleura2
  ProgExe='EMPTY'
706 1 pfleura2
  Runmode='SERIAL'
707 1 pfleura2
  Input='EMPTY'
708 1 pfleura2
  Poscar="POSCAR"
709 1 pfleura2
  AutoCart=.FALSE.
710 1 pfleura2
  VMD=.TRUE.
711 1 pfleura2
  WriteVASP=.FALSE.
712 1 pfleura2
  DynMaxStep=.FALSE.
713 1 pfleura2
  CalcName="EMPTY"
714 1 pfleura2
  ISuffix="EMPTY"
715 1 pfleura2
  OSuffix="EMPTY"
716 1 pfleura2
  NGintMax=10
717 1 pfleura2
  Align=.TRUE.
718 1 pfleura2
719 1 pfleura2
  ! Number of point used to compute the length of the path,
720 1 pfleura2
  ! and to reparameterize the path
721 1 pfleura2
  NMaxPtPath=-1
722 1 pfleura2
  FrozTol=1e-4
723 1 pfleura2
724 1 pfleura2
  ! Read the namelist
725 1 pfleura2
  read (IOIN,path)
726 1 pfleura2
727 5 pfleura2
  debug=valid("Main").or.valid("Path")
728 5 pfleura2
729 1 pfleura2
  FCalcHess=valid("CalcHess")
730 1 pfleura2
  PathName=AdjustL(PathName)
731 1 pfleura2
  Coord=AdjustL(Coord)
732 1 pfleura2
  CALL UpCase(coord)
733 1 pfleura2
734 1 pfleura2
  Step_method=AdjustL(Step_method)
735 1 pfleura2
  CALL UpCase(Step_method)
736 1 pfleura2
737 1 pfleura2
  Prog=AdjustL(Prog)
738 1 pfleura2
  CALL UpCase(prog)
739 1 pfleura2
740 1 pfleura2
  Input=AdjustL(Input)
741 1 pfleura2
  CALL UpCase(input)
742 1 pfleura2
743 1 pfleura2
  Runmode=AdjustL(Runmode)
744 1 pfleura2
  Call UpCase(Runmode)
745 1 pfleura2
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"
746 1 pfleura2
747 1 pfleura2
  if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN
748 1 pfleura2
     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"
749 1 pfleura2
     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"
750 1 pfleura2
     CalcEReac=.FALSE.
751 1 pfleura2
     CalcEProd=.FALSE.
752 1 pfleura2
  END IF
753 1 pfleura2
754 1 pfleura2
  If (IReparamT<=0) IReparamT=IReparam
755 1 pfleura2
756 1 pfleura2
! We take care of HesUpd flag
757 1 pfleura2
! this is needed because some users might still use it
758 1 pfleura2
  if (HesUpd.NE."EMPTY") THEN
759 1 pfleura2
! We are pityless:
760 1 pfleura2
     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"
761 1 pfleura2
     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"
762 1 pfleura2
     STOP
763 1 pfleura2
  END IF
764 1 pfleura2
765 5 pfleura2
  IF (HUpdate=='EMPTY') THEN
766 5 pfleura2
     IF (OptGeom>=1) THEN
767 5 pfleura2
        HupDate="BFGS"
768 5 pfleura2
     ELSE
769 5 pfleura2
        HUpdate="MS"
770 5 pfleura2
     END IF
771 5 pfleura2
  END IF
772 1 pfleura2
773 7 pfleura2
  If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN
774 5 pfleura2
     COORD='CART'
775 5 pfleura2
  END IF
776 5 pfleura2
777 1 pfleura2
  IF (COORD.EQ.'CART') THEN
778 1 pfleura2
     Renum=.FALSE.
779 5 pfleura2
     IF (OptGeom>0) THEN
780 5 pfleura2
        IGeomRef=OptGeom
781 5 pfleura2
     ELSE
782 5 pfleura2
        IGeomRef=1
783 5 pfleura2
     END IF
784 1 pfleura2
  END IF
785 1 pfleura2
786 5 pfleura2
787 1 pfleura2
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
788 5 pfleura2
  if (Prog.EQ.'2D') Prog="TEST2D"
789 3 pfleura2
  if (Prog.EQ.'TEST3') Prog="LEPS"
790 1 pfleura2
791 1 pfleura2
792 1 pfleura2
  Select Case (Prog)
793 5 pfleura2
    CASE ("EMPTY","G09","GAUSSIAN")
794 5 pfleura2
     Prog='GAUSSIAN'
795 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="g09"
796 5 pfleura2
     if (CalcName=="EMPTY") CalcName="Path"
797 5 pfleura2
     if (ISuffix=="EMPTY")  ISuffix="com"
798 5 pfleura2
     if (OSuffix=="EMPTY")  OSuffix="log"
799 8 pfleura2
     UnitProg="au"
800 8 pfleura2
     ConvE=au2kcal
801 5 pfleura2
    CASE ("G03")
802 5 pfleura2
     Prog='GAUSSIAN'
803 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="g03"
804 5 pfleura2
     if (CalcName=="EMPTY") CalcName="Path"
805 5 pfleura2
     if (ISuffix=="EMPTY")  ISuffix="com"
806 5 pfleura2
     if (OSuffix=="EMPTY")  OSuffix="log"
807 8 pfleura2
     UnitProg="au"
808 8 pfleura2
     ConvE=au2kcal
809 1 pfleura2
    CASE ("EXT")
810 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
811 1 pfleura2
     if (CalcName=="EMPTY") CalcName="Ext"
812 1 pfleura2
     if (ISuffix=="EMPTY")  ISuffix="in"
813 1 pfleura2
     if (OSuffix=="EMPTY")  OSuffix="out"
814 8 pfleura2
     UnitProg="au"
815 8 pfleura2
     ConvE=au2kcal
816 1 pfleura2
    CASE ('MOPAC')
817 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="mopac"
818 1 pfleura2
     if (CalcName=="EMPTY") CalcName="Path"
819 1 pfleura2
     if (ISuffix=="EMPTY")  ISuffix="mop"
820 1 pfleura2
     if (OSuffix=="EMPTY")  OSuffix="out"
821 8 pfleura2
     UnitProg="au"
822 8 pfleura2
     ConvE=au2kcal
823 5 pfleura2
    CASE ("SIESTA")
824 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="siesta"
825 1 pfleura2
     if (CalcName=="EMPTY") CalcName="Path"
826 5 pfleura2
     if (ISuffix=="EMPTY")  ISuffix="fdf"
827 1 pfleura2
     if (OSuffix=="EMPTY")  OSuffix="out"
828 8 pfleura2
     UnitProg="eV"
829 8 pfleura2
     ConvE=eV2kcal
830 5 pfleura2
    CASE ("VASP")
831 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="VASP"
832 8 pfleura2
     UnitProg="eV"
833 8 pfleura2
     ConvE=eV2kcal
834 5 pfleura2
    CASE ("TM","TURBOMOLE")
835 5 pfleura2
     Prog="TURBOMOLE"
836 5 pfleura2
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
837 8 pfleura2
     UnitProg="au"
838 8 pfleura2
     ConvE=au2kcal
839 5 pfleura2
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
840 5 pfleura2
     ProgExe=""
841 8 pfleura2
     UnitProg="au"
842 8 pfleura2
     ConvE=au2kcal
843 5 pfleura2
    CASE DEFAULT
844 5 pfleura2
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
845 5 pfleura2
     STOP
846 1 pfleura2
  END SELECT
847 1 pfleura2
848 1 pfleura2
  if (Input.EQ.'EMPTY') THEN
849 1 pfleura2
     Select Case (Prog)
850 1 pfleura2
       CASE ("VASP")
851 1 pfleura2
          Input=Prog
852 1 pfleura2
       CASE ("CHAMFRE")
853 1 pfleura2
          Input="VASP"
854 1 pfleura2
       CASE DEFAULT
855 1 pfleura2
          Input='XYZ'
856 1 pfleura2
     END SELECT
857 1 pfleura2
    WRITE(*,*) "Input has been set to the default: ",INPUT
858 1 pfleura2
 END IF
859 1 pfleura2
860 1 pfleura2
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
861 1 pfleura2
862 1 pfleura2
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
863 1 pfleura2
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
864 1 pfleura2
     STOP
865 1 pfleura2
  END IF
866 1 pfleura2
867 1 pfleura2
  if (OptGeom.GE.1) THEN
868 5 pfleura2
     FPrintGeom=.FALSE.
869 5 pfleura2
     If (GeomFile/='EMPTY') THEN
870 5 pfleura2
        FPrintGeom=.TRUE.
871 5 pfleura2
     END IF
872 1 pfleura2
     IF (IGeomRef.LE.0) THEN
873 1 pfleura2
        IGeomRef=OptGeom
874 1 pfleura2
     ELSE
875 1 pfleura2
        IF (IGeomRef.NE.OptGeom) THEN
876 1 pfleura2
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
877 1 pfleura2
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
878 1 pfleura2
879 1 pfleura2
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
880 1 pfleura2
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
881 1 pfleura2
           STOP
882 1 pfleura2
        END IF
883 1 pfleura2
     END IF
884 1 pfleura2
  END IF
885 1 pfleura2
886 5 pfleura2
  ALLOCATE(Cart(Nat))
887 5 pfleura2
  Cart=0
888 1 pfleura2
889 5 pfleura2
  ALLOCATE(Frozen(Nat))
890 5 pfleura2
  Frozen=0
891 1 pfleura2
892 5 pfleura2
  IF (FCart) THEN
893 5 pfleura2
! We rewind the file to be sure to browse all of it to read the Cart namelist
894 5 pfleura2
     REWIND(IOIN)
895 5 pfleura2
     List=0
896 5 pfleura2
     READ(IOIN,cartlist)
897 5 pfleura2
     IF (COORD.NE.'CART') THEN
898 5 pfleura2
        Cart=List(1:Nat)
899 5 pfleura2
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
900 1 pfleura2
     ELSE
901 5 pfleura2
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
902 5 pfleura2
        WRITE(*,*) "I will discard the cartlist"
903 1 pfleura2
        Cart=0
904 1 pfleura2
     END IF
905 1 pfleura2
  END IF
906 1 pfleura2
907 5 pfleura2
  if (FFrozen) THEN
908 1 pfleura2
909 5 pfleura2
     REWIND(IOIN)
910 5 pfleura2
     List=0
911 5 pfleura2
     READ(IOIN,Frozenlist)
912 5 pfleura2
     Frozen=List(1:Nat)
913 5 pfleura2
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
914 5 pfleura2
  END IF
915 5 pfleura2
916 1 pfleura2
  If (FCart) Call Expand(Cart,NCart,Nat)
917 1 pfleura2
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
918 1 pfleura2
919 1 pfleura2
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
920 1 pfleura2
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
921 1 pfleura2
922 7 pfleura2
  if (AnaGeom) THEN
923 7 pfleura2
     REWIND(IOIN)
924 7 pfleura2
     READ(IOIN,AnaList,IOSTAT=IoS)
925 7 pfleura2
     IF (IOS==0) THEN ! we have a AnaList namelist
926 7 pfleura2
        Call ReadAnaList
927 8 pfleura2
     ELSE
928 8 pfleura2
! No AnaList namelist, we will print only the energy
929 8 pfleura2
        FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)'
930 7 pfleura2
     END IF
931 7 pfleura2
     IF (AnaFile=="EMPTY") THEN
932 7 pfleura2
        AnaFile=TRIM(PathName) // '.dat'
933 7 pfleura2
     END IF
934 7 pfleura2
     OPEN(IODat,File=AnaFile)
935 8 pfleura2
     IF (GplotFile=="EMPTY") THEN
936 8 pfleura2
        GplotFile=TRIM(PathName) // '.gplot'
937 8 pfleura2
     END IF
938 8 pfleura2
939 8 pfleura2
     OPEN(IODat,File=AnaFile)
940 8 pfleura2
     If (OptGeom<=0) THEN
941 8 pfleura2
        OPEN(IOGplot,File=GPlotFile)
942 8 pfleura2
        WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist"
943 8 pfleura2
        WRITE(IoGplot,'(A)') " set pointsize 2"
944 8 pfleura2
        WRITE(IoGplot,'(A)') " xcol=1"
945 8 pfleura2
        WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3
946 8 pfleura2
     END IF
947 7 pfleura2
  END IF
948 7 pfleura2
949 7 pfleura2
950 7 pfleura2
951 1 pfleura2
  if (NMaxPtPath.EQ.-1) then
952 1 pfleura2
     NMaxPtPath=min(50*NGeomF,2000)
953 1 pfleura2
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
954 1 pfleura2
  end if
955 1 pfleura2
956 1 pfleura2
  !  NMaxPtPath=10000
957 1 pfleura2
958 1 pfleura2
  ! We ensure that FTan is in [0.:1.]
959 1 pfleura2
  FTan=min(abs(FTan),1.d0)
960 1 pfleura2
961 1 pfleura2
! PFL 2011 Mar 14 ->
962 1 pfleura2
! Added some printing for analyses with Anapath
963 1 pfleura2
  if (debug) THEN
964 1 pfleura2
     WRITE(IOOUT,path)
965 1 pfleura2
  ELSE
966 1 pfleura2
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC
967 1 pfleura2
! and PATHNAME
968 1 pfleura2
     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat
969 1 pfleura2
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI
970 1 pfleura2
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF
971 1 pfleura2
     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC
972 1 pfleura2
     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)
973 1 pfleura2
  END IF
974 1 pfleura2
975 1 pfleura2
  FTan=1.-FTan
976 1 pfleura2
977 1 pfleura2
  !Thresholds...
978 1 pfleura2
  if (SThresh.LE.0) SThresh=0.01d0
979 1 pfleura2
  If (GThresh.LE.0) GThresh=SThresh/4.
980 1 pfleura2
  SThreshRms=Sthresh*2./3.
981 1 pfleura2
  GThreshRms=Gthresh*2./3.
982 1 pfleura2
983 1 pfleura2
  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
984 1 pfleura2
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
985 1 pfleura2
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
986 1 pfleura2
  ALLOCATE(MassAt(NAt),OrderInv(Nat))
987 1 pfleura2
988 7 pfleura2
  AtName=""
989 7 pfleura2
  MassAt=1.
990 7 pfleura2
991 8 pfleura2
! As we have used rewind many times, the position in the input file
992 8 pfleura2
! is blurry (at least !) we thus have to go back to the Geometries description
993 8 pfleura2
! (TO DO: use a better Input parser !)
994 8 pfleura2
995 8 pfleura2
  REWIND(IOIN)
996 8 pfleura2
! We search the '&path' namelist
997 8 pfleura2
  READ(IOIN,'(A)',Iostat=Ios) Line
998 8 pfleura2
  Fini=.FALSE.
999 8 pfleura2
  DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0)
1000 8 pfleura2
     if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line)
1001 8 pfleura2
     READ(IOIN,'(A)',Iostat=Ios) Line
1002 8 pfleura2
  END DO
1003 8 pfleura2
  if (debug) WRITE(*,*) "Path found:",TRIM(LINE)
1004 8 pfleura2
  ! We are on the &path line, we move to the end
1005 8 pfleura2
  Call NoComment(Line)
1006 8 pfleura2
  DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
1007 8 pfleura2
     if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
1008 8 pfleura2
     READ(IOIN,'(A)',Iostat=Ios) Line
1009 8 pfleura2
     Call NoComment(Line)
1010 8 pfleura2
     Call noString(Line)
1011 8 pfleura2
  END DO
1012 8 pfleura2
! We are at then end of &Path / namelist
1013 8 pfleura2
! We scan for other Namelists
1014 8 pfleura2
  READ(IOIN,'(A)',Iostat=Ios) Line
1015 8 pfleura2
  if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
1016 8 pfleura2
  DO WHILE  ((IoS==0).AND.(Index(Line,'&')/=0))
1017 8 pfleura2
     Call NoComment(Line)
1018 8 pfleura2
     if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
1019 8 pfleura2
     Idx=InString(Line,'Analist')
1020 8 pfleura2
     DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
1021 8 pfleura2
        if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
1022 8 pfleura2
        READ(IOIN,'(A)',Iostat=Ios) Line
1023 8 pfleura2
        Call NoComment(Line)
1024 8 pfleura2
     END DO
1025 8 pfleura2
     If (Idx/=0) THEN
1026 8 pfleura2
! we have just read &Analist namelist, we have to skip the variables description
1027 8 pfleura2
        DO I=1,Nb
1028 8 pfleura2
           READ(IOIN,'(A)',Iostat=Ios) Line
1029 8 pfleura2
           if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE)
1030 8 pfleura2
        END DO
1031 8 pfleura2
     END IF
1032 8 pfleura2
     READ(IOIN,'(A)',Iostat=Ios) Line
1033 8 pfleura2
     if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE)
1034 8 pfleura2
  END DO
1035 8 pfleura2
  BACKSPACE(IOIN)
1036 8 pfleura2
1037 1 pfleura2
  ! We read the initial geometries
1038 1 pfleura2
  Call Read_Geom(input)
1039 1 pfleura2
1040 1 pfleura2
  ! We convert atome names into atom mass number
1041 1 pfleura2
  DO I=1,NAt
1042 1 pfleura2
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
1043 1 pfleura2
     Atome(I)=ConvertNumAt(AtName(I))
1044 1 pfleura2
  END DO
1045 1 pfleura2
1046 7 pfleura2
 If (AnaGeom) THEN
1047 8 pfleura2
    If (NbVar>0) THEN
1048 7 pfleura2
! We print the description of the varialbes in AnaFile
1049 8 pfleura2
       Call PrintAnaList(IoDat)
1050 8 pfleura2
       If (OptGeom<=0) THEN
1051 8 pfleura2
          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
1052 8 pfleura2
       ELSE
1053 8 pfleura2
          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
1054 8 pfleura2
       END IF
1055 8 pfleura2
    ELSE
1056 8 pfleura2
       If (OptGeom<=0) THEN
1057 8 pfleura2
          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
1058 8 pfleura2
       ELSE
1059 8 pfleura2
          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
1060 8 pfleura2
       END IF
1061 8 pfleura2
    END IF
1062 7 pfleura2
 END IF
1063 1 pfleura2
1064 1 pfleura2
  ! PFL 23 Sep 2008 ->
1065 1 pfleura2
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
1066 1 pfleura2
  ! between the first geometry and the others.
1067 1 pfleura2
  IF (NFroz.GT.0) THEN
1068 1 pfleura2
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
1069 1 pfleura2
     ALLOCATE(x2(nat),y2(nat),z2(nat))
1070 1 pfleura2
     ListAtFroz=.FALSE.
1071 1 pfleura2
1072 1 pfleura2
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
1073 1 pfleura2
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
1074 1 pfleura2
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
1075 1 pfleura2
1076 1 pfleura2
     SELECT CASE (NFroz)
1077 1 pfleura2
        ! We have Frozen atoms
1078 1 pfleura2
        ! PFL 28 Oct 2008 ->
1079 1 pfleura2
        ! It might happen that the geometries are translated/rotated.
1080 1 pfleura2
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
1081 1 pfleura2
        ! three frozen atoms.
1082 1 pfleura2
1083 1 pfleura2
     CASE (3:)
1084 1 pfleura2
        DO I=1,NFroz
1085 1 pfleura2
           ListAtFroz(Frozen(I))=.TRUE.
1086 1 pfleura2
        END DO
1087 1 pfleura2
     CASE (2)
1088 1 pfleura2
        IAt1=Frozen(1)
1089 1 pfleura2
        IAt2=Frozen(2)
1090 1 pfleura2
        ListAtFroz(Iat1)=.TRUE.
1091 1 pfleura2
        ListAtFroz(Iat2)=.TRUE.
1092 1 pfleura2
        x2(1)=x1(Iat2)-x1(Iat1)
1093 1 pfleura2
        y2(1)=y1(Iat2)-y1(Iat1)
1094 1 pfleura2
        z2(1)=z1(Iat2)-z1(Iat1)
1095 1 pfleura2
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1096 1 pfleura2
        Dist=-1.
1097 1 pfleura2
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1098 1 pfleura2
        DO I=1,Nat
1099 1 pfleura2
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
1100 1 pfleura2
           x2(2)=x1(Iat2)-x1(I)
1101 1 pfleura2
           y2(2)=y1(Iat2)-y1(I)
1102 1 pfleura2
           z2(2)=z1(Iat2)-z1(I)
1103 1 pfleura2
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
1104 1 pfleura2
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
1105 1 pfleura2
           if (abs(ca)<0.9) THEN
1106 1 pfleura2
              IF (Norm2>Dist) THEN
1107 1 pfleura2
                 Iat3=I
1108 1 pfleura2
                 Dist=Norm2
1109 1 pfleura2
              END IF
1110 1 pfleura2
           END IF
1111 1 pfleura2
        END DO
1112 1 pfleura2
        ListAtFroz(Iat3)=.TRUE.
1113 1 pfleura2
        if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3
1114 1 pfleura2
     CASE (1)
1115 1 pfleura2
        IAt1=Frozen(1)
1116 1 pfleura2
        ListAtFroz(Iat1)=.TRUE.
1117 1 pfleura2
        Dist=-1.
1118 1 pfleura2
        ! We scan all atoms to find one that is further from At1
1119 1 pfleura2
        DO I=1,Nat
1120 1 pfleura2
           if (I.EQ.Iat1) Cycle
1121 1 pfleura2
           x2(1)=x1(Iat1)-x1(I)
1122 1 pfleura2
           y2(1)=y1(Iat1)-y1(I)
1123 1 pfleura2
           z2(1)=z1(Iat1)-z1(I)
1124 1 pfleura2
           Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1125 1 pfleura2
           IF (Norm1>Dist) THEN
1126 1 pfleura2
              Iat2=I
1127 1 pfleura2
              Dist=Norm1
1128 1 pfleura2
           END IF
1129 1 pfleura2
        END DO
1130 1 pfleura2
        ListAtFroz(Iat2)=.TRUE.
1131 1 pfleura2
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
1132 1 pfleura2
        x2(1)=x1(Iat2)-x1(Iat1)
1133 1 pfleura2
        y2(1)=y1(Iat2)-y1(Iat1)
1134 1 pfleura2
        z2(1)=z1(Iat2)-z1(Iat1)
1135 1 pfleura2
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1136 1 pfleura2
        Dist=-1.
1137 1 pfleura2
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1138 8 pfleura2
        Iat3=1
1139 8 pfleura2
        DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat))
1140 8 pfleura2
           Iat3=Iat3+1
1141 8 pfleura2
        END DO
1142 1 pfleura2
        DO I=1,Nat
1143 1 pfleura2
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
1144 1 pfleura2
           x2(2)=x1(Iat2)-x1(I)
1145 1 pfleura2
           y2(2)=y1(Iat2)-y1(I)
1146 1 pfleura2
           z2(2)=z1(Iat2)-z1(I)
1147 1 pfleura2
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
1148 1 pfleura2
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
1149 1 pfleura2
           if (abs(ca)<0.9) THEN
1150 1 pfleura2
              IF (Norm2>Dist) THEN
1151 1 pfleura2
                 Iat3=I
1152 1 pfleura2
                 Dist=Norm2
1153 1 pfleura2
              END IF
1154 1 pfleura2
           END IF
1155 1 pfleura2
        END DO
1156 1 pfleura2
        ListAtFroz(Iat3)=.TRUE.
1157 1 pfleura2
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3
1158 1 pfleura2
     END SELECT
1159 1 pfleura2
1160 1 pfleura2
     DO I=2,NGeomI
1161 1 pfleura2
        ! First, we align the Ith geometry on the first one, using only
1162 1 pfleura2
        ! the positions of the "frozen" atoms. (see above: it might be that
1163 8 pfleura2
        ! ListAtFroz contains non frozen atoms used to align the geometries
1164 1 pfleura2
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
1165 1 pfleura2
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
1166 1 pfleura2
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
1167 1 pfleura2
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
1168 1 pfleura2
1169 1 pfleura2
        FChkFrozen=.FALSE.
1170 1 pfleura2
        DO J=1,NFroz
1171 1 pfleura2
           IAt=Frozen(J)
1172 1 pfleura2
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
1173 1 pfleura2
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
1174 1 pfleura2
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
1175 1 pfleura2
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
1176 1 pfleura2
              FChkFrozen=.TRUE.
1177 1 pfleura2
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
1178 1 pfleura2
              WRITE(*,*) "from geometry 1 to geometry ",I
1179 1 pfleura2
           END IF
1180 1 pfleura2
        END DO
1181 1 pfleura2
     END DO
1182 1 pfleura2
     IF (FChkFrozen) THEN
1183 1 pfleura2
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
1184 1 pfleura2
        STOP
1185 1 pfleura2
     END IF
1186 1 pfleura2
  END IF
1187 1 pfleura2
1188 1 pfleura2
1189 1 pfleura2
  N3at=3*Nat
1190 1 pfleura2
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
1191 1 pfleura2
1192 1 pfleura2
1193 5 pfleura2
  Call ReadInput
1194 1 pfleura2
1195 1 pfleura2
  IF (COORD=="MIXED") Call TestCart(AutoCart)
1196 1 pfleura2
1197 1 pfleura2
  SELECT CASE(NFroz)
1198 1 pfleura2
  CASE (2)
1199 1 pfleura2
     IntFroz=1
1200 1 pfleura2
  CASE (3)
1201 1 pfleura2
     IntFroz=3
1202 1 pfleura2
  CASE (4:)
1203 1 pfleura2
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
1204 1 pfleura2
  END SELECT
1205 1 pfleura2
1206 1 pfleura2
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
1207 1 pfleura2
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
1208 1 pfleura2
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
1209 1 pfleura2
1210 1 pfleura2
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
1211 1 pfleura2
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
1212 1 pfleura2
     COORD="ZMAT"
1213 1 pfleura2
  END IF
1214 1 pfleura2
1215 1 pfleura2
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
1216 1 pfleura2
1217 1 pfleura2
  SELECT CASE(COORD)
1218 1 pfleura2
  CASE ('ZMAT')
1219 1 pfleura2
     NCoord=Nfree
1220 1 pfleura2
     ALLOCATE(dzdc(3,nat,3,nat))
1221 1 pfleura2
  CASE ('MIXED')
1222 1 pfleura2
     SELECT CASE (NCart+NFroz)
1223 1 pfleura2
     CASE (1)
1224 1 pfleura2
        NCoord=N3At-3
1225 1 pfleura2
     CASE (2)
1226 1 pfleura2
        NCoord=N3At-1
1227 1 pfleura2
     CASE (3:)
1228 1 pfleura2
        NCoord=N3At
1229 1 pfleura2
     END SELECT
1230 1 pfleura2
     ALLOCATE(dzdc(3,nat,3,nat))
1231 1 pfleura2
  CASE ('BAKER')
1232 1 pfleura2
     Symmetry_elimination=0
1233 1 pfleura2
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
1234 1 pfleura2
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
1235 1 pfleura2
     ALLOCATE(BTransInv(NCoord,3*Nat))
1236 1 pfleura2
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
1237 1 pfleura2
  CASE ('HYBRID')
1238 1 pfleura2
     NCoord=N3at
1239 1 pfleura2
  CASE ('CART')
1240 1 pfleura2
     NCoord=N3at
1241 1 pfleura2
  END SELECT
1242 1 pfleura2
1243 1 pfleura2
  if (debug) WRITE(*,*) "DBG Path, L826, Coord, NCart,NCoord, N3At",TRIM(COORD),Ncart, NCoord,N3at
1244 1 pfleura2
1245 1 pfleura2
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
1246 1 pfleura2
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
1247 1 pfleura2
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
1248 1 pfleura2
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
1249 1 pfleura2
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
1250 1 pfleura2
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
1251 1 pfleura2
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
1252 1 pfleura2
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
1253 1 pfleura2
1254 1 pfleura2
  ZeroVec=0.d0
1255 1 pfleura2
  DO I =1, NGeomF
1256 1 pfleura2
     IntTangent(I,:)=0.d0
1257 1 pfleura2
  END DO
1258 1 pfleura2
1259 1 pfleura2
  if (debug) THEN
1260 1 pfleura2
     Print *, 'L886, IntTangent(1,:)='
1261 1 pfleura2
     Print *, IntTangent(1,:)
1262 1 pfleura2
  END IF
1263 1 pfleura2
1264 1 pfleura2
  if (.NOT.OptReac) Energies(1)=Ereac
1265 1 pfleura2
  if (.NOT.OptProd) Energies(NGeomF)=EProd
1266 1 pfleura2
  MaxStep=SMax
1267 1 pfleura2
1268 1 pfleura2
  ! We initialize the mass array
1269 1 pfleura2
  if (MW) THEN
1270 1 pfleura2
     WRITE(*,*) 'Working in MW coordinates'
1271 1 pfleura2
     DO I=1,Nat
1272 1 pfleura2
        massAt(I)=Mass(Atome(I))
1273 1 pfleura2
     END DO
1274 1 pfleura2
  ELSE
1275 1 pfleura2
     DO I=1,Nat
1276 1 pfleura2
        MassAt(I)=1.d0
1277 1 pfleura2
     END DO
1278 1 pfleura2
  END IF
1279 1 pfleura2
1280 1 pfleura2
  WRITE(*,*) "Prog=",TRIM(PROG)
1281 1 pfleura2
1282 1 pfleura2
  ALLOCATE(FrozAtoms(Nat))
1283 1 pfleura2
1284 1 pfleura2
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
1285 1 pfleura2
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
1286 1 pfleura2
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
1287 1 pfleura2
  IF (NFroz.EQ.0) THEN
1288 1 pfleura2
     FrozAtoms=.TRUE.
1289 1 pfleura2
  ELSE
1290 1 pfleura2
     I=1
1291 1 pfleura2
     NFr=0
1292 1 pfleura2
     FrozAtoms=.False.
1293 1 pfleura2
     DO WHILE (Frozen(I).NE.0)
1294 1 pfleura2
        IF (Frozen(I).LT.0) THEN
1295 1 pfleura2
           DO J=Frozen(I-1),abs(Frozen(I))
1296 1 pfleura2
              IF (.NOT.FrozAtoms(J)) THEN
1297 1 pfleura2
                 NFr=NFr+1
1298 1 pfleura2
                 FrozAtoms(J)=.TRUE.
1299 1 pfleura2
              END IF
1300 1 pfleura2
           END DO
1301 1 pfleura2
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1302 1 pfleura2
           FrozAtoms(Frozen(I))=.TRUE.
1303 1 pfleura2
           NFr=NFr+1
1304 1 pfleura2
        END IF
1305 1 pfleura2
        I=I+1
1306 1 pfleura2
     END DO
1307 1 pfleura2
     IF (NFr.NE.NFroz) THEN
1308 1 pfleura2
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1309 1 pfleura2
        STOP
1310 1 pfleura2
     END IF
1311 1 pfleura2
  END IF
1312 1 pfleura2
1313 5 pfleura2
!  if (debug) THEN
1314 1 pfleura2
     !Print *, 'L968, IntTangent(1,:)='
1315 1 pfleura2
     !Print *, IntTangent(1,:)
1316 5 pfleura2
!  END IF
1317 1 pfleura2
1318 1 pfleura2
  ! We have read everything we needed... time to work :-)
1319 1 pfleura2
  IOpt=0
1320 1 pfleura2
  FirstTimePathCreate = .TRUE. ! Initialized.
1321 1 pfleura2
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1322 1 pfleura2
  FirstTimePathCreate = .FALSE.
1323 1 pfleura2
1324 1 pfleura2
  if (debug) THEN
1325 1 pfleura2
     Print *, 'L980, IntTangent(1,:)='
1326 1 pfleura2
     Print *, IntTangent(1,:)
1327 1 pfleura2
  END IF
1328 1 pfleura2
1329 1 pfleura2
  ! PFL 30 Mar 2008 ->
1330 1 pfleura2
  ! The following is now allocated in Calc_Baker. This is more logical to me
1331 1 pfleura2
  !   IF (COORD .EQ. "BAKER") Then
1332 1 pfleura2
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1333 1 pfleura2
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1334 1 pfleura2
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1335 1 pfleura2
  !   END IF
1336 1 pfleura2
  ! <- PFL 30 mar 2008
1337 1 pfleura2
1338 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1339 1 pfleura2
  !
1340 1 pfleura2
  !
1341 1 pfleura2
  ! OptGeom is the index of the geometry which is to be optimized.
1342 1 pfleura2
  IF (Optgeom.GT.0) THEN
1343 1 pfleura2
     Flag_Opt_Geom=.TRUE.
1344 1 pfleura2
     SELECT CASE(Coord)
1345 1 pfleura2
     CASE ('CART','HYBRID')
1346 1 pfleura2
!!! CAUTION : PFL 29.AUG.2008 ->
1347 1 pfleura2
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1348 1 pfleura2
        ! This should be made  consistent !!!!!!!!!!!!!!!
1349 5 pfleura2
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1350 1 pfleura2
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1351 1 pfleura2
!!! <- CAUTION : PFL 29.AUG.2008
1352 1 pfleura2
     CASE ('ZMAT','MIXED')
1353 1 pfleura2
        !GeomTmp=IntCoordF(OptGeom,:)
1354 1 pfleura2
        GeomTmp=IntCoordI(OptGeom,:)
1355 1 pfleura2
     CASE ('BAKER')
1356 1 pfleura2
        !GeomTmp=IntCoordF(OptGeom,:)
1357 1 pfleura2
        GeomTmp=IntCoordI(OptGeom,:)
1358 1 pfleura2
     CASE DEFAULT
1359 1 pfleura2
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1360 1 pfleura2
        STOP
1361 1 pfleura2
     END SELECT
1362 1 pfleura2
1363 1 pfleura2
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1364 1 pfleura2
     Flag_Opt_Geom = .TRUE.
1365 1 pfleura2
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1366 1 pfleura2
1367 1 pfleura2
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1368 9 pfleura2
     if (debug) WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1369 1 pfleura2
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1370 1 pfleura2
     STOP
1371 1 pfleura2
  END IF ! IF (Optgeom.GT.0) THEN
1372 1 pfleura2
1373 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1374 1 pfleura2
  ! End of GeomOpt
1375 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1376 1 pfleura2
1377 1 pfleura2
  IF (PathOnly) THEN
1378 9 pfleura2
     Call Header("PathOnly=.T. , Stop here")
1379 1 pfleura2
     Call Write_path(-1)
1380 1 pfleura2
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1381 1 pfleura2
     STOP
1382 1 pfleura2
  END IF
1383 1 pfleura2
1384 1 pfleura2
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1385 1 pfleura2
1386 1 pfleura2
  DEALLOCATE(XyzGeomI,IntCoordI)
1387 1 pfleura2
  NGeomI=NGeomF
1388 1 pfleura2
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1389 1 pfleura2
1390 1 pfleura2
  IF (DEBUG) WRITE(*,*) ' Back from PathCreate, L965, Path.f90'
1391 1 pfleura2
  ! we print this path, which is the first one :-) called Iteration0
1392 1 pfleura2
  ! we now have to calculate the gradient for each point (and the energy
1393 1 pfleura2
  ! if we work at 0K)
1394 1 pfleura2
1395 1 pfleura2
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1396 1 pfleura2
     DO IGeom=1,NGeomF
1397 1 pfleura2
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1398 1 pfleura2
        DO I=1,3*Nat-6
1399 1 pfleura2
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1400 1 pfleura2
        END DO
1401 1 pfleura2
     END DO
1402 1 pfleura2
  END IF
1403 1 pfleura2
1404 1 pfleura2
  ! To calculate B^T^-1 for all extrapolated final geometries:
1405 1 pfleura2
  ! PFL 18 Jan 2008 ->
1406 1 pfleura2
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1407 1 pfleura2
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1408 1 pfleura2
  ! <-- PFL 18 Jan 2008
1409 1 pfleura2
1410 1 pfleura2
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1411 1 pfleura2
     DO IGeom=1,NGeomF
1412 1 pfleura2
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
1413 1 pfleura2
        DO I=1,3*Nat-6
1414 1 pfleura2
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1415 1 pfleura2
        END DO
1416 1 pfleura2
     END DO
1417 1 pfleura2
  END IF
1418 1 pfleura2
1419 1 pfleura2
  IOpt=0
1420 2 pfleura2
  Call EgradPath(IOpt)
1421 1 pfleura2
1422 1 pfleura2
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1423 1 pfleura2
  Call Write_path(IOpt)
1424 7 pfleura2
! Shall we analyze the geometries ?
1425 7 pfleura2
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1426 7 pfleura2
1427 1 pfleura2
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1428 1 pfleura2
1429 1 pfleura2
  ! We have the geometries, the first gradients... we will generate the first hessian matrices :-)
1430 1 pfleura2
  ! ... v3.71 Only if IniHUp=.TRUE. which is not the default. In this version, we take the hessian
1431 1 pfleura2
  ! of end points as Unity matrices, and we update them along the path.
1432 1 pfleura2
1433 1 pfleura2
  ! By default, Hess=Id
1434 1 pfleura2
  Hess=0.
1435 1 pfleura2
  IF(Coord .EQ. "ZMAT") Then
1436 1 pfleura2
     ! We use the fact that we know that approximate good hessian values
1437 1 pfleura2
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1438 1 pfleura2
     DO IGeom=1,NGeomF
1439 1 pfleura2
        Hess(IGeom,1,1)=0.5d0
1440 1 pfleura2
        Hess(IGeom,2,2)=0.5d0
1441 1 pfleura2
        Hess(IGeom,3,3)=0.1d0
1442 1 pfleura2
        DO J=1,Nat-3
1443 1 pfleura2
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1444 1 pfleura2
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1445 1 pfleura2
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1446 1 pfleura2
        END DO
1447 1 pfleura2
        IF (HInv) THEN
1448 1 pfleura2
           DO I=1,NCoord
1449 1 pfleura2
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1450 1 pfleura2
           END DO
1451 1 pfleura2
        END IF
1452 1 pfleura2
     END DO
1453 1 pfleura2
  ELSE
1454 1 pfleura2
     DO I=1,NCoord
1455 1 pfleura2
        DO J=1,NGeomF
1456 1 pfleura2
           Hess(J,I,I)=1.d0
1457 1 pfleura2
        END DO
1458 1 pfleura2
     END DO
1459 1 pfleura2
  END IF
1460 1 pfleura2
1461 1 pfleura2
  IF (COORD .EQ. "BAKER") THEN
1462 1 pfleura2
     ! UMat(NPrim,NCoord)
1463 1 pfleura2
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1464 1 pfleura2
     Hprim=0.d0
1465 1 pfleura2
     ScanCoord=>Coordinate
1466 1 pfleura2
     I=0
1467 1 pfleura2
     DO WHILE (Associated(ScanCoord%next))
1468 1 pfleura2
        I=I+1
1469 1 pfleura2
        SELECT CASE (ScanCoord%Type)
1470 1 pfleura2
        CASE ('BOND')
1471 1 pfleura2
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1472 1 pfleura2
           Hprim(I,I) = 0.5d0
1473 1 pfleura2
           !END DO
1474 1 pfleura2
        CASE ('ANGLE')
1475 1 pfleura2
           Hprim(I,I) = 0.2d0
1476 1 pfleura2
        CASE ('DIHEDRAL')
1477 1 pfleura2
           Hprim(I,I) = 0.1d0
1478 1 pfleura2
        END SELECT
1479 1 pfleura2
        ScanCoord => ScanCoord%next
1480 1 pfleura2
     END DO
1481 1 pfleura2
1482 1 pfleura2
     ! HprimU = H^prim U:
1483 1 pfleura2
     HprimU=0.d0
1484 1 pfleura2
     DO I=1, 3*Nat-6
1485 1 pfleura2
        DO J=1,NPrim
1486 1 pfleura2
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1487 1 pfleura2
        END DO
1488 1 pfleura2
     END DO
1489 1 pfleura2
1490 1 pfleura2
     ! Hess = U^T HprimU = U^T H^prim U:
1491 1 pfleura2
     Hess=0.d0
1492 1 pfleura2
     DO I=1, 3*Nat-6
1493 1 pfleura2
        DO J=1,NPrim
1494 1 pfleura2
           ! UMat^T is needed below.
1495 1 pfleura2
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1496 1 pfleura2
        END DO
1497 1 pfleura2
     END DO
1498 1 pfleura2
     DO K=2,NGeomF
1499 1 pfleura2
        Hess(K,:,:)=Hess(1,:,:)
1500 1 pfleura2
     END DO
1501 1 pfleura2
     DEALLOCATE(Hprim,HprimU)
1502 1 pfleura2
  end if !  ends IF COORD=BAKER
1503 1 pfleura2
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1504 1 pfleura2
  ALLOCATE(HessTmp(NCoord,NCoord))
1505 1 pfleura2
1506 1 pfleura2
  IF (IniHUp) THEN
1507 1 pfleura2
     IF (FCalcHess) THEN
1508 1 pfleura2
        ! The numerical calculation of the Hessian has been switched on
1509 1 pfleura2
        DO IGeom=1,NGeomF
1510 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1511 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
1512 1 pfleura2
           ELSE
1513 1 pfleura2
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1514 1 pfleura2
           END IF
1515 1 pfleura2
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1516 1 pfleura2
        END DO
1517 1 pfleura2
     ELSE
1518 1 pfleura2
        ! We generate a forward hessian and a backward one and we mix them.
1519 1 pfleura2
        ! First, the forward way:
1520 1 pfleura2
        DO IGeom=2,NGeomF-1
1521 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1522 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1523 1 pfleura2
           ELSE
1524 1 pfleura2
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1525 1 pfleura2
           END IF
1526 1 pfleura2
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1527 1 pfleura2
           HessTmp=Hess(IGeom-1,:,:)
1528 1 pfleura2
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1529 1 pfleura2
           Hess(IGeom,:,:)=HessTmp
1530 1 pfleura2
        END DO
1531 1 pfleura2
1532 1 pfleura2
        ! Now, the backward way:
1533 1 pfleura2
        HessTmp=Hess(NGeomF,:,:)
1534 1 pfleura2
        DO IGeom=NGeomF-1,2,-1
1535 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1536 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1537 1 pfleura2
           ELSE
1538 1 pfleura2
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1539 1 pfleura2
           END IF
1540 1 pfleura2
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1541 1 pfleura2
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1542 1 pfleura2
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1543 1 pfleura2
           ca=cos(alpha)
1544 1 pfleura2
           sa=sin(alpha)
1545 1 pfleura2
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1546 1 pfleura2
        END DO
1547 1 pfleura2
     END IF ! matches IF (FCalcHess)
1548 1 pfleura2
  END IF ! matches IF (IniHUp) THEN
1549 1 pfleura2
1550 1 pfleura2
  ! For debug purposes, we diagonalize the Hessian matrices
1551 1 pfleura2
  IF (debug) THEN
1552 1 pfleura2
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1553 1 pfleura2
     DO I=1,NGeomF
1554 1 pfleura2
        WRITE(*,*) "DBG Main - Point ",I
1555 1 pfleura2
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1556 1 pfleura2
        DO J=1,NCoord
1557 1 pfleura2
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1558 1 pfleura2
        END DO
1559 1 pfleura2
     END DO
1560 1 pfleura2
  END IF ! matches IF (debug) THEN
1561 1 pfleura2
1562 1 pfleura2
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1563 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1564 1 pfleura2
  !
1565 1 pfleura2
  ! Main LOOP : Optimization of the Path
1566 1 pfleura2
  !
1567 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1568 1 pfleura2
  if (debug)   Print *, 'NGeomF=', NGeomF
1569 1 pfleura2
  allocation_flag = .TRUE.
1570 1 pfleura2
1571 1 pfleura2
  Fini=.FALSE.
1572 1 pfleura2
1573 1 pfleura2
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1574 1 pfleura2
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1575 1 pfleura2
     IOpt=IOpt+1
1576 1 pfleura2
1577 1 pfleura2
     SELECT CASE (COORD)
1578 1 pfleura2
     CASE ('ZMAT','MIXED','BAKER')
1579 1 pfleura2
        GeomOld_all=IntCoordF
1580 1 pfleura2
     CASE ('CART','HYBRID')
1581 1 pfleura2
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1582 1 pfleura2
     CASE DEFAULT
1583 1 pfleura2
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1584 1 pfleura2
        STOP
1585 1 pfleura2
     END SELECT
1586 1 pfleura2
1587 1 pfleura2
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1588 1 pfleura2
          valid("printtangent")) THEN
1589 1 pfleura2
        WRITE(*,*) "Visualization of tangents"
1590 1 pfleura2
        Idx=min(12,NCoord)
1591 1 pfleura2
        DO I=1,NGeomF
1592 1 pfleura2
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1593 1 pfleura2
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1594 1 pfleura2
           WRITE(*,*)
1595 1 pfleura2
        END DO
1596 1 pfleura2
        WRITE(*,*) "END of tangents"
1597 1 pfleura2
     END IF
1598 1 pfleura2
1599 1 pfleura2
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1600 1 pfleura2
          valid("printgrad")) THEN
1601 1 pfleura2
        WRITE(*,*) "Visualization of gradients"
1602 1 pfleura2
        Idx=min(12,NCoord)
1603 1 pfleura2
        DO I=1,NGeomF
1604 1 pfleura2
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1605 1 pfleura2
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1606 1 pfleura2
           WRITe(*,*)
1607 1 pfleura2
        END DO
1608 1 pfleura2
        WRITE(*,*) "END of gradients"
1609 1 pfleura2
     END IF
1610 1 pfleura2
1611 1 pfleura2
     Fini=.TRUE.
1612 1 pfleura2
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1613 1 pfleura2
        IGeom=1
1614 1 pfleura2
        if (debug) WRITE(*,*) '**************************************'
1615 1 pfleura2
        if (debug) WRITE(*,*) 'Optimizing reactant'
1616 1 pfleura2
        if (debug) WRITE(*,*) '**************************************'
1617 1 pfleura2
        SELECT CASE (COORD)
1618 1 pfleura2
        CASE ('ZMAT','MIXED','BAKER')
1619 1 pfleura2
           SELECT CASE (Step_method)
1620 1 pfleura2
           CASE ('RFO')
1621 1 pfleura2
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
1622 1 pfleura2
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1623 1 pfleura2
              Print *, GradTmp
1624 1 pfleura2
           CASE ('GDIIS')
1625 1 pfleura2
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1626 1 pfleura2
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1627 1 pfleura2
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1628 1 pfleura2
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1629 1 pfleura2
              Print *, GradTmp
1630 1 pfleura2
           CASE ('GEDIIS')
1631 1 pfleura2
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1632 1 pfleura2
              ! Energies are updated in EgradPath.f90
1633 1 pfleura2
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1634 1 pfleura2
                   NCoord,allocation_flag,ZeroVec)
1635 1 pfleura2
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1636 1 pfleura2
              !GradTmp = GeomTmp - IntCoordF(1,:)
1637 1 pfleura2
              !Print *, 'Old Geometry:'
1638 1 pfleura2
              !Print *, IntCoordF(1,:)
1639 1 pfleura2
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1640 1 pfleura2
              Print *, GradTmp
1641 1 pfleura2
              !Print *, 'GeomTmp:'
1642 1 pfleura2
              !Print *, GeomTmp
1643 1 pfleura2
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1644 1 pfleura2
           CASE DEFAULT
1645 1 pfleura2
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1646 1 pfleura2
              STOP
1647 1 pfleura2
           END SELECT
1648 1 pfleura2
        CASE ('HYBRID','CART')
1649 1 pfleura2
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1650 1 pfleura2
        CASE DEFAULT
1651 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1652 1 pfleura2
           STOP
1653 1 pfleura2
        END SELECT
1654 1 pfleura2
1655 1 pfleura2
        IF (debug) THEN
1656 1 pfleura2
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1657 1 pfleura2
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1658 1 pfleura2
        END IF
1659 1 pfleura2
1660 1 pfleura2
        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
1661 1 pfleura2
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
1662 1 pfleura2
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
1663 1 pfleura2
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1664 1 pfleura2
        OptReac=(NormStep.GT.SThresh)
1665 1 pfleura2
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
1666 1 pfleura2
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
1667 1 pfleura2
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1668 1 pfleura2
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
1669 1 pfleura2
        !Print *, 'Grad(1,:):'
1670 1 pfleura2
        !Print *, Grad(1,:)
1671 1 pfleura2
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac
1672 1 pfleura2
1673 1 pfleura2
        GradTmp=GradTmp*FactStep
1674 1 pfleura2
1675 1 pfleura2
        ! we take care that frozen atoms do not move.
1676 1 pfleura2
        IF (NFroz.NE.0) THEN
1677 1 pfleura2
           SELECT CASE (COORD)
1678 1 pfleura2
           CASE ('ZMAT','MIXED')
1679 1 pfleura2
              GradTmp(1:IntFroz)=0.d0
1680 1 pfleura2
           CASE ('CART','HYBRID')
1681 1 pfleura2
              DO I=1,Nat
1682 1 pfleura2
                 IF (FrozAtoms(I)) THEN
1683 1 pfleura2
                    Iat=Order(I)
1684 1 pfleura2
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1685 1 pfleura2
                 END IF
1686 1 pfleura2
              END DO
1687 1 pfleura2
           CASE('BAKER')
1688 1 pfleura2
              GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..)
1689 1 pfleura2
           CASE DEFAULT
1690 1 pfleura2
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP'
1691 1 pfleura2
              STOP
1692 1 pfleura2
           END SELECT
1693 1 pfleura2
        END IF ! matches IF (NFroz.NE.0) THEN
1694 1 pfleura2
1695 1 pfleura2
        IF (debug) THEN
1696 1 pfleura2
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1697 1 pfleura2
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1698 1 pfleura2
        END IF
1699 1 pfleura2
1700 1 pfleura2
        IF (debug) THEN
1701 1 pfleura2
           !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1'
1702 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1703 1 pfleura2
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1704 1 pfleura2
           ELSE
1705 1 pfleura2
              DO Iat=1,Nat
1706 1 pfleura2
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1707 1 pfleura2
                      XyzGeomF(IGeom,1:3,Iat)
1708 1 pfleura2
              END DO
1709 1 pfleura2
           END IF
1710 1 pfleura2
        END IF ! matches IF (debug) THEN
1711 1 pfleura2
1712 1 pfleura2
        SELECT CASE (COORD)
1713 1 pfleura2
        CASE ('ZMAT','MIXED','BAKER')
1714 1 pfleura2
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1715 1 pfleura2
        CASE ('HYBRID','CART')
1716 1 pfleura2
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1717 1 pfleura2
        CASE DEFAULT
1718 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1719 1 pfleura2
           STOP
1720 1 pfleura2
        END SELECT
1721 1 pfleura2
1722 1 pfleura2
        IF (debug) THEN
1723 1 pfleura2
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1724 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1725 1 pfleura2
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1726 1 pfleura2
           ELSE
1727 1 pfleura2
              DO Iat=1,Nat
1728 1 pfleura2
1729 1 pfleura2
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1730 1 pfleura2
                 ! It was with indzmat. IndZmat cannot appear here...
1731 1 pfleura2
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1732 1 pfleura2
                      XyzGeomF(IGeom,1:3,Iat)
1733 1 pfleura2
              END DO
1734 1 pfleura2
           END IF
1735 1 pfleura2
        END IF ! matches IF (debug) THEN
1736 1 pfleura2
1737 1 pfleura2
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1738 1 pfleura2
     ELSE ! Optreac
1739 1 pfleura2
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1740 1 pfleura2
     END IF ! matches IF (OptReac) THEN
1741 1 pfleura2
1742 1 pfleura2
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1743 1 pfleura2
        IGeom=NGeomF
1744 1 pfleura2
        if (debug) WRITE(*,*) '******************'
1745 1 pfleura2
        if (debug) WRITE(*,*) 'Optimizing product'
1746 1 pfleura2
        if (debug) WRITE(*,*) '******************'
1747 1 pfleura2
        SELECT CASE (COORD)
1748 1 pfleura2
        CASE ('ZMAT','MIXED','BAKER')
1749 1 pfleura2
           Print *, 'Optimizing product'
1750 1 pfleura2
           SELECT CASE (Step_method)
1751 1 pfleura2
           CASE ('RFO')
1752 1 pfleura2
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1753 1 pfleura2
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1754 1 pfleura2
              Print *, GradTmp
1755 1 pfleura2
           CASE ('GDIIS')
1756 1 pfleura2
              ! GradTmp becomes "step" below and has INTENT(OUT):
1757 1 pfleura2
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1758 1 pfleura2
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1759 1 pfleura2
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1760 1 pfleura2
              Print *, GradTmp
1761 1 pfleura2
           CASE ('GEDIIS')
1762 1 pfleura2
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1763 1 pfleura2
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1764 1 pfleura2
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1765 1 pfleura2
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1766 1 pfleura2
              Print *, GradTmp
1767 1 pfleura2
           CASE DEFAULT
1768 1 pfleura2
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1769 1 pfleura2
              STOP
1770 1 pfleura2
           END SELECT
1771 1 pfleura2
        CASE ('HYBRID','CART')
1772 1 pfleura2
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1773 1 pfleura2
        CASE DEFAULT
1774 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1775 1 pfleura2
           STOP
1776 1 pfleura2
        END SELECT
1777 1 pfleura2
1778 1 pfleura2
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1779 1 pfleura2
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
1780 1 pfleura2
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1781 1 pfleura2
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1782 1 pfleura2
        OptProd=.NOT.(NormStep.LE.SThresh)
1783 1 pfleura2
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1784 1 pfleura2
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1785 1 pfleura2
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1786 1 pfleura2
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1787 1 pfleura2
1788 1 pfleura2
        GradTmp=GradTmp*FactStep
1789 1 pfleura2
1790 1 pfleura2
        ! we take care that frozen atoms do not move
1791 1 pfleura2
        IF (NFroz.NE.0) THEN
1792 1 pfleura2
           SELECT CASE (COORD)
1793 1 pfleura2
           CASE ('ZMAT','MIXED','BAKER')
1794 1 pfleura2
              GradTmp(1:IntFroz)=0.d0
1795 1 pfleura2
           CASE ('CART','HYBRID')
1796 1 pfleura2
              DO I=1,Nat
1797 1 pfleura2
                 IF (FrozAtoms(I)) THEN
1798 1 pfleura2
                    Iat=Order(I)
1799 1 pfleura2
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1800 1 pfleura2
                 END IF
1801 1 pfleura2
              END DO
1802 1 pfleura2
           CASE DEFAULT
1803 1 pfleura2
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1804 1 pfleura2
              STOP
1805 1 pfleura2
           END SELECT
1806 1 pfleura2
        END IF ! matches IF (NFroz.NE.0) THEN
1807 1 pfleura2
1808 1 pfleura2
        IF (debug) THEN
1809 1 pfleura2
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
1810 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1811 1 pfleura2
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1812 1 pfleura2
           ELSE
1813 1 pfleura2
              DO Iat=1,Nat
1814 1 pfleura2
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1815 1 pfleura2
                      XyzGeomF(IGeom,1:3,Iat)
1816 1 pfleura2
              END DO
1817 1 pfleura2
           END IF
1818 1 pfleura2
        END IF
1819 1 pfleura2
1820 1 pfleura2
        SELECT CASE (COORD)
1821 1 pfleura2
        CASE ('ZMAT','MIXED','BAKER')
1822 1 pfleura2
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1823 1 pfleura2
        CASE ('HYBRID','CART')
1824 1 pfleura2
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1825 1 pfleura2
        CASE DEFAULT
1826 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1827 1 pfleura2
           STOP
1828 1 pfleura2
        END SELECT
1829 1 pfleura2
1830 1 pfleura2
        IF (debug) THEN
1831 1 pfleura2
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1832 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1833 1 pfleura2
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1834 1 pfleura2
           ELSE
1835 1 pfleura2
              DO Iat=1,Nat
1836 1 pfleura2
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1837 1 pfleura2
                      XyzGeomF(IGeom,1:3,Iat)
1838 1 pfleura2
              END DO
1839 1 pfleura2
           END IF
1840 1 pfleura2
        END IF
1841 1 pfleura2
     ELSE ! Optprod
1842 1 pfleura2
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1843 1 pfleura2
     END IF !matches IF (OptProd) THEN
1844 1 pfleura2
1845 1 pfleura2
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1846 1 pfleura2
1847 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1848 1 pfleura2
     !
1849 1 pfleura2
     !  Optimizing other geometries
1850 1 pfleura2
     !
1851 1 pfleura2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1852 1 pfleura2
1853 1 pfleura2
1854 1 pfleura2
1855 1 pfleura2
     DO IGeom=2,NGeomF-1 ! matches at L1556
1856 1 pfleura2
        if (debug) WRITE(*,*) '****************************'
1857 1 pfleura2
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1858 1 pfleura2
        if (debug) WRITE(*,*) '****************************'
1859 1 pfleura2
1860 1 pfleura2
        SELECT CASE (COORD)
1861 1 pfleura2
        CASE ('ZMAT','BAKER','MIXED')
1862 1 pfleura2
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
1863 1 pfleura2
           ! PFL 6 Apr 2008 ->
1864 1 pfleura2
           ! Special case, if FTan=0. we do not consider Tangent at all.
1865 1 pfleura2
           ! To DO: add the full treatment of FTan
1866 1 pfleura2
           if (FTan==0.) GradTmp2=ZeroVec
1867 1 pfleura2
           ! <- PFL 6 Apr 2008
1868 1 pfleura2
           if (debug) THEN
1869 1 pfleura2
              Print *, 'L1500, IntTangent(',IGeom,',:)='
1870 1 pfleura2
              Print *, IntTangent(IGeom,:)
1871 1 pfleura2
           END IF
1872 1 pfleura2
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
1873 1 pfleura2
           SELECT CASE (Step_method)
1874 1 pfleura2
           CASE ('RFO')
1875 1 pfleura2
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1876 1 pfleura2
           CASE ('GDIIS')
1877 1 pfleura2
              !The following has no effect at IOpt==1
1878 1 pfleura2
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
1879 1 pfleura2
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&
1880 1 pfleura2
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1881 1 pfleura2
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
1882 1 pfleura2
              Print *, GradTmp
1883 1 pfleura2
           CASE ('GEDIIS')
1884 1 pfleura2
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1885 1 pfleura2
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
1886 1 pfleura2
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1887 1 pfleura2
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
1888 1 pfleura2
              Print *, GradTmp
1889 1 pfleura2
           CASE DEFAULT
1890 1 pfleura2
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
1891 1 pfleura2
              STOP
1892 1 pfleura2
           END SELECT
1893 1 pfleura2
        CASE ('HYBRID','CART')
1894 1 pfleura2
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
1895 1 pfleura2
           ! so we change it
1896 1 pfleura2
           DO I=1,Nat
1897 1 pfleura2
              DO J=1,3
1898 1 pfleura2
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
1899 1 pfleura2
              END DO
1900 1 pfleura2
           END DO
1901 1 pfleura2
           !           GradTmp2=XyzTangent(IGeom,:)
1902 1 pfleura2
           ! PFL 6 Apr 2008 ->
1903 1 pfleura2
           ! Special case, if FTan=0. we do not consider Tangent at all.
1904 1 pfleura2
           ! To DO: add the full treatment of FTan
1905 1 pfleura2
           if (FTan==0.) GradTmp2=ZeroVec
1906 1 pfleura2
           ! <- PFL 6 Apr 2008
1907 1 pfleura2
1908 1 pfleura2
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1909 1 pfleura2
        CASE DEFAULT
1910 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
1911 1 pfleura2
           STOP
1912 1 pfleura2
        END SELECT
1913 1 pfleura2
1914 1 pfleura2
        ! we take care that frozen atoms do not move
1915 1 pfleura2
        IF (NFroz.NE.0) THEN
1916 1 pfleura2
           SELECT CASE (COORD)
1917 1 pfleura2
           CASE ('ZMAT','MIXED','BAKER')
1918 1 pfleura2
              IF (debug) THEN
1919 1 pfleura2
                 WRITE(*,*) 'Step computed. Coord=',Coord
1920 1 pfleura2
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
1921 1 pfleura2
              END IF
1922 1 pfleura2
              GradTmp(1:IntFroz)=0.d0
1923 1 pfleura2
           CASE ('CART','HYBRID')
1924 1 pfleura2
              if (debug) THEN
1925 1 pfleura2
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
1926 1 pfleura2
                 DO I=1,Nat
1927 1 pfleura2
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
1928 1 pfleura2
                 END DO
1929 1 pfleura2
              END IF
1930 1 pfleura2
              DO I=1,Nat
1931 1 pfleura2
                 IF (FrozAtoms(I)) THEN
1932 1 pfleura2
                    if (debug) THEN
1933 1 pfleura2
                       write(*,*) 'Step Computed. Coord=',Coord
1934 1 pfleura2
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
1935 1 pfleura2
                    END IF
1936 1 pfleura2
                    Iat=Order(I)
1937 1 pfleura2
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1938 1 pfleura2
                 END IF
1939 1 pfleura2
              END DO
1940 1 pfleura2
1941 1 pfleura2
                 if (debug) THEN
1942 1 pfleura2
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
1943 1 pfleura2
                    DO I=1,Nat
1944 1 pfleura2
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
1945 1 pfleura2
                    END DO
1946 1 pfleura2
                 END IF
1947 1 pfleura2
1948 1 pfleura2
           CASE DEFAULT
1949 1 pfleura2
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
1950 1 pfleura2
              STOP
1951 1 pfleura2
           END SELECT
1952 1 pfleura2
        END IF !matches IF (NFroz.NE.0) THEN
1953 1 pfleura2
1954 1 pfleura2
        IF (debug) THEN
1955 1 pfleura2
           SELECT CASE (COORD)
1956 1 pfleura2
           CASE ('ZMAT','MIXED','BAKER')
1957 1 pfleura2
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
1958 1 pfleura2
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1959 1 pfleura2
           CASE ('CART','HYBRID')
1960 1 pfleura2
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
1961 1 pfleura2
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1962 1 pfleura2
           CASE DEFAULT
1963 1 pfleura2
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
1964 1 pfleura2
              STOP
1965 1 pfleura2
           END SELECT
1966 1 pfleura2
        END IF ! matches if (debug) THEN
1967 1 pfleura2
1968 1 pfleura2
        ! we check that the maximum displacement is not greater than Smax
1969 1 pfleura2
        If (debug) WRITE(*,*) "Fini before test:",Fini
1970 1 pfleura2
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))
1971 1 pfleura2
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
1972 1 pfleura2
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1973 1 pfleura2
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
1974 1 pfleura2
1975 1 pfleura2
        GradTmp=GraDTmp*FactStep
1976 1 pfleura2
1977 1 pfleura2
        if (debug) WRITE(*,*) "DBG Main, L1475, FactStep=",FactStep
1978 1 pfleura2
        if (debug.AND.(COORD.EQ.'ZMAT')) WRITE(*,'(A,12(1X,F10.6))') 'Scaled Step:',GradTmp(1:min(12,NFree))
1979 1 pfleura2
1980 1 pfleura2
        IF (debug) THEN
1981 1 pfleura2
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
1982 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1983 1 pfleura2
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1984 1 pfleura2
           ELSE
1985 1 pfleura2
              DO Iat=1,Nat
1986 1 pfleura2
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1987 1 pfleura2
                      XyzGeomF(IGeom,1:3,Iat)
1988 1 pfleura2
              END DO
1989 1 pfleura2
           END IF
1990 1 pfleura2
        END IF ! matches if (debug) THEN
1991 1 pfleura2
1992 1 pfleura2
        SELECT CASE (COORD)
1993 1 pfleura2
        CASE ('ZMAT')
1994 1 pfleura2
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1995 1 pfleura2
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1996 1 pfleura2
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1997 1 pfleura2
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1998 1 pfleura2
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1999 1 pfleura2
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
2000 1 pfleura2
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2001 1 pfleura2
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
2002 1 pfleura2
           Idx=4
2003 1 pfleura2
           DO Iat=4,Nat
2004 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2005 1 pfleura2
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2006 1 pfleura2
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2007 1 pfleura2
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2008 1 pfleura2
              Idx=Idx+3
2009 1 pfleura2
           END DO
2010 1 pfleura2
2011 1 pfleura2
        CASE ('BAKER')
2012 1 pfleura2
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2013 1 pfleura2
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2014 1 pfleura2
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2015 1 pfleura2
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
2016 1 pfleura2
2017 1 pfleura2
        CASE ('MIXED')
2018 1 pfleura2
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2019 1 pfleura2
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2020 1 pfleura2
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2021 1 pfleura2
           DO Iat=1,NCart
2022 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
2023 1 pfleura2
           END DO
2024 1 pfleura2
           Idx=3*NCart+1
2025 1 pfleura2
           IF (NCart.EQ.1) THEN
2026 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2027 1 pfleura2
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
2028 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2029 1 pfleura2
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
2030 1 pfleura2
2031 1 pfleura2
              Idx=7
2032 1 pfleura2
           END IF
2033 1 pfleura2
           IF (NCart.EQ.2) THEN
2034 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2035 1 pfleura2
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
2036 1 pfleura2
              Idx=9
2037 1 pfleura2
           END IF
2038 1 pfleura2
2039 1 pfleura2
2040 1 pfleura2
           DO Iat=max(NCart+1,4),Nat
2041 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2042 1 pfleura2
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2043 1 pfleura2
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2044 1 pfleura2
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2045 1 pfleura2
              Idx=Idx+3
2046 1 pfleura2
           END DO
2047 1 pfleura2
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2048 1 pfleura2
2049 1 pfleura2
        CASE ('HYBRID','CART')
2050 1 pfleura2
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
2051 1 pfleura2
        CASE DEFAULT
2052 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
2053 1 pfleura2
           STOP
2054 1 pfleura2
        END SELECT
2055 1 pfleura2
2056 1 pfleura2
        IF (debug) THEN
2057 1 pfleura2
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
2058 1 pfleura2
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2059 1 pfleura2
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
2060 1 pfleura2
           ELSE
2061 1 pfleura2
              DO Iat=1,Nat
2062 1 pfleura2
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
2063 1 pfleura2
                      XyzGeomF(IGeom,1:3,Iat)
2064 1 pfleura2
              END DO
2065 1 pfleura2
           END IF
2066 1 pfleura2
        END IF ! matches if (debug) THEN
2067 1 pfleura2
2068 1 pfleura2
        ! We project out the tangential part of the gradient to know if we are done
2069 1 pfleura2
        GradTmp=Grad(IGeom,:)
2070 1 pfleura2
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
2071 1 pfleura2
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
2072 1 pfleura2
        SELECT CASE (COORD)
2073 1 pfleura2
        CASE ('ZMAT','MIXED','BAKER')
2074 1 pfleura2
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2075 1 pfleura2
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
2076 1 pfleura2
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2077 1 pfleura2
           Ca=S/(sqrt(Norm)*NormGrad)
2078 1 pfleura2
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2079 1 pfleura2
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2080 1 pfleura2
        CASE ('CART','HYBRID')
2081 1 pfleura2
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2082 1 pfleura2
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
2083 1 pfleura2
           Ca=S/(sqrt(Norm)*NormGrad)
2084 1 pfleura2
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2085 1 pfleura2
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2086 1 pfleura2
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2087 1 pfleura2
        CASE DEFAULT
2088 1 pfleura2
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
2089 1 pfleura2
           STOP
2090 1 pfleura2
        END SELECT
2091 1 pfleura2
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
2092 1 pfleura2
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
2093 1 pfleura2
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
2094 1 pfleura2
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
2095 1 pfleura2
2096 1 pfleura2
     END DO ! matches DO IGeom=2,NGeomF-1
2097 1 pfleura2
     ! we save the old gradients
2098 1 pfleura2
     GradOld=Grad
2099 1 pfleura2
     EPathOld=Energies
2100 1 pfleura2
2101 7 pfleura2
2102 7 pfleura2
2103 1 pfleura2
     ! Shall we re-parameterize the path ?
2104 1 pfleura2
     ! PFL 2007/Apr/10 ->
2105 1 pfleura2
     ! We call PathCreate in all cases, it will take care of the
2106 1 pfleura2
     ! reparameterization, as well as calculating the tangents.
2107 1 pfleura2
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
2108 1 pfleura2
     XyzGeomI=XyzGeomF
2109 1 pfleura2
     IntCoordI=IntCoordF
2110 1 pfleura2
2111 1 pfleura2
     Call PathCreate(IOpt)
2112 1 pfleura2
     !     END IF
2113 1 pfleura2
     ! <- PFL 2007/Apr/10
2114 1 pfleura2
2115 1 pfleura2
     if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)')  'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini
2116 1 pfleura2
     IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN
2117 1 pfleura2
2118 1 pfleura2
        ! We have the new path, we calculate its energy and gradient
2119 1 pfleura2
2120 2 pfleura2
        Call EgradPath(IOpt)
2121 1 pfleura2
        !IF(IOPT .EQ. 10) Then
2122 1 pfleura2
        !         Print *, 'Stopping at my checkpoint.'
2123 1 pfleura2
        !           STOP !This is my temporary checkpoint.
2124 1 pfleura2
        !ENDIF
2125 1 pfleura2
2126 1 pfleura2
        ! PFL 31 Mar 2008 ->
2127 1 pfleura2
        ! Taken from v3.99 but we added a flag...
2128 1 pfleura2
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
2129 1 pfleura2
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1
2130 1 pfleura2
        ! else it is decreased by 0.8
2131 1 pfleura2
2132 1 pfleura2
        if (DynMaxStep) THEN
2133 1 pfleura2
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
2134 1 pfleura2
           if (OptReac) THEN
2135 1 pfleura2
              If (Energies(1)<EPathOld(1)) Then
2136 1 pfleura2
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
2137 1 pfleura2
              ELSE
2138 1 pfleura2
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
2139 1 pfleura2
              END IF
2140 1 pfleura2
           END IF
2141 1 pfleura2
2142 1 pfleura2
           if (OptProd) THEN
2143 1 pfleura2
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
2144 1 pfleura2
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
2145 1 pfleura2
              ELSE
2146 1 pfleura2
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
2147 1 pfleura2
              END IF
2148 1 pfleura2
           END IF
2149 1 pfleura2
2150 1 pfleura2
           DO IGeom=2,NGeomF-1
2151 1 pfleura2
              If (Energies(IGeom)<EPathOld(IGeom)) Then
2152 1 pfleura2
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
2153 1 pfleura2
              ELSE
2154 1 pfleura2
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
2155 1 pfleura2
              END IF
2156 1 pfleura2
           END DO
2157 1 pfleura2
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2158 1 pfleura2
        END IF ! test on DynMaxStep
2159 1 pfleura2
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2160 1 pfleura2
        ! <- PFL 31 MAr 2008
2161 1 pfleura2
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
2162 1 pfleura2
        !  It should get updated for other cases too.
2163 1 pfleura2
2164 1 pfleura2
        DO IGeom=1,NGeomF
2165 1 pfleura2
           SELECT CASE (COORD)
2166 1 pfleura2
           CASE('ZMAT')
2167 1 pfleura2
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2168 1 pfleura2
              GeomTmp=IntCoordF(IGeom,:)
2169 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2170 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2171 1 pfleura2
                   OrderInv(indzmat(2,2)),Geomtmp(1)
2172 1 pfleura2
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2173 1 pfleura2
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
2174 1 pfleura2
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
2175 1 pfleura2
              Idx=4
2176 1 pfleura2
              DO Iat=4,Nat
2177 1 pfleura2
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2178 1 pfleura2
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
2179 1 pfleura2
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
2180 1 pfleura2
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
2181 1 pfleura2
                 Idx=Idx+3
2182 1 pfleura2
              END DO
2183 1 pfleura2
           CASE('BAKER')
2184 1 pfleura2
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2185 1 pfleura2
              !GeomTmp=IntCoordF(IGeom,:)
2186 1 pfleura2
           CASE('CART','HYBRID','MIXED')
2187 1 pfleura2
              WRITE(IOOUT,'(1X,I5)') Nat
2188 1 pfleura2
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2189 1 pfleura2
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2190 1 pfleura2
              DO I=1,Nat
2191 1 pfleura2
                 Iat=I
2192 1 pfleura2
                 If (renum) Iat=OrderInv(I)
2193 1 pfleura2
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
2194 1 pfleura2
              END DO
2195 1 pfleura2
           CASE DEFAULT
2196 1 pfleura2
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
2197 1 pfleura2
              STOP
2198 1 pfleura2
           END SELECT
2199 1 pfleura2
        END DO ! matches DO IGeom=1,NGeomF
2200 1 pfleura2
2201 1 pfleura2
        Call Write_path(IOpt)
2202 7 pfleura2
! Shall we analyze the geometries ?
2203 7 pfleura2
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
2204 1 pfleura2
2205 1 pfleura2
        ! We update the Hessian matrices
2206 1 pfleura2
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
2207 1 pfleura2
        ! First using the displacement from the old path
2208 1 pfleura2
        IGeom0=2
2209 1 pfleura2
        IGeomF=NGeomF-1
2210 1 pfleura2
        IF (OptReac) IGeom0=1
2211 1 pfleura2
        IF (OptProd) IGeomF=NGeomF
2212 1 pfleura2
2213 1 pfleura2
        IF (FCalcHess) THEN
2214 1 pfleura2
           DO IGeom=IGeom0,IGeomF
2215 1 pfleura2
              SELECT CASE (COORD)
2216 1 pfleura2
              CASE ('ZMAT','MIXED','BAKER')
2217 1 pfleura2
                 GeomTmp2=IntCoordF(IGeom,:)
2218 1 pfleura2
              CASE ('CART','HYBRID')
2219 1 pfleura2
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2220 1 pfleura2
              CASE DEFAULT
2221 1 pfleura2
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2222 1 pfleura2
                 STOP
2223 1 pfleura2
              END SELECT
2224 1 pfleura2
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
2225 1 pfleura2
           END DO
2226 1 pfleura2
        ELSE
2227 1 pfleura2
           DO IGeom=IGeom0,IGeomF
2228 1 pfleura2
              GeomTmp=GeomOld_all(IGeom,:)
2229 1 pfleura2
              SELECT CASE (COORD)
2230 1 pfleura2
              CASE ('ZMAT','MIXED','BAKER')
2231 1 pfleura2
                 GeomTmp2=IntCoordF(IGeom,:)
2232 1 pfleura2
              CASE ('CART','HYBRID')
2233 1 pfleura2
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2234 1 pfleura2
              CASE DEFAULT
2235 1 pfleura2
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2236 1 pfleura2
                 STOP
2237 1 pfleura2
              END SELECT
2238 1 pfleura2
              GeomTmp=GeomTmp2-GeomTmp
2239 1 pfleura2
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
2240 1 pfleura2
              HessTmp=Hess(IGeom,:,:)
2241 1 pfleura2
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2242 1 pfleura2
              Hess(IGeom,:,:)=HessTmp
2243 1 pfleura2
           END DO ! matches DO IGeom=IGeom0,IGeomF
2244 1 pfleura2
2245 1 pfleura2
           ! Second using the neighbour points
2246 1 pfleura2
           IF (HupNeighbour) THEN
2247 1 pfleura2
              ! Only one point for end points.
2248 1 pfleura2
              IF (OptReac) THEN
2249 1 pfleura2
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2250 1 pfleura2
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
2251 1 pfleura2
                 ELSE
2252 1 pfleura2
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
2253 1 pfleura2
                 END IF
2254 1 pfleura2
                 GradTmp=Grad(1,:)-Grad(2,:)
2255 1 pfleura2
                 HessTmp=Hess(1,:,:)
2256 1 pfleura2
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2257 1 pfleura2
                 Hess(1,:,:)=HessTmp
2258 1 pfleura2
              END IF
2259 1 pfleura2
2260 1 pfleura2
              IF (OptProd) THEN
2261 1 pfleura2
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2262 1 pfleura2
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
2263 1 pfleura2
                 ELSE
2264 1 pfleura2
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
2265 1 pfleura2
                 END IF
2266 1 pfleura2
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
2267 1 pfleura2
                 HessTmp=Hess(NGeomF,:,:)
2268 1 pfleura2
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2269 1 pfleura2
                 Hess(NGeomF,:,:)=HessTmp
2270 1 pfleura2
              END IF
2271 1 pfleura2
2272 1 pfleura2
              ! Two points for the rest of the path
2273 1 pfleura2
              DO IGeom=2,NGeomF-1
2274 1 pfleura2
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
2275 1 pfleura2
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
2276 1 pfleura2
                 ELSE
2277 1 pfleura2
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
2278 1 pfleura2
                 END IF
2279 1 pfleura2
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
2280 1 pfleura2
                 HessTmp=Hess(IGeom,:,:)
2281 1 pfleura2
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2282 1 pfleura2
2283 1 pfleura2
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2284 1 pfleura2
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
2285 1 pfleura2
                 ELSE
2286 1 pfleura2
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
2287 1 pfleura2
                 END IF
2288 1 pfleura2
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
2289 1 pfleura2
2290 1 pfleura2
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2291 1 pfleura2
                 Hess(IGeom,:,:)=HessTmp
2292 1 pfleura2
              END DO ! matches DO IGeom=2,NGeomF-1
2293 1 pfleura2
           END IF !matches IF (HupNeighbour) THEN
2294 1 pfleura2
        END IF ! matches IF (FCalcHess)
2295 1 pfleura2
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
2296 1 pfleura2
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
2297 1 pfleura2
2298 2 pfleura2
!   IF (PROG=="GAUSSIAN") THEN
2299 2 pfleura2
!      DEALLOCATE(Gauss_paste)
2300 2 pfleura2
!      DEALLOCATE(Gauss_root)
2301 2 pfleura2
!      DEALLOCATE(Gauss_comment)
2302 2 pfleura2
!      DEALLOCATE(Gauss_end)
2303 2 pfleura2
!   END IF
2304 2 pfleura2
2305 1 pfleura2
  DEALLOCATE(XyzGeomF, IntCoordF)
2306 1 pfleura2
  DEALLOCATE(XyzGeomI, IntCoordI)
2307 1 pfleura2
  DEALLOCATE(XyzTangent,IntTangent)
2308 1 pfleura2
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2309 1 pfleura2
  DEALLOCATE(GeomOld_all)
2310 1 pfleura2
2311 1 pfleura2
  if (PROG=="GAUSSIAN") THEN
2312 1 pfleura2
     ! We de-allocate the Gauss_XX lists
2313 1 pfleura2
     ! transverse the list and deallocate each node
2314 1 pfleura2
     previous => Gauss_root ! make current point to head of list
2315 1 pfleura2
     DO
2316 1 pfleura2
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2317 1 pfleura2
        current => previous%next ! make list point to next node of head
2318 1 pfleura2
        DEALLOCATE(previous) ! deallocate current head node
2319 1 pfleura2
        previous => current  ! make current point to new head
2320 1 pfleura2
     END DO
2321 1 pfleura2
2322 1 pfleura2
     previous => Gauss_comment ! make current point to head of list
2323 1 pfleura2
     DO
2324 1 pfleura2
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2325 1 pfleura2
        current => previous%next ! make list point to next node of head
2326 1 pfleura2
        DEALLOCATE(previous) ! deallocate current head node
2327 1 pfleura2
        previous => current  ! make current point to new head
2328 1 pfleura2
     END DO
2329 1 pfleura2
2330 1 pfleura2
     previous => Gauss_end ! make current point to head of list
2331 1 pfleura2
     DO
2332 1 pfleura2
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2333 1 pfleura2
        current => previous%next ! make list point to next node of head
2334 1 pfleura2
        DEALLOCATE(previous) ! deallocate current head node
2335 1 pfleura2
        previous => current  ! make current point to new head
2336 1 pfleura2
     END DO
2337 1 pfleura2
2338 2 pfleura2
2339 1 pfleura2
  END IF
2340 1 pfleura2
2341 1 pfleura2
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2342 1 pfleura2
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2343 1 pfleura2
  IF (COORD.EQ.'BAKER') THEN
2344 1 pfleura2
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2345 1 pfleura2
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2346 1 pfleura2
  END IF
2347 1 pfleura2
2348 1 pfleura2
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2349 1 pfleura2
2350 7 pfleura2
  Close(IOIN)
2351 7 pfleura2
  Close(IOOUT)
2352 7 pfleura2
  IF (AnaGeom) Close(IODAT)
2353 8 pfleura2
  IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot)
2354 7 pfleura2
2355 1 pfleura2
END PROGRAM PathOpt