Statistiques
| Révision :

root / src / Path.f90 @ 8

Historique | Voir | Annoter | Télécharger (83,71 ko)

1
! This programs generates and optimizes a reaction path
2
! between at least two endpoints.
3
!
4
! It uses NAMELIST for input, see below !
5
!
6
! v3.0
7
! First version entirely in Fortran.
8
! It uses routines from Cart (path creation and parameterisation) 
9
! and from IRC06, especially routines for point optimization.
10
!
11
! TO DO:
12
! 1) Possibility to read and/or calculate some hessian structures
13
! 2) Use of internal coordinate to estimate the hessian for stable
14
! species in a much better way...
15
!
16
!!!!!!!!!!!!!!!
17
! v3.1
18
! The Hessian update includes neighbour points.
19
!
20
! v3.2
21
! We add the option Hinv that uses the BFGS update of the Hessian inverse
22
!
23
! v3.3
24
! We add the full option for ZMAT coordinates, where all is done using
25
! internal coordinates.
26
! v3.31
27
! The step size is controlled using the step norm and not the maximum step
28
! component.
29
! v3.5
30
! Big modifications of Calc_zmat_const_frag in order to be able
31
! to treat VASP geometries. 
32
! For now, only XYZ geometries are read and written. VASP interface
33
! has to be written.
34
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
35
! v3.6
36
! Some minor modif in Calc_zmat_const_frag to have nicer programming
37
! and a new mode added: Coord=MIXED where part of the system
38
! is extrapolated in cartesian coordinates and the rest of it in zmat.
39
! it might be Useful for VASP, or when lots of atoms are frozen
40
! by default, when Coord=MIXED, all frozen atoms are described in cartesian
41
! coordinates.
42
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
43
! v3.61
44
! We move further and further to VASP program. New variables added to the 'path' namelist:
45
! input: what is the program used for the geometries ?
46
!     Possible values: xyz (or cart) for Xmol format; VASP. Xyz should be used for Gaussian
47
! prog: what is the program used for the calculations ?
48
!     Possible values for now: gaussian, vasp.
49
! For now (02/2007), when prog='vasp' is used, only the initial path is created.
50
! That means that we have only added input/output subroutines :-)
51
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
52
! v3.7
53
! We try to improve the optimization steps. For that, we first modify many routines
54
! so that the new step is calculated in the 'free' degrees of freedom, as it is done
55
! in IRC07/ADF for example. That will allow us to impose constraints in an easier way.
56
!
57
! v3.701
58
! Add a new option for the initial extrapolation of the Hessian
59
!
60
!v3.71
61
! The optimization step has been improved... but the vfree part has not yet been included.
62
! This has still to be done :-)
63
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
64
! v3.8
65
! I have included part of the vfree procedure: for now freemv returns the Id
66
! matrix, that is used to project out the tangent vector.
67
! This improves the optimization a lot, but I will have to implement
68
! a real freemw subroutine to have constraints there !
69
! v3.81 Technical version
70
! I have added a new program 'Ext'. When this one is used, Path creates a Ext.in file,
71
! calls Ext.exe file to get a Ext.out file that contains the Energy and the gradient.
72
! Format:
73
! Ext.in : Xmol format. The second line (comment in Xmol) contains the COORD name.
74
! Ext.out: Energy on firt line (1X,D25.15), then gradient in CARTESIAN coordinates
75
!  (3(1X,D25.15)) format. Natoms lines. 
76
! TO DO:  Gradient could be given in COORD format, ie cartesian for COORD=CART, XYZ or HYBRID
77
! ZMAT for CORD=ZMAT (in that case, 6 zeros are there at the begining !)
78
! For MIXED, it is equivalent to CART :-)
79
! v3.811
80
! I have added a prog="TEST" option, where egradient calls egrad_test subroutine.
81
! It then calculates the energy and gradient as it wants and it returns the cartesian 
82
! gradient. The purpose is to have an analytical approximation of the PES of the system
83
! under study to gain some time in testing the program. This is not documented :-)
84
!
85
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
86
!  v3.9
87
! In order to be fully interfaced with VASP, I have to change the architecture
88
! of the program. In particular, with VASP, one can use a Para8+NEB routine
89
! of VASP to calculate energy+forces for all points of the path.
90
! v3.91
91
! One small modification: in the cart list, one can freeze a whole sequence by
92
! entering a negative number for the last atom of the sequence
93
! for example : cart=1 -160 165 -170 will define atoms from 1 to 160 and from
94
! 165 to 170 (included) as cartesian atoms.  
95
! Of course, Same applies to Frozen
96
! v3.92 / v3.93
97
! Slight modifications. One noticeable: when using Vasp, the program analyses the
98
! displacements to suggest which atoms could be described in CART.
99
! v3.94 
100
! VaspMode=Para now implemented. One needs to be careful in defining 'ProgExe'.
101
! In particular, one must put the mpirun command in ProgExe !!! There is NO test.
102
! v3.95
103
! When using internal coordinates (zmat, hybrid or mixes), Path calculates the number
104
! of fragments for each geometry. The reference one (ie the one used to compute the internal
105
! coordinates) is now the one with the least fragments.
106
! user can also choose one geometry by specifying IGeomRef in the Path namelist.
107
! v3.96
108
! I have added a test for the geometry step: it does not allow valence angle to 
109
! be negative or bigger than 180. I hope that this might help geometry optimization
110
! when dealing with nearly linear molecules.
111
! v3.961
112
! The options cart and frozen are now LOGICAL:
113
! Fcart=T indicates that a &cartlist namelist should be read !
114
! Ffrozen=T indicates that a &frozenlist namelist should be read !
115
! TO DO: Autocart=T indicates that the atoms described in cartesian coordiantes should
116
! be found automatically
117
! v3.97
118
! Baker has been implemented by P. Dayal: we are interpolating Xint for now.
119
! TO do: interpolate Umat ?
120
! -----
121
! v3.971
122
! AutoCart=T is being implemented. Goal: having a 'black box' version for Vasp users.
123
! Vmd: True if user wants to watch the path using VMD. Only used for VASP for now.
124
!------------
125
! v3.98
126
! Calc_zmat has been replaced by Calc_zmat_frag which is more robust.
127
! Still missing: linear molecules and some tests.
128
!
129
! v3.981
130
! * Added a non document flag in .valid file that switch on the numerical
131
! calculation of the hessian instead of the update.
132
! * Added HesUpd variable in the namelist. Default value is MS for path
133
!   optimization because it does not imposes a positive hessian, and BFGS
134
!  for geometry optimization
135
! v 4.0
136
! * This version has been made public within the CARTE project.
137
! * Added:
138
! - Step_method to choose the optimization method: RFO or Diis
139
! - DynMaxStep: if TRUE, the maximum size of a step is updated at each step. 
140
!             If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. 
141
!             It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]
142
!
143
!!!!!!!!!!!!!!!!!!!!!!
144
! v4.1
145
! * Para mode has been partly implemented for Gaussian.
146
!  VaspMode is thus now 'RunMode' and can be Serial or Para
147
! Only Gaussian and VASP are supported for Para mode now.
148
! v4.12 
149
! Some bugs in unconstrained zmat construction are corrected.
150
! v4.13
151
! Prakash has implemented the GEDIIS optimization procedure.
152
! v4.14
153
! There is a small problem for some inputs in VASP, where the
154
! last geometry is not the one given by the user. 
155
! It seems to come from the fact that some people try to freeze
156
! some unfrozen atoms
157
! So some tests have been added to check that frozen atoms do not move.
158
!!!
159
! v4.15
160
! There were some bugs when reading and expanding cartlist and frozenlist
161
! I have changed the way frozen atoms are compared by adding partial alignment.
162
! v4.16
163
! Some bugs were corrected.
164
! v4.17
165
! Added MOPAC as a program.
166
! v4.175
167
! Added Three more parameters for defining the input and output names for
168
! the calculations.
169
! CalcName: Prefix for the files used for the energy and gradient calculations
170
! ISuffix: Suffix for the input file
171
! OSuffix: suffix for the output file.
172
! This maters for Gaussian for example.
173
!
174
! v4.177 (P) 2010 Nov 22
175
! - We add TurboMole as a new external code.
176
! - Subtle change for Input: the default is XYZ whatever the code, except
177
! for VASP. 
178
!
179
! v4.178 (P) 2010 Nov 28
180
! - We add the possibility to call CHAMFRE reactive force field using 
181
! prog=TEST2 or CHAMFRE
182
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183
! OpenPath v1.4  2011 May-June
184
! - Also added two flags: CalcEReac and CalcEProd. If .TRUE. Path will
185
! compute the energy of the reactant
186
! (and/or Product) at the begining of the calculations.
187
! This is useful if you start from non optimized points.
188
!
189
! - Also added more methods for Hessian updates. The following methods are
190
! now availables:
191
!   For inverse Hessian: BFGS, DFP, MS. 
192
!   For Hessian: BFGS, DFP, MS, Bofill
193
! 
194
! -Small change: HesUpd is replaced by Hupdate
195
!
196
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197
! TO DO
198
! - Trying to improve the output of Path. In particular, we now create
199
! a .datl file containing the image index, curvilinear coordinate and energy
200
! in the spirit of the AnaPath tool.
201
! AnaPath is still usefull as it can extract geometrical parameters on top
202
! of the energy
203
! - When using Hessian update, it is more efficient to use Cholesky 
204
! decomposition, and to update this decomposition rather than the Hessian.
205
! See Nocedal&Wright p141.
206
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
207
!
208
! Project Name has changed to Opt'n Path
209
! OpenPaht 2011 May-June -> Opt'n Path v1.4 
210
!
211
! 2013 Jan 18 Opt'n Path v1.47 (going to 1.5 slowly)
212
! We add Siesta (and start thinking of add CP2K)
213
!  as a new "energy engine" code
214
! TO DO to go to 1.5: finish cleaning...
215
!
216
! 2013 Feb Opt'n Path v1.48
217
! We add some keyword for more output for Geometry Optimizations
218
! GeomFile: name of the file to print out the geometries and their energy
219
! as Xyz file. (only format for now)
220
!
221
! 2013 Feb opt'n Path v1.49
222
! We add the possibility to analyze the geometries in terms
223
! of Bonds, Angles, dihedrals
224
! with the option to add CenterOfMass atoms for the analysis (a la Xyz2Path)
225
! This is basicaly the merging of AnaPath into Opt'n Path... at last !
226
! This is done by adding a Flag into the Path namelist:
227
! AnaGeom: logical, if true Path looks for the AnaList namelist
228
! AnaList: list of variables, if present and if AnaGeom=T then Opt'n Path
229
! will read it and print the values of the variable in a Path.dat file
230
! If AnaGeom is T but Analist is not given, then Path.dat just contains the energy. 
231

    
232
PROGRAM PathOpt
233

    
234
  use VarTypes
235
  use Path_module
236
  use Io_module
237

    
238
  IMPLICIT NONE
239

    
240
  CHARACTER(132) :: FileIn, FileOut, Line
241

    
242
  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF
243
  INTEGER(KINT) :: IOpt
244
  INTEGER(KINT) :: Idx
245
  INTEGER(KINT) :: N3at, NFree, Nfr
246

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

    
249
  REAL(KREAL) :: E,FactStep !E=Energy
250
  REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad
251
  REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E
252
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord
253
  REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord)
254
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord
255
  REAL(KREAL), ALLOCATABLE :: HessTmp(:,:)
256
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
257
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
258
  LOGICAL, SAVE :: allocation_flag
259

    
260
  ! Flag to see if frozen atoms are frozen (see below)
261
  LOGICAL :: FChkFrozen
262

    
263
  ! Energies for all points along the path at the previous iteration
264
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF
265
  ! Maximum step allowed for this geometry
266
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF
267

    
268
! these are used to read temporary coordinates
269
  LOGICAL :: FFrozen,FCart
270

    
271
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272
  !
273
  ! To analyse the frozen atoms
274
  !
275
  REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat
276
  REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat
277
  REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist
278
  LOGICAL, ALLOCATABLE :: ListAtFroz(:)
279
  INTEGER(KINT) :: Iat1, Iat2, Iat3
280

    
281
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282

    
283
  LOGICAL :: Debug, Fini,Flag_Opt_Geom
284

    
285
!  INTEGER(KINT), EXTERNAL :: Iargc
286
  INTEGER(KINT), EXTERNAL ::  ConvertNumAt
287
  INTEGER(KINT) :: IoS
288

    
289
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
290
  !
291
  ! For Debugging purposes
292
  !
293
  LOGICAL :: FCalcHess
294
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
295

    
296

    
297
  INTERFACE
298
     function valid(string) result (isValid)
299
       CHARACTER(*), intent(in) :: string
300
       logical                  :: isValid
301
     END function VALID
302

    
303
     SUBROUTINE Read_Geom(Input)
304
       CHARACTER(32), INTENT(IN) :: input
305
     END SUBROUTINE READ_GEOM
306

    
307

    
308
     subroutine hupdate_all (n, dx, dgrad, HessOld)
309

    
310
       IMPLICIT NONE
311

    
312
       INTEGER, PARAMETER :: KINT=KIND(1)
313
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
314

    
315
       INTEGER(KINT), INTENT(IN) :: n
316
       real(KREAL) :: dx(n), dgrad(n), HessOld(n*n)
317
     END subroutine hupdate_all
318

    
319
     SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess)
320

    
321
       IMPLICIT NONE
322
       
323
       INTEGER, PARAMETER :: KINT=KIND(1)
324
       INTEGER, PARAMETER :: KREAL=KIND(1.D0)
325
       
326
       INTEGER(KINT), INTENT(IN) :: NFREE
327
       REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree)
328
       REAL(KREAL), INTENT(IN) :: ds(Nfree)
329
       REAL(KREAL), INTENT(IN) :: dGrad(Nfree)
330
       
331
     END subroutine Hinvup_BFGS_new
332

    
333

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

    
336
      Use VarTypes
337

    
338
      implicit none
339
! Input
340
      CHARACTER(*), INTENT(IN) :: Line
341
      CHARACTER(*), INTENT(IN) :: String
342
      LOGICAL, OPTIONAL, INTENT(IN) :: Case
343
      LOGICAL, OPTIONAL, INTENT(IN) :: Back
344
      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean
345

    
346
! Output
347
! the position of String in Line (the first one) unless Back is present
348
      INTEGER(KINT) :: Pos
349
    END FUNCTION InString
350

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

    
353
      Use VarTypes
354
      Use io_module
355

    
356
      implicit none
357
      character(len=*), intent(in)           :: routine, msg
358
      character(len=*), intent(in), optional :: file
359
      integer(KINT), intent(in), optional      :: line, unit
360

    
361
    END SUBROUTINE die
362

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

    
365
      Use VarTypes
366
      Use io_module
367

    
368
      implicit none
369

    
370
      character(len=*), intent(in)           :: routine, msg
371
      character(len=*), intent(in), optional :: file
372
      integer(KINT), intent(in), optional      :: line, unit
373

    
374
    END SUBROUTINE Warning
375

    
376

    
377
  END INTERFACE
378

    
379
  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
380
       MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, &
381
       SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, &
382
       input,prog, ProgExe,RunMode, AtTypes,poscar,PathName,  &
383
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
384
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
385
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
386
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile,GplotFile
387

    
388
  Namelist/cartlist/list
389
  Namelist/frozenlist/list
390
  Namelist/analist/nb
391

    
392

    
393
  Flag_Opt_Geom = .FALSE. ! Initialized.
394

    
395
!!!!!!!!!!!!!!!!!!!!!!!!!
396
  ! 
397
  ! We print the Version info in file
398
  !
399
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
400

    
401
  ! We read the files  names
402
!  SELECT CASE (iargc())
403
  SELECT CASE (command_argument_count())
404
  CASE (2)
405
     call getarg(1,FileIn)
406
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
407
     call getarg(2,FileOut)
408
     OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN')
409
  CASE (1)
410
     call getarg(1,FileIn)
411
     IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
412
        WRITE(*,*) "Path mini-help"
413
        WRITE(*,*) "--------------"
414
        WRITE(*,*) ""
415
        WRITE(*,*) "Use: Path Input_file Output_file"
416
        WRITE(*,*) "Input_file starts with a Namelist called path"
417
        WRITE(*,*) ""
418
        WRITE(*,*) "Compulsory variables are:"
419
        WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
420
        WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
421
        WRITE(*,*) "Nat   : Number of atoms"
422
        WRITE(*,*) ""
423
        WRITE(*,*) ""
424
        WRITE(*,*) "Other options"
425
        WRITE(*,*) "Input: string that indicates the type of the input geometries."
426
        WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
427
        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
428
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, TurboMole, Siesta, Test or Ext"
429
        WRITE(*,*) "   *   In case of a Gaussian or Mopac calculations, input must be set to Cart."
430
        WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
431
        WRITE(*,*) "     input file. See example file Test_G03.path or Test_Mopac.path"
432
        WRITE(*,*) "   *   In the case of a VASP calculation, if input is set to Cart, then"
433
        WRITE(*,*) "     the preamble of a VASP calculation should be added at the end of"
434
        WRITE(*,*) "     the input file. See example file Test_VASP_cart.path"
435
        WRITE(*,*) "       In the case of a VASP calculation, one should also give value of the "
436
        WRITE(*,*) "    Runmode variable"
437
        WRITE(*,*) "   *   In the case of a SIESTA calculation, an example of a Siesta input file"
438
        WRITE(*,*) "     should be added at the end of the input file. See Test_Siesta.path"
439
        WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy"
440
        WRITE(*,*) "       and gradient of the whole path or not. If one wants to use VASP,"
441
        WRITE(*,*) "       Runmode must be set to PARA."
442
        WRITE(*,*) "PathName: Prefix used to save the path. Default is Path"
443
        WRITE(*,*) "Poscar: string that will be used as the prefix for the different "
444
        WRITE(*,*) "        POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.,"
445
        WRITE(*,*) "        not used for internal calculations."
446
        WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations"
447
        WRITE(*,*) " ISuffix: Suffix for the input file."
448
        WRITE(*,*) " OSuffix: suffix for the output file."
449
        WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates."
450
        WRITE(*,*) "      Valid only for Coord=Zmat, Hybrid or Mixed"
451
        WRITE(*,*) "Fact: REAL used to define if two atoms are linked."
452
        WRITE(*,*) "      If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked."
453
        WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
454
        WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
455
        WRITE(*,*) "       - CART (or Xyz): works in cartesian"
456
        WRITE(*,*) "       - Zmat: works in internal coordinates (Zmat)"
457
        WRITE(*,*) "       - Mixed: frozen atoms, as well as atoms defined by the "
458
        WRITE(*,*) "       'cart' array(see below) are describe in CARTESIAN, whereas the"
459
        WRITE(*,*) "       others are described in Zmat" 
460
        WRITE(*,*) "       - Baker: use of Baker coordinates, also called delocalized internal coordinates"
461
        WRITE(*,*) "       - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
462
        WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
463
        WRITE(*,*) "       - RFO: Rational function optimization"
464
        WRITE(*,*) "       - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
465
!        WRITE(*,*) " HesUpd: Deprecated. Use Hupdate."
466
        WRITE(*,*) " HUpdate: method to update the hessian. By default, it is Murtagh-Sargent"
467
        WRITE(*,*) "      Except for geometry optimization where it is BFGS."
468
        WRITE(*,*) " GeomFile: Name of the file to print the geometries and their energy"
469
        WRITE(*,*) "      during a geometry optimization. If this variable is not given"
470
        WRITE(*,*) "      then nothing is printed."
471
        WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
472
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
473
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
474
        WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
475
        WRITE(*,*) "FTan: When moving the path, this gives the proportion of the displacement tangent to the path"
476
        WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
477
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path. Default: Ftan=0."
478
        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
479
        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
480
        WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
481
        WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, "
482
        WRITE(*,*) "        it might happen that an atom moves out of the unit cell."
483
        WRITE(*,*) "        Path detects this by comparing the displacement to boxtol:"
484
        WRITE(*,*) "        if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5"
485
        WRITE(*,*) ""
486
        WRITE(*,*) "Arrays:"
487
        WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements."
488
        WRITE(*,*) "      You can modify it using, rcov(6)=0.8."
489
        WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements."
490
        WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations."
491
        WRITe(*,*) "If not given, Path will read the POTCAR file."
492
        WRITE(*,*) ""
493
        WRITE(*,*) "Flags:"
494
        WRITE(*,*) "MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
495
        WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time."
496
        WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
497
        WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
498
        WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para"
499
        WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para"
500
        WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
501
        WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
502
        WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
503
        WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path."
504
        WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
505
        WRITE(*,*) "         If True, a &frozenlist namelist containing the list of frozen atoms must be given."
506
        WRITE(*,*) "          If VASP is used, and frozen is not given"
507
        WRITE(*,*) " here, the program will use the F flags of the input geometry"
508
        WRITE(*,*) "FCart:  True if one wants to describe some atoms using cartesian coordinates. "
509
        WRITE(*,*) "        *** Only used in 'mixed' calculations. ***"
510
        WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
511
        WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
512
        WRITE(*,*) ""
513
        WRITE(*,*) "AnaGEom: If TRUE, Opt'n Path will create a file .dat for geometries analysis."
514
        WRITE(*,*) "        If True, Opt'n Path will look for the AnaList namelist after the Path Namelist."
515
        WRITE(*,*) "AnaList: list of variables for geometry analysis. If present and if AnaGeom=T then"
516
        WRITE(*,*) "     Opt'n Path will read it and print the values of the variable in a .dat file"
517
        WRITE(*,*) "    If AnaGeom is T but Analist is not given, then Path.dat just contains the energy." 
518
        WRITE(*,*) "Analist contains the number of variables to monitor: Nb, and is followed by the description"
519
        WRITE(*,*) "of the variables among:"
520
        WRITE(*,*) "b(ond) At1 At2"
521
        WRITE(*,*) "a(ngle) At1 At2 At3"
522
        WRITE(*,*) "d(ihedral) At1 At2 At3 At4"
523
        WRITE(*,*) "c NbAt At1 At2 At3 At4 At5... to create a center of mass. The centers of mass are added"
524
        WRITE(*,*) "at the end of the real atoms of the system"
525
        WRITE(*,*) ""
526
        WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
527
        WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
528
        WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
529
        WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
530
        WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now"
531
        WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files."
532
        WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP."
533
        WRITE(*,*) ""
534
        STOP
535
     END IF
536

    
537
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
538
     IOOUT=6
539
  CASE (0)
540
     IOIN=5
541
     IOOUT=6
542
  END SELECT
543

    
544

    
545
  ! We initiliaze variables
546
  Pi=dacos(-1.0d0)
547
  Nat=1
548
  Fact=1.1
549
  IGeomRef=-1
550
  NGeomI=1
551
  NGeomF=8
552
  IReparam=5
553
  IReparamT=-1
554
  ISpline=5
555
  debug=.False.
556
  MW=.TRUE.
557
  bohr=.False.
558
  Coord='ZMAT'
559
  Step_method='RFO'
560
  DebugFile='Path.valid'
561
  PathName="Path"
562
  MaxCyc=50
563
  SMax=0.1D0
564
  SThresh=-1.
565
  GThresh=-1.
566
  FTan=0.
567
  Hinv=.TRUE.
568
  IniHUp=.FALSE.
569
  HupNeighbour=.FALSE.
570
  HesUpd="EMPTY"
571
  HUpdate="EMPTY"
572
  FFrozen=.FALSE.
573
  FCart=.FALSE.
574
  List=0
575
  renum=.TRUE.
576
  OptReac=.FALSE.
577
  OptProd=.FALSE.
578
  CalcEReac=.FALSE.
579
  CalcEProd=.FALSE.
580
  EReac=0.d0
581
  EProd=0.d0 
582
  OptGeom=-1
583
  GeomFile="EMPTY"
584
  AnaGeom=.FALSE.
585
  AnaFile="EMPTY"
586
  GplotFile="EMPTY"
587
  Nb=0
588
  NbCom=0
589
  PathOnly=.FALSE.
590
  Prog='EMPTY'
591
  ProgExe='EMPTY'
592
  Runmode='SERIAL'
593
  Input='EMPTY'
594
  Poscar="POSCAR"
595
  AutoCart=.FALSE.
596
  VMD=.TRUE.
597
  WriteVASP=.FALSE.
598
  DynMaxStep=.FALSE.
599
  CalcName="EMPTY"
600
  ISuffix="EMPTY"
601
  OSuffix="EMPTY"
602
  NGintMax=10
603
  Align=.TRUE.
604

    
605
  ! Number of point used to compute the length of the path,
606
  ! and to reparameterize the path
607
  NMaxPtPath=-1 
608
  FrozTol=1e-4
609

    
610
  ! Read the namelist
611
  read (IOIN,path)
612

    
613
  debug=valid("Main").or.valid("Path")
614

    
615
  FCalcHess=valid("CalcHess")
616
  PathName=AdjustL(PathName)
617
  Coord=AdjustL(Coord)
618
  CALL UpCase(coord)
619

    
620
  Step_method=AdjustL(Step_method)
621
  CALL UpCase(Step_method)
622

    
623
  Prog=AdjustL(Prog)
624
  CALL UpCase(prog)
625

    
626
  Input=AdjustL(Input)
627
  CALL UpCase(input)
628

    
629
  Runmode=AdjustL(Runmode)
630
  Call UpCase(Runmode)
631
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"
632

    
633
  if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN
634
     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"
635
     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"
636
     CalcEReac=.FALSE.
637
     CalcEProd=.FALSE.
638
  END IF
639

    
640
  If (IReparamT<=0) IReparamT=IReparam
641

    
642
! We take care of HesUpd flag
643
! this is needed because some users might still use it
644
  if (HesUpd.NE."EMPTY") THEN
645
! We are pityless:
646
     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"
647
     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"
648
     STOP
649
  END IF
650

    
651
  IF (HUpdate=='EMPTY') THEN
652
     IF (OptGeom>=1) THEN
653
        HupDate="BFGS"
654
     ELSE
655
        HUpdate="MS"
656
     END IF
657
  END IF
658

    
659
  If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 
660
     COORD='CART'
661
  END IF
662

    
663
  IF (COORD.EQ.'CART') THEN
664
     Renum=.FALSE.
665
     IF (OptGeom>0) THEN
666
        IGeomRef=OptGeom
667
     ELSE
668
        IGeomRef=1
669
     END IF
670
  END IF
671

    
672
  
673
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
674
  if (Prog.EQ.'2D') Prog="TEST2D"
675
  if (Prog.EQ.'TEST3') Prog="LEPS"
676

    
677

    
678
  Select Case (Prog)
679
    CASE ("EMPTY","G09","GAUSSIAN")
680
     Prog='GAUSSIAN'
681
     If (ProgExe=="EMPTY") ProgExe="g09"
682
     if (CalcName=="EMPTY") CalcName="Path"
683
     if (ISuffix=="EMPTY")  ISuffix="com"
684
     if (OSuffix=="EMPTY")  OSuffix="log"
685
     UnitProg="au"
686
     ConvE=au2kcal
687
    CASE ("G03")
688
     Prog='GAUSSIAN'
689
     If (ProgExe=="EMPTY") ProgExe="g03"
690
     if (CalcName=="EMPTY") CalcName="Path"
691
     if (ISuffix=="EMPTY")  ISuffix="com"
692
     if (OSuffix=="EMPTY")  OSuffix="log"
693
     UnitProg="au"
694
     ConvE=au2kcal
695
    CASE ("EXT")
696
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
697
     if (CalcName=="EMPTY") CalcName="Ext"
698
     if (ISuffix=="EMPTY")  ISuffix="in"
699
     if (OSuffix=="EMPTY")  OSuffix="out"
700
     UnitProg="au"
701
     ConvE=au2kcal
702
    CASE ('MOPAC')
703
     If (ProgExe=="EMPTY") ProgExe="mopac"
704
     if (CalcName=="EMPTY") CalcName="Path"
705
     if (ISuffix=="EMPTY")  ISuffix="mop"
706
     if (OSuffix=="EMPTY")  OSuffix="out"
707
     UnitProg="au"
708
     ConvE=au2kcal
709
    CASE ("SIESTA")
710
     If (ProgExe=="EMPTY") ProgExe="siesta"
711
     if (CalcName=="EMPTY") CalcName="Path"
712
     if (ISuffix=="EMPTY")  ISuffix="fdf"
713
     if (OSuffix=="EMPTY")  OSuffix="out"
714
     UnitProg="eV"
715
     ConvE=eV2kcal
716
    CASE ("VASP")
717
     If (ProgExe=="EMPTY") ProgExe="VASP"
718
     UnitProg="eV"
719
     ConvE=eV2kcal
720
    CASE ("TM","TURBOMOLE")
721
     Prog="TURBOMOLE"
722
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
723
     UnitProg="au"
724
     ConvE=au2kcal
725
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
726
     ProgExe=""
727
     UnitProg="au"
728
     ConvE=au2kcal
729
    CASE DEFAULT
730
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
731
     STOP
732
  END SELECT 
733

    
734
  if (Input.EQ.'EMPTY') THEN
735
     Select Case (Prog)
736
       CASE ("VASP")
737
          Input=Prog
738
       CASE ("CHAMFRE")
739
          Input="VASP"
740
       CASE DEFAULT
741
          Input='XYZ'
742
     END SELECT
743
    WRITE(*,*) "Input has been set to the default: ",INPUT
744
 END IF
745

    
746
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
747

    
748
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
749
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
750
     STOP
751
  END IF
752

    
753
  if (OptGeom.GE.1) THEN
754
     FPrintGeom=.FALSE.
755
     If (GeomFile/='EMPTY') THEN
756
        FPrintGeom=.TRUE.
757
     END IF
758
     IF (IGeomRef.LE.0) THEN
759
        IGeomRef=OptGeom
760
     ELSE
761
        IF (IGeomRef.NE.OptGeom) THEN
762
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
763
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
764

    
765
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
766
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
767
           STOP
768
        END IF
769
     END IF
770
  END IF
771

    
772
  ALLOCATE(Cart(Nat))
773
  Cart=0
774

    
775
  ALLOCATE(Frozen(Nat))
776
  Frozen=0
777

    
778
  IF (FCart) THEN
779
! We rewind the file to be sure to browse all of it to read the Cart namelist
780
     REWIND(IOIN)
781
     List=0
782
     READ(IOIN,cartlist)
783
     IF (COORD.NE.'CART') THEN
784
        Cart=List(1:Nat)
785
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
786
     ELSE
787
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
788
        WRITE(*,*) "I will discard the cartlist"
789
        Cart=0
790
     END IF
791
  END IF
792

    
793
  if (FFrozen) THEN
794

    
795
     REWIND(IOIN)
796
     List=0
797
     READ(IOIN,Frozenlist)
798
     Frozen=List(1:Nat)
799
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
800
  END IF
801

    
802
  If (FCart) Call Expand(Cart,NCart,Nat)
803
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
804

    
805
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
806
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
807

    
808
  if (AnaGeom) THEN
809
     REWIND(IOIN)
810
     READ(IOIN,AnaList,IOSTAT=IoS)
811
     IF (IOS==0) THEN ! we have a AnaList namelist 
812
        Call ReadAnaList
813
     ELSE
814
! No AnaList namelist, we will print only the energy
815
        FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)'
816
     END IF
817
     IF (AnaFile=="EMPTY") THEN
818
        AnaFile=TRIM(PathName) // '.dat'
819
     END IF
820
     OPEN(IODat,File=AnaFile)
821
     IF (GplotFile=="EMPTY") THEN
822
        GplotFile=TRIM(PathName) // '.gplot'
823
     END IF
824
     
825
     OPEN(IODat,File=AnaFile)
826
     If (OptGeom<=0) THEN
827
        OPEN(IOGplot,File=GPlotFile)
828
        WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist"
829
        WRITE(IoGplot,'(A)') " set pointsize 2"
830
        WRITE(IoGplot,'(A)') " xcol=1"
831
        WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3
832
     END IF
833
  END IF
834

    
835
        
836

    
837
  if (NMaxPtPath.EQ.-1) then
838
     NMaxPtPath=min(50*NGeomF,2000)
839
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
840
  end if
841

    
842
  !  NMaxPtPath=10000
843

    
844
  ! We ensure that FTan is in [0.:1.]
845
  FTan=min(abs(FTan),1.d0)
846

    
847
! PFL 2011 Mar 14 ->
848
! Added some printing for analyses with Anapath
849
  if (debug) THEN
850
     WRITE(IOOUT,path)
851
  ELSE
852
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC
853
! and PATHNAME 
854
     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat
855
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI
856
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF
857
     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC
858
     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)
859
  END IF
860

    
861
  FTan=1.-FTan
862

    
863
  !Thresholds...
864
  if (SThresh.LE.0) SThresh=0.01d0
865
  If (GThresh.LE.0) GThresh=SThresh/4.
866
  SThreshRms=Sthresh*2./3.
867
  GThreshRms=Gthresh*2./3.
868

    
869
  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
870
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
871
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
872
  ALLOCATE(MassAt(NAt),OrderInv(Nat))
873

    
874
  AtName=""
875
  MassAt=1.
876

    
877
! As we have used rewind many times, the position in the input file
878
! is blurry (at least !) we thus have to go back to the Geometries description
879
! (TO DO: use a better Input parser !)
880

    
881
  REWIND(IOIN)
882
! We search the '&path' namelist
883
  READ(IOIN,'(A)',Iostat=Ios) Line
884
  Fini=.FALSE.
885
  DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0)
886
     if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line)
887
     READ(IOIN,'(A)',Iostat=Ios) Line
888
  END DO
889
  if (debug) WRITE(*,*) "Path found:",TRIM(LINE)
890
  ! We are on the &path line, we move to the end
891
  Call NoComment(Line)
892
  DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
893
     if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
894
     READ(IOIN,'(A)',Iostat=Ios) Line
895
     Call NoComment(Line)
896
     Call noString(Line)
897
  END DO
898
! We are at then end of &Path / namelist
899
! We scan for other Namelists
900
  READ(IOIN,'(A)',Iostat=Ios) Line
901
  if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
902
  DO WHILE  ((IoS==0).AND.(Index(Line,'&')/=0))
903
     Call NoComment(Line)
904
     if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
905
     Idx=InString(Line,'Analist')
906
     DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
907
        if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
908
        READ(IOIN,'(A)',Iostat=Ios) Line
909
        Call NoComment(Line)
910
     END DO
911
     If (Idx/=0) THEN
912
! we have just read &Analist namelist, we have to skip the variables description
913
        DO I=1,Nb
914
           READ(IOIN,'(A)',Iostat=Ios) Line
915
           if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE)
916
        END DO
917
     END IF
918
     READ(IOIN,'(A)',Iostat=Ios) Line
919
     if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE)
920
  END DO
921
  BACKSPACE(IOIN)
922

    
923
  ! We read the initial geometries
924
  Call Read_Geom(input)
925

    
926
  ! We convert atome names into atom mass number
927
  DO I=1,NAt
928
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
929
     Atome(I)=ConvertNumAt(AtName(I))
930
  END DO
931

    
932
 If (AnaGeom) THEN
933
    If (NbVar>0) THEN
934
! We print the description of the varialbes in AnaFile
935
       Call PrintAnaList(IoDat)
936
       If (OptGeom<=0) THEN
937
          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
938
       ELSE
939
          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
940
       END IF
941
    ELSE
942
       If (OptGeom<=0) THEN
943
          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
944
       ELSE
945
          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"       
946
       END IF
947
    END IF
948
 END IF
949

    
950
  ! PFL 23 Sep 2008 ->
951
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
952
  ! between the first geometry and the others.
953
  IF (NFroz.GT.0) THEN
954
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
955
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
956
     ListAtFroz=.FALSE.
957

    
958
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
959
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
960
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
961

    
962
     SELECT CASE (NFroz)
963
        ! We have Frozen atoms
964
        ! PFL 28 Oct 2008 ->
965
        ! It might happen that the geometries are translated/rotated.
966
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
967
        ! three frozen atoms.
968

    
969
     CASE (3:)
970
        DO I=1,NFroz
971
           ListAtFroz(Frozen(I))=.TRUE.
972
        END DO
973
     CASE (2)
974
        IAt1=Frozen(1)
975
        IAt2=Frozen(2)
976
        ListAtFroz(Iat1)=.TRUE.
977
        ListAtFroz(Iat2)=.TRUE.
978
        x2(1)=x1(Iat2)-x1(Iat1)
979
        y2(1)=y1(Iat2)-y1(Iat1)
980
        z2(1)=z1(Iat2)-z1(Iat1)
981
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
982
        Dist=-1.
983
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
984
        DO I=1,Nat
985
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
986
           x2(2)=x1(Iat2)-x1(I)
987
           y2(2)=y1(Iat2)-y1(I)
988
           z2(2)=z1(Iat2)-z1(I)
989
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
990
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
991
           if (abs(ca)<0.9) THEN
992
              IF (Norm2>Dist) THEN
993
                 Iat3=I
994
                 Dist=Norm2
995
              END IF
996
           END IF
997
        END DO
998
        ListAtFroz(Iat3)=.TRUE.
999
        if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3
1000
     CASE (1)
1001
        IAt1=Frozen(1)
1002
        ListAtFroz(Iat1)=.TRUE.
1003
        Dist=-1.
1004
        ! We scan all atoms to find one that is further from At1
1005
        DO I=1,Nat
1006
           if (I.EQ.Iat1) Cycle
1007
           x2(1)=x1(Iat1)-x1(I)
1008
           y2(1)=y1(Iat1)-y1(I)
1009
           z2(1)=z1(Iat1)-z1(I)
1010
           Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1011
           IF (Norm1>Dist) THEN
1012
              Iat2=I
1013
              Dist=Norm1
1014
           END IF
1015
        END DO
1016
        ListAtFroz(Iat2)=.TRUE.
1017
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
1018
        x2(1)=x1(Iat2)-x1(Iat1)
1019
        y2(1)=y1(Iat2)-y1(Iat1)
1020
        z2(1)=z1(Iat2)-z1(Iat1)
1021
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1022
        Dist=-1.
1023
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1024
        Iat3=1
1025
        DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat))
1026
           Iat3=Iat3+1
1027
        END DO
1028
        DO I=1,Nat
1029
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
1030
           x2(2)=x1(Iat2)-x1(I)
1031
           y2(2)=y1(Iat2)-y1(I)
1032
           z2(2)=z1(Iat2)-z1(I)
1033
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
1034
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
1035
           if (abs(ca)<0.9) THEN
1036
              IF (Norm2>Dist) THEN
1037
                 Iat3=I
1038
                 Dist=Norm2
1039
              END IF
1040
           END IF
1041
        END DO
1042
        ListAtFroz(Iat3)=.TRUE.
1043
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3
1044
     END SELECT
1045

    
1046
     DO I=2,NGeomI
1047
        ! First, we align the Ith geometry on the first one, using only
1048
        ! the positions of the "frozen" atoms. (see above: it might be that
1049
        ! ListAtFroz contains non frozen atoms used to align the geometries
1050
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
1051
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
1052
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
1053
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
1054

    
1055
        FChkFrozen=.FALSE.
1056
        DO J=1,NFroz
1057
           IAt=Frozen(J)
1058
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
1059
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
1060
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
1061
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
1062
              FChkFrozen=.TRUE.
1063
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
1064
              WRITE(*,*) "from geometry 1 to geometry ",I
1065
           END IF
1066
        END DO
1067
     END DO
1068
     IF (FChkFrozen) THEN
1069
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
1070
        STOP
1071
     END IF
1072
  END IF
1073

    
1074

    
1075
  N3at=3*Nat
1076
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
1077

    
1078

    
1079
  Call ReadInput
1080

    
1081
  IF (COORD=="MIXED") Call TestCart(AutoCart)
1082

    
1083
  SELECT CASE(NFroz)
1084
  CASE (2)
1085
     IntFroz=1
1086
  CASE (3)
1087
     IntFroz=3
1088
  CASE (4:)
1089
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
1090
  END SELECT
1091

    
1092
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
1093
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
1094
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
1095

    
1096
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
1097
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
1098
     COORD="ZMAT"
1099
  END IF
1100

    
1101
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
1102

    
1103
  SELECT CASE(COORD)
1104
  CASE ('ZMAT')
1105
     NCoord=Nfree
1106
     ALLOCATE(dzdc(3,nat,3,nat))
1107
  CASE ('MIXED')
1108
     SELECT CASE (NCart+NFroz)
1109
     CASE (1) 
1110
        NCoord=N3At-3
1111
     CASE (2)
1112
        NCoord=N3At-1
1113
     CASE (3:)
1114
        NCoord=N3At
1115
     END SELECT
1116
     ALLOCATE(dzdc(3,nat,3,nat))
1117
  CASE ('BAKER')
1118
     Symmetry_elimination=0
1119
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
1120
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
1121
     ALLOCATE(BTransInv(NCoord,3*Nat))
1122
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
1123
  CASE ('HYBRID')
1124
     NCoord=N3at
1125
  CASE ('CART')
1126
     NCoord=N3at
1127
  END SELECT
1128

    
1129
  if (debug) WRITE(*,*) "DBG Path, L826, Coord, NCart,NCoord, N3At",TRIM(COORD),Ncart, NCoord,N3at
1130

    
1131
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
1132
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
1133
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
1134
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
1135
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
1136
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
1137
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
1138
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
1139

    
1140
  ZeroVec=0.d0
1141
  DO I =1, NGeomF
1142
     IntTangent(I,:)=0.d0
1143
  END DO
1144

    
1145
  if (debug) THEN
1146
     Print *, 'L886, IntTangent(1,:)='
1147
     Print *, IntTangent(1,:)
1148
  END IF
1149

    
1150
  if (.NOT.OptReac) Energies(1)=Ereac
1151
  if (.NOT.OptProd) Energies(NGeomF)=EProd
1152
  MaxStep=SMax
1153

    
1154
  ! We initialize the mass array
1155
  if (MW) THEN
1156
     WRITE(*,*) 'Working in MW coordinates'
1157
     DO I=1,Nat
1158
        massAt(I)=Mass(Atome(I))
1159
     END DO
1160
  ELSE
1161
     DO I=1,Nat
1162
        MassAt(I)=1.d0
1163
     END DO
1164
  END IF
1165

    
1166
  WRITE(*,*) "Prog=",TRIM(PROG)
1167

    
1168
  ALLOCATE(FrozAtoms(Nat))
1169

    
1170
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
1171
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
1172
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
1173
  IF (NFroz.EQ.0) THEN
1174
     FrozAtoms=.TRUE.
1175
  ELSE
1176
     I=1
1177
     NFr=0
1178
     FrozAtoms=.False.
1179
     DO WHILE (Frozen(I).NE.0)
1180
        IF (Frozen(I).LT.0) THEN
1181
           DO J=Frozen(I-1),abs(Frozen(I))
1182
              IF (.NOT.FrozAtoms(J)) THEN
1183
                 NFr=NFr+1
1184
                 FrozAtoms(J)=.TRUE.
1185
              END IF
1186
           END DO
1187
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1188
           FrozAtoms(Frozen(I))=.TRUE.
1189
           NFr=NFr+1
1190
        END IF
1191
        I=I+1
1192
     END DO
1193
     IF (NFr.NE.NFroz) THEN
1194
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1195
        STOP
1196
     END IF
1197
  END IF
1198

    
1199
!  if (debug) THEN
1200
     !Print *, 'L968, IntTangent(1,:)='
1201
     !Print *, IntTangent(1,:)
1202
!  END IF
1203

    
1204
  ! We have read everything we needed... time to work :-)
1205
  IOpt=0
1206
  FirstTimePathCreate = .TRUE. ! Initialized.
1207
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1208
  FirstTimePathCreate = .FALSE.
1209

    
1210
  if (debug) THEN
1211
     Print *, 'L980, IntTangent(1,:)='
1212
     Print *, IntTangent(1,:)
1213
  END IF
1214

    
1215
  ! PFL 30 Mar 2008 ->
1216
  ! The following is now allocated in Calc_Baker. This is more logical to me
1217
  !   IF (COORD .EQ. "BAKER") Then
1218
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1219
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1220
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1221
  !   END IF
1222
  ! <- PFL 30 mar 2008
1223

    
1224
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1225
  !
1226
  ! For Debugging purposes
1227
  !
1228
  ! OptGeom is the index of the geometry which is to be optimized.
1229
  IF (Optgeom.GT.0) THEN
1230
     Flag_Opt_Geom=.TRUE.
1231
     SELECT CASE(Coord)
1232
     CASE ('CART','HYBRID')
1233
!!! CAUTION : PFL 29.AUG.2008 ->
1234
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1235
        ! This should be made  consistent !!!!!!!!!!!!!!!
1236
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1237
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1238
!!! <- CAUTION : PFL 29.AUG.2008
1239
     CASE ('ZMAT','MIXED')
1240
        !GeomTmp=IntCoordF(OptGeom,:)
1241
        GeomTmp=IntCoordI(OptGeom,:)
1242
     CASE ('BAKER')
1243
        !GeomTmp=IntCoordF(OptGeom,:)
1244
        GeomTmp=IntCoordI(OptGeom,:)
1245
     CASE DEFAULT
1246
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1247
        STOP
1248
     END SELECT
1249

    
1250
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1251
     Flag_Opt_Geom = .TRUE.
1252
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1253

    
1254
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1255
     WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1256
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1257
     STOP
1258
  END IF ! IF (Optgeom.GT.0) THEN
1259

    
1260
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1261
  ! End of GeomOpt
1262
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1263

    
1264
  IF (PathOnly) THEN
1265
     WRITE(*,*) "PathOnly=.T. , Stop here"
1266
     Call Write_path(-1)
1267
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1268
     STOP
1269
  END IF
1270

    
1271
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1272

    
1273
  DEALLOCATE(XyzGeomI,IntCoordI)
1274
  NGeomI=NGeomF
1275
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1276

    
1277
  IF (DEBUG) WRITE(*,*) ' Back from PathCreate, L965, Path.f90'
1278
  ! we print this path, which is the first one :-) called Iteration0
1279
  ! we now have to calculate the gradient for each point (and the energy
1280
  ! if we work at 0K)
1281

    
1282
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1283
     DO IGeom=1,NGeomF
1284
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1285
        DO I=1,3*Nat-6
1286
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1287
        END DO
1288
     END DO
1289
  END IF
1290

    
1291
  ! To calculate B^T^-1 for all extrapolated final geometries:
1292
  ! PFL 18 Jan 2008 ->
1293
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1294
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1295
  ! <-- PFL 18 Jan 2008
1296

    
1297
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1298
     DO IGeom=1,NGeomF
1299
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
1300
        DO I=1,3*Nat-6
1301
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1302
        END DO
1303
     END DO
1304
  END IF
1305

    
1306
  IOpt=0
1307
  Call EgradPath(IOpt)
1308

    
1309
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1310
  Call Write_path(IOpt)
1311
! Shall we analyze the geometries ?
1312
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1313

    
1314
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1315

    
1316
  ! We have the geometries, the first gradients... we will generate the first hessian matrices :-)
1317
  ! ... v3.71 Only if IniHUp=.TRUE. which is not the default. In this version, we take the hessian
1318
  ! of end points as Unity matrices, and we update them along the path.
1319

    
1320
  ! By default, Hess=Id
1321
  Hess=0.
1322
  IF(Coord .EQ. "ZMAT") Then
1323
     ! We use the fact that we know that approximate good hessian values
1324
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1325
     DO IGeom=1,NGeomF
1326
        Hess(IGeom,1,1)=0.5d0
1327
        Hess(IGeom,2,2)=0.5d0
1328
        Hess(IGeom,3,3)=0.1d0
1329
        DO J=1,Nat-3
1330
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1331
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1332
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1333
        END DO
1334
        IF (HInv) THEN
1335
           DO I=1,NCoord
1336
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1337
           END DO
1338
        END IF
1339
     END DO
1340
  ELSE
1341
     DO I=1,NCoord
1342
        DO J=1,NGeomF
1343
           Hess(J,I,I)=1.d0
1344
        END DO
1345
     END DO
1346
  END IF
1347

    
1348
  IF (COORD .EQ. "BAKER") THEN
1349
     ! UMat(NPrim,NCoord)
1350
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1351
     Hprim=0.d0
1352
     ScanCoord=>Coordinate
1353
     I=0
1354
     DO WHILE (Associated(ScanCoord%next))
1355
        I=I+1
1356
        SELECT CASE (ScanCoord%Type)
1357
        CASE ('BOND')
1358
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1359
           Hprim(I,I) = 0.5d0
1360
           !END DO
1361
        CASE ('ANGLE')
1362
           Hprim(I,I) = 0.2d0
1363
        CASE ('DIHEDRAL')
1364
           Hprim(I,I) = 0.1d0
1365
        END SELECT
1366
        ScanCoord => ScanCoord%next
1367
     END DO
1368

    
1369
     ! HprimU = H^prim U:
1370
     HprimU=0.d0
1371
     DO I=1, 3*Nat-6
1372
        DO J=1,NPrim
1373
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1374
        END DO
1375
     END DO
1376

    
1377
     ! Hess = U^T HprimU = U^T H^prim U:
1378
     Hess=0.d0
1379
     DO I=1, 3*Nat-6
1380
        DO J=1,NPrim
1381
           ! UMat^T is needed below.
1382
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1383
        END DO
1384
     END DO
1385
     DO K=2,NGeomF
1386
        Hess(K,:,:)=Hess(1,:,:)
1387
     END DO
1388
     DEALLOCATE(Hprim,HprimU)
1389
  end if !  ends IF COORD=BAKER
1390
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1391
  ALLOCATE(HessTmp(NCoord,NCoord))
1392

    
1393
  IF (IniHUp) THEN
1394
     IF (FCalcHess) THEN
1395
        ! The numerical calculation of the Hessian has been switched on
1396
        DO IGeom=1,NGeomF
1397
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1398
              GeomTmp=IntCoordF(IGeom,:)
1399
           ELSE
1400
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1401
           END IF
1402
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1403
        END DO
1404
     ELSE
1405
        ! We generate a forward hessian and a backward one and we mix them.
1406
        ! First, the forward way:
1407
        DO IGeom=2,NGeomF-1
1408
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1409
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1410
           ELSE
1411
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1412
           END IF
1413
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1414
           HessTmp=Hess(IGeom-1,:,:)
1415
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1416
           Hess(IGeom,:,:)=HessTmp
1417
        END DO
1418

    
1419
        ! Now, the backward way:
1420
        HessTmp=Hess(NGeomF,:,:)
1421
        DO IGeom=NGeomF-1,2,-1
1422
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1423
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1424
           ELSE
1425
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1426
           END IF
1427
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1428
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1429
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1430
           ca=cos(alpha)
1431
           sa=sin(alpha)
1432
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1433
        END DO
1434
     END IF ! matches IF (FCalcHess)
1435
  END IF ! matches IF (IniHUp) THEN
1436

    
1437
  ! For debug purposes, we diagonalize the Hessian matrices
1438
  IF (debug) THEN
1439
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1440
     DO I=1,NGeomF
1441
        WRITE(*,*) "DBG Main - Point ",I
1442
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1443
        DO J=1,NCoord
1444
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1445
        END DO
1446
     END DO
1447
  END IF ! matches IF (debug) THEN
1448

    
1449
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1450
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1451
  !
1452
  ! Main LOOP : Optimization of the Path
1453
  !
1454
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1455
  if (debug)   Print *, 'NGeomF=', NGeomF
1456
  allocation_flag = .TRUE.
1457

    
1458
  Fini=.FALSE.
1459

    
1460
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1461
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1462
     IOpt=IOpt+1
1463

    
1464
     SELECT CASE (COORD)
1465
     CASE ('ZMAT','MIXED','BAKER')
1466
        GeomOld_all=IntCoordF
1467
     CASE ('CART','HYBRID')
1468
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1469
     CASE DEFAULT
1470
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1471
        STOP
1472
     END SELECT
1473

    
1474
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1475
          valid("printtangent")) THEN
1476
        WRITE(*,*) "Visualization of tangents"
1477
        Idx=min(12,NCoord)
1478
        DO I=1,NGeomF
1479
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1480
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1481
           WRITE(*,*) 
1482
        END DO
1483
        WRITE(*,*) "END of tangents"
1484
     END IF
1485

    
1486
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1487
          valid("printgrad")) THEN
1488
        WRITE(*,*) "Visualization of gradients"
1489
        Idx=min(12,NCoord)
1490
        DO I=1,NGeomF
1491
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1492
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1493
           WRITe(*,*) 
1494
        END DO
1495
        WRITE(*,*) "END of gradients"
1496
     END IF
1497

    
1498
     Fini=.TRUE.
1499
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1500
        IGeom=1
1501
        if (debug) WRITE(*,*) '**************************************'
1502
        if (debug) WRITE(*,*) 'Optimizing reactant'
1503
        if (debug) WRITE(*,*) '**************************************'
1504
        SELECT CASE (COORD)
1505
        CASE ('ZMAT','MIXED','BAKER')
1506
           SELECT CASE (Step_method)
1507
           CASE ('RFO')
1508
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
1509
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1510
              Print *, GradTmp
1511
           CASE ('GDIIS')
1512
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1513
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1514
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1515
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1516
              Print *, GradTmp
1517
           CASE ('GEDIIS')
1518
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1519
              ! Energies are updated in EgradPath.f90
1520
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1521
                   NCoord,allocation_flag,ZeroVec)
1522
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1523
              !GradTmp = GeomTmp - IntCoordF(1,:)
1524
              !Print *, 'Old Geometry:'
1525
              !Print *, IntCoordF(1,:)
1526
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1527
              Print *, GradTmp
1528
              !Print *, 'GeomTmp:'
1529
              !Print *, GeomTmp
1530
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1531
           CASE DEFAULT
1532
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1533
              STOP
1534
           END SELECT
1535
        CASE ('HYBRID','CART')
1536
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1537
        CASE DEFAULT
1538
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1539
           STOP
1540
        END SELECT
1541

    
1542
        IF (debug) THEN
1543
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1544
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1545
        END IF
1546

    
1547
        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
1548
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1549
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
1550
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1551
        OptReac=(NormStep.GT.SThresh)
1552
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
1553
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
1554
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1555
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
1556
        !Print *, 'Grad(1,:):'
1557
        !Print *, Grad(1,:)
1558
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac
1559

    
1560
        GradTmp=GradTmp*FactStep
1561

    
1562
        ! we take care that frozen atoms do not move.
1563
        IF (NFroz.NE.0) THEN
1564
           SELECT CASE (COORD)
1565
           CASE ('ZMAT','MIXED')
1566
              GradTmp(1:IntFroz)=0.d0
1567
           CASE ('CART','HYBRID')
1568
              DO I=1,Nat
1569
                 IF (FrozAtoms(I)) THEN
1570
                    Iat=Order(I)
1571
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1572
                 END IF
1573
              END DO
1574
           CASE('BAKER')
1575
              GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..)
1576
           CASE DEFAULT
1577
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP'
1578
              STOP
1579
           END SELECT
1580
        END IF ! matches IF (NFroz.NE.0) THEN
1581

    
1582
        IF (debug) THEN
1583
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1584
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1585
        END IF
1586

    
1587
        IF (debug) THEN
1588
           !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1'
1589
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1590
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1591
           ELSE
1592
              DO Iat=1,Nat
1593
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1594
                      XyzGeomF(IGeom,1:3,Iat)
1595
              END DO
1596
           END IF
1597
        END IF ! matches IF (debug) THEN
1598

    
1599
        SELECT CASE (COORD)
1600
        CASE ('ZMAT','MIXED','BAKER')
1601
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1602
        CASE ('HYBRID','CART')
1603
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1604
        CASE DEFAULT
1605
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1606
           STOP
1607
        END SELECT
1608

    
1609
        IF (debug) THEN
1610
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1611
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1612
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1613
           ELSE
1614
              DO Iat=1,Nat
1615

    
1616
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1617
                 ! It was with indzmat. IndZmat cannot appear here...
1618
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1619
                      XyzGeomF(IGeom,1:3,Iat)
1620
              END DO
1621
           END IF
1622
        END IF ! matches IF (debug) THEN
1623

    
1624
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1625
     ELSE ! Optreac
1626
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1627
     END IF ! matches IF (OptReac) THEN
1628

    
1629
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1630
        IGeom=NGeomF
1631
        if (debug) WRITE(*,*) '******************'
1632
        if (debug) WRITE(*,*) 'Optimizing product'
1633
        if (debug) WRITE(*,*) '******************'
1634
        SELECT CASE (COORD)
1635
        CASE ('ZMAT','MIXED','BAKER')
1636
           Print *, 'Optimizing product'
1637
           SELECT CASE (Step_method)
1638
           CASE ('RFO')
1639
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1640
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1641
              Print *, GradTmp
1642
           CASE ('GDIIS')
1643
              ! GradTmp becomes "step" below and has INTENT(OUT):
1644
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1645
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1646
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1647
              Print *, GradTmp
1648
           CASE ('GEDIIS')
1649
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1650
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1651
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1652
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1653
              Print *, GradTmp
1654
           CASE DEFAULT
1655
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1656
              STOP
1657
           END SELECT
1658
        CASE ('HYBRID','CART')
1659
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1660
        CASE DEFAULT
1661
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1662
           STOP
1663
        END SELECT
1664

    
1665
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1666
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1667
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1668
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1669
        OptProd=.NOT.(NormStep.LE.SThresh)
1670
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1671
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1672
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1673
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1674

    
1675
        GradTmp=GradTmp*FactStep
1676

    
1677
        ! we take care that frozen atoms do not move
1678
        IF (NFroz.NE.0) THEN
1679
           SELECT CASE (COORD)
1680
           CASE ('ZMAT','MIXED','BAKER')
1681
              GradTmp(1:IntFroz)=0.d0
1682
           CASE ('CART','HYBRID')
1683
              DO I=1,Nat
1684
                 IF (FrozAtoms(I)) THEN
1685
                    Iat=Order(I)
1686
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1687
                 END IF
1688
              END DO
1689
           CASE DEFAULT
1690
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1691
              STOP
1692
           END SELECT
1693
        END IF ! matches IF (NFroz.NE.0) THEN
1694

    
1695
        IF (debug) THEN
1696
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
1697
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1698
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1699
           ELSE
1700
              DO Iat=1,Nat
1701
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1702
                      XyzGeomF(IGeom,1:3,Iat)
1703
              END DO
1704
           END IF
1705
        END IF
1706

    
1707
        SELECT CASE (COORD)
1708
        CASE ('ZMAT','MIXED','BAKER')
1709
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1710
        CASE ('HYBRID','CART')
1711
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1712
        CASE DEFAULT
1713
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1714
           STOP
1715
        END SELECT
1716

    
1717
        IF (debug) THEN
1718
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1719
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1720
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1721
           ELSE
1722
              DO Iat=1,Nat
1723
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1724
                      XyzGeomF(IGeom,1:3,Iat)
1725
              END DO
1726
           END IF
1727
        END IF
1728
     ELSE ! Optprod
1729
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1730
     END IF !matches IF (OptProd) THEN 
1731

    
1732
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1733

    
1734
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1735
     !
1736
     !  Optimizing other geometries
1737
     !
1738
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1739

    
1740

    
1741

    
1742
     DO IGeom=2,NGeomF-1 ! matches at L1556
1743
        if (debug) WRITE(*,*) '****************************'
1744
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1745
        if (debug) WRITE(*,*) '****************************'
1746

    
1747
        SELECT CASE (COORD)
1748
        CASE ('ZMAT','BAKER','MIXED')
1749
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
1750
           ! PFL 6 Apr 2008 ->
1751
           ! Special case, if FTan=0. we do not consider Tangent at all.
1752
           ! To DO: add the full treatment of FTan
1753
           if (FTan==0.) GradTmp2=ZeroVec
1754
           ! <- PFL 6 Apr 2008
1755
           if (debug) THEN
1756
              Print *, 'L1500, IntTangent(',IGeom,',:)='
1757
              Print *, IntTangent(IGeom,:)
1758
           END IF
1759
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
1760
           SELECT CASE (Step_method)
1761
           CASE ('RFO')
1762
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1763
           CASE ('GDIIS')
1764
              !The following has no effect at IOpt==1
1765
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
1766
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&        
1767
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1768
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
1769
              Print *, GradTmp
1770
           CASE ('GEDIIS')
1771
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1772
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
1773
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1774
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
1775
              Print *, GradTmp
1776
           CASE DEFAULT
1777
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
1778
              STOP
1779
           END SELECT
1780
        CASE ('HYBRID','CART')
1781
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
1782
           ! so we change it
1783
           DO I=1,Nat
1784
              DO J=1,3
1785
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
1786
              END DO
1787
           END DO
1788
           !           GradTmp2=XyzTangent(IGeom,:)
1789
           ! PFL 6 Apr 2008 ->
1790
           ! Special case, if FTan=0. we do not consider Tangent at all.
1791
           ! To DO: add the full treatment of FTan
1792
           if (FTan==0.) GradTmp2=ZeroVec
1793
           ! <- PFL 6 Apr 2008
1794

    
1795
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1796
        CASE DEFAULT
1797
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
1798
           STOP
1799
        END SELECT
1800

    
1801
        ! we take care that frozen atoms do not move
1802
        IF (NFroz.NE.0) THEN
1803
           SELECT CASE (COORD)
1804
           CASE ('ZMAT','MIXED','BAKER')
1805
              IF (debug) THEN 
1806
                 WRITE(*,*) 'Step computed. Coord=',Coord
1807
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
1808
              END IF
1809
              GradTmp(1:IntFroz)=0.d0
1810
           CASE ('CART','HYBRID')
1811
              if (debug) THEN
1812
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
1813
                 DO I=1,Nat
1814
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
1815
                 END DO
1816
              END IF
1817
              DO I=1,Nat
1818
                 IF (FrozAtoms(I)) THEN
1819
                    if (debug) THEN
1820
                       write(*,*) 'Step Computed. Coord=',Coord
1821
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
1822
                    END IF
1823
                    Iat=Order(I)
1824
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1825
                 END IF
1826
              END DO
1827

    
1828
                 if (debug) THEN
1829
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
1830
                    DO I=1,Nat
1831
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
1832
                    END DO
1833
                 END IF
1834

    
1835
           CASE DEFAULT
1836
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
1837
              STOP
1838
           END SELECT
1839
        END IF !matches IF (NFroz.NE.0) THEN
1840

    
1841
        IF (debug) THEN
1842
           SELECT CASE (COORD)
1843
           CASE ('ZMAT','MIXED','BAKER')
1844
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
1845
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1846
           CASE ('CART','HYBRID')
1847
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
1848
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1849
           CASE DEFAULT
1850
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
1851
              STOP
1852
           END SELECT
1853
        END IF ! matches if (debug) THEN
1854

    
1855
        ! we check that the maximum displacement is not greater than Smax
1856
        If (debug) WRITE(*,*) "Fini before test:",Fini
1857
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1858
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
1859
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1860
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
1861

    
1862
        GradTmp=GraDTmp*FactStep
1863

    
1864
        if (debug) WRITE(*,*) "DBG Main, L1475, FactStep=",FactStep
1865
        if (debug.AND.(COORD.EQ.'ZMAT')) WRITE(*,'(A,12(1X,F10.6))') 'Scaled Step:',GradTmp(1:min(12,NFree))
1866

    
1867
        IF (debug) THEN
1868
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
1869
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1870
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1871
           ELSE
1872
              DO Iat=1,Nat
1873
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1874
                      XyzGeomF(IGeom,1:3,Iat)
1875
              END DO
1876
           END IF
1877
        END IF ! matches if (debug) THEN
1878

    
1879
        SELECT CASE (COORD)
1880
        CASE ('ZMAT')
1881
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1882
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1883
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1884
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1885
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1886
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
1887
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1888
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
1889
           Idx=4
1890
           DO Iat=4,Nat
1891
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1892
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
1893
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
1894
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
1895
              Idx=Idx+3
1896
           END DO
1897

    
1898
        CASE ('BAKER')
1899
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1900
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1901
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1902
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
1903

    
1904
        CASE ('MIXED')
1905
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1906
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1907
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1908
           DO Iat=1,NCart
1909
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
1910
           END DO
1911
           Idx=3*NCart+1
1912
           IF (NCart.EQ.1) THEN
1913
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1914
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
1915
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1916
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
1917

    
1918
              Idx=7
1919
           END IF
1920
           IF (NCart.EQ.2) THEN
1921
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1922
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
1923
              Idx=9
1924
           END IF
1925

    
1926
           
1927
           DO Iat=max(NCart+1,4),Nat
1928
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1929
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
1930
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
1931
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
1932
              Idx=Idx+3
1933
           END DO
1934
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1935

    
1936
        CASE ('HYBRID','CART')
1937
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1938
        CASE DEFAULT
1939
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
1940
           STOP
1941
        END SELECT
1942

    
1943
        IF (debug) THEN
1944
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
1945
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1946
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1947
           ELSE
1948
              DO Iat=1,Nat
1949
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1950
                      XyzGeomF(IGeom,1:3,Iat)
1951
              END DO
1952
           END IF
1953
        END IF ! matches if (debug) THEN
1954

    
1955
        ! We project out the tangential part of the gradient to know if we are done
1956
        GradTmp=Grad(IGeom,:)
1957
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
1958
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
1959
        SELECT CASE (COORD)
1960
        CASE ('ZMAT','MIXED','BAKER')
1961
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
1962
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
1963
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
1964
           Ca=S/(sqrt(Norm)*NormGrad)
1965
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
1966
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
1967
        CASE ('CART','HYBRID')
1968
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
1969
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
1970
           Ca=S/(sqrt(Norm)*NormGrad)
1971
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
1972
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
1973
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
1974
        CASE DEFAULT
1975
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
1976
           STOP
1977
        END SELECT
1978
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
1979
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
1980
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1981
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
1982

    
1983
     END DO ! matches DO IGeom=2,NGeomF-1
1984
     ! we save the old gradients
1985
     GradOld=Grad
1986
     EPathOld=Energies
1987

    
1988

    
1989

    
1990
     ! Shall we re-parameterize the path ?
1991
     ! PFL 2007/Apr/10 ->
1992
     ! We call PathCreate in all cases, it will take care of the 
1993
     ! reparameterization, as well as calculating the tangents.
1994
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
1995
     XyzGeomI=XyzGeomF
1996
     IntCoordI=IntCoordF
1997

    
1998
     Call PathCreate(IOpt)
1999
     !     END IF
2000
     ! <- PFL 2007/Apr/10
2001

    
2002
     if (debug) WRITE(*,'(1X,A,I5,A,I5,A,L2)')  'Path.f90, L1415, IOpt=', IOpt, 'MaxCyc=', MaxCyc,' Fini=', Fini
2003
     IF ((.NOT.Fini).AND.(IOpt.LT.MaxCyc)) THEN
2004

    
2005
        ! We have the new path, we calculate its energy and gradient
2006

    
2007
        Call EgradPath(IOpt)
2008
        !IF(IOPT .EQ. 10) Then
2009
        !         Print *, 'Stopping at my checkpoint.'
2010
        !           STOP !This is my temporary checkpoint.
2011
        !ENDIF
2012

    
2013
        ! PFL 31 Mar 2008 ->
2014
        ! Taken from v3.99 but we added a flag...
2015
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
2016
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2017
        ! else it is decreased by 0.8
2018

    
2019
        if (DynMaxStep) THEN
2020
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
2021
           if (OptReac) THEN
2022
              If (Energies(1)<EPathOld(1)) Then
2023
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
2024
              ELSE
2025
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
2026
              END IF
2027
           END IF
2028

    
2029
           if (OptProd) THEN
2030
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
2031
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
2032
              ELSE
2033
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
2034
              END IF
2035
           END IF
2036

    
2037
           DO IGeom=2,NGeomF-1
2038
              If (Energies(IGeom)<EPathOld(IGeom)) Then
2039
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
2040
              ELSE
2041
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
2042
              END IF
2043
           END DO
2044
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2045
        END IF ! test on DynMaxStep
2046
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2047
        ! <- PFL 31 MAr 2008
2048
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
2049
        !  It should get updated for other cases too.
2050

    
2051
        DO IGeom=1,NGeomF
2052
           SELECT CASE (COORD)
2053
           CASE('ZMAT')
2054
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2055
              GeomTmp=IntCoordF(IGeom,:)
2056
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2057
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2058
                   OrderInv(indzmat(2,2)),Geomtmp(1)
2059
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2060
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
2061
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
2062
              Idx=4
2063
              DO Iat=4,Nat
2064
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2065
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
2066
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
2067
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
2068
                 Idx=Idx+3
2069
              END DO
2070
           CASE('BAKER')
2071
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2072
              !GeomTmp=IntCoordF(IGeom,:)
2073
           CASE('CART','HYBRID','MIXED')
2074
              WRITE(IOOUT,'(1X,I5)') Nat
2075
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2076
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2077
              DO I=1,Nat
2078
                 Iat=I
2079
                 If (renum) Iat=OrderInv(I)
2080
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
2081
              END DO
2082
           CASE DEFAULT
2083
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
2084
              STOP
2085
           END SELECT
2086
        END DO ! matches DO IGeom=1,NGeomF
2087

    
2088
        Call Write_path(IOpt)
2089
! Shall we analyze the geometries ?
2090
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
2091

    
2092
        ! We update the Hessian matrices
2093
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
2094
        ! First using the displacement from the old path
2095
        IGeom0=2
2096
        IGeomF=NGeomF-1
2097
        IF (OptReac) IGeom0=1
2098
        IF (OptProd) IGeomF=NGeomF
2099

    
2100
        IF (FCalcHess) THEN
2101
           DO IGeom=IGeom0,IGeomF
2102
              SELECT CASE (COORD)
2103
              CASE ('ZMAT','MIXED','BAKER')
2104
                 GeomTmp2=IntCoordF(IGeom,:)
2105
              CASE ('CART','HYBRID')
2106
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2107
              CASE DEFAULT
2108
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2109
                 STOP
2110
              END SELECT
2111
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
2112
           END DO
2113
        ELSE
2114
           DO IGeom=IGeom0,IGeomF
2115
              GeomTmp=GeomOld_all(IGeom,:)
2116
              SELECT CASE (COORD)
2117
              CASE ('ZMAT','MIXED','BAKER')
2118
                 GeomTmp2=IntCoordF(IGeom,:)
2119
              CASE ('CART','HYBRID')
2120
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2121
              CASE DEFAULT
2122
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2123
                 STOP
2124
              END SELECT
2125
              GeomTmp=GeomTmp2-GeomTmp
2126
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
2127
              HessTmp=Hess(IGeom,:,:)
2128
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2129
              Hess(IGeom,:,:)=HessTmp
2130
           END DO ! matches DO IGeom=IGeom0,IGeomF
2131

    
2132
           ! Second using the neighbour points
2133
           IF (HupNeighbour) THEN
2134
              ! Only one point for end points.
2135
              IF (OptReac) THEN
2136
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2137
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
2138
                 ELSE
2139
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
2140
                 END IF
2141
                 GradTmp=Grad(1,:)-Grad(2,:)
2142
                 HessTmp=Hess(1,:,:)
2143
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2144
                 Hess(1,:,:)=HessTmp
2145
              END IF
2146

    
2147
              IF (OptProd) THEN
2148
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2149
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
2150
                 ELSE
2151
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
2152
                 END IF
2153
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
2154
                 HessTmp=Hess(NGeomF,:,:)
2155
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2156
                 Hess(NGeomF,:,:)=HessTmp
2157
              END IF
2158

    
2159
              ! Two points for the rest of the path
2160
              DO IGeom=2,NGeomF-1
2161
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
2162
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
2163
                 ELSE
2164
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
2165
                 END IF
2166
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
2167
                 HessTmp=Hess(IGeom,:,:)
2168
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2169

    
2170
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2171
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
2172
                 ELSE
2173
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
2174
                 END IF
2175
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
2176

    
2177
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2178
                 Hess(IGeom,:,:)=HessTmp
2179
              END DO ! matches DO IGeom=2,NGeomF-1
2180
           END IF !matches IF (HupNeighbour) THEN
2181
        END IF ! matches IF (FCalcHess)
2182
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
2183
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
2184

    
2185
!   IF (PROG=="GAUSSIAN") THEN
2186
!      DEALLOCATE(Gauss_paste)
2187
!      DEALLOCATE(Gauss_root)
2188
!      DEALLOCATE(Gauss_comment)
2189
!      DEALLOCATE(Gauss_end)
2190
!   END IF
2191

    
2192
  DEALLOCATE(XyzGeomF, IntCoordF)
2193
  DEALLOCATE(XyzGeomI, IntCoordI)
2194
  DEALLOCATE(XyzTangent,IntTangent)
2195
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2196
  DEALLOCATE(GeomOld_all)
2197

    
2198
  if (PROG=="GAUSSIAN") THEN
2199
     ! We de-allocate the Gauss_XX lists
2200
     ! transverse the list and deallocate each node
2201
     previous => Gauss_root ! make current point to head of list
2202
     DO
2203
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2204
        current => previous%next ! make list point to next node of head
2205
        DEALLOCATE(previous) ! deallocate current head node
2206
        previous => current  ! make current point to new head
2207
     END DO
2208

    
2209
     previous => Gauss_comment ! make current point to head of list
2210
     DO
2211
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2212
        current => previous%next ! make list point to next node of head
2213
        DEALLOCATE(previous) ! deallocate current head node
2214
        previous => current  ! make current point to new head
2215
     END DO
2216

    
2217
     previous => Gauss_end ! make current point to head of list
2218
     DO
2219
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2220
        current => previous%next ! make list point to next node of head
2221
        DEALLOCATE(previous) ! deallocate current head node
2222
        previous => current  ! make current point to new head
2223
     END DO
2224

    
2225

    
2226
  END IF
2227

    
2228
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2229
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2230
  IF (COORD.EQ.'BAKER') THEN
2231
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2232
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2233
  END IF
2234

    
2235
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2236

    
2237
  Close(IOIN)
2238
  Close(IOOUT)
2239
  IF (AnaGeom) Close(IODAT)
2240
  IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot)
2241

    
2242
END PROGRAM PathOpt