Statistics
| Revision:

root / src / Path.f90 @ 10

History | View | Annotate | Download (89.9 kB)

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

    
651
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
652
     IOOUT=6
653
  CASE (0)
654
     IOIN=5
655
     IOOUT=6
656
  END SELECT
657

    
658

    
659
  ! We initiliaze variables
660
  Pi=dacos(-1.0d0)
661
  Nat=1
662
  Fact=1.1
663
  IGeomRef=-1
664
  NGeomI=1
665
  NGeomF=8
666
  IReparam=5
667
  IReparamT=-1
668
  ISpline=5
669
  debug=.False.
670
  MW=.TRUE.
671
  bohr=.False.
672
  Coord='ZMAT'
673
  Step_method='RFO'
674
  DebugFile='Path.valid'
675
  PathName="Path"
676
  MaxCyc=50
677
  SMax=0.1D0
678
  SThresh=-1.
679
  GThresh=-1.
680
  FTan=0.
681
  Hinv=.TRUE.
682
  IniHUp=.FALSE.
683
  HupNeighbour=.FALSE.
684
  HesUpd="EMPTY"
685
  HUpdate="EMPTY"
686
  FFrozen=.FALSE.
687
  FCart=.FALSE.
688
  List=0
689
  renum=.TRUE.
690
  OptReac=.FALSE.
691
  OptProd=.FALSE.
692
  CalcEReac=.FALSE.
693
  CalcEProd=.FALSE.
694
  EReac=0.d0
695
  EProd=0.d0 
696
  OptGeom=-1
697
  GeomFile="EMPTY"
698
  AnaGeom=.FALSE.
699
  AnaFile="EMPTY"
700
  GplotFile="EMPTY"
701
  Nb=0
702
  NbCom=0
703
  PathOnly=.FALSE.
704
  Prog='EMPTY'
705
  ProgExe='EMPTY'
706
  Runmode='SERIAL'
707
  Input='EMPTY'
708
  Poscar="POSCAR"
709
  AutoCart=.FALSE.
710
  VMD=.TRUE.
711
  WriteVASP=.FALSE.
712
  DynMaxStep=.FALSE.
713
  CalcName="EMPTY"
714
  ISuffix="EMPTY"
715
  OSuffix="EMPTY"
716
  NGintMax=10
717
  Align=.TRUE.
718

    
719
  ! Number of point used to compute the length of the path,
720
  ! and to reparameterize the path
721
  NMaxPtPath=-1 
722
  FrozTol=1e-4
723

    
724
  ! Read the namelist
725
  read (IOIN,path)
726

    
727
  debug=valid("Main").or.valid("Path")
728

    
729
  FCalcHess=valid("CalcHess")
730
  PathName=AdjustL(PathName)
731
  Coord=AdjustL(Coord)
732
  CALL UpCase(coord)
733

    
734
  Step_method=AdjustL(Step_method)
735
  CALL UpCase(Step_method)
736

    
737
  Prog=AdjustL(Prog)
738
  CALL UpCase(prog)
739

    
740
  Input=AdjustL(Input)
741
  CALL UpCase(input)
742

    
743
  Runmode=AdjustL(Runmode)
744
  Call UpCase(Runmode)
745
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"
746

    
747
  if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN
748
     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"
749
     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"
750
     CalcEReac=.FALSE.
751
     CalcEProd=.FALSE.
752
  END IF
753

    
754
  If (IReparamT<=0) IReparamT=IReparam
755

    
756
! We take care of HesUpd flag
757
! this is needed because some users might still use it
758
  if (HesUpd.NE."EMPTY") THEN
759
! We are pityless:
760
     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"
761
     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"
762
     STOP
763
  END IF
764

    
765
  IF (HUpdate=='EMPTY') THEN
766
     IF (OptGeom>=1) THEN
767
        HupDate="BFGS"
768
     ELSE
769
        HUpdate="MS"
770
     END IF
771
  END IF
772

    
773
  If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 
774
     COORD='CART'
775
  END IF
776

    
777
  IF (COORD.EQ.'CART') THEN
778
     Renum=.FALSE.
779
     IF (OptGeom>0) THEN
780
        IGeomRef=OptGeom
781
     ELSE
782
        IGeomRef=1
783
     END IF
784
  END IF
785

    
786
  
787
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
788
  if (Prog.EQ.'2D') Prog="TEST2D"
789
  if (Prog.EQ.'TEST3') Prog="LEPS"
790

    
791

    
792
  Select Case (Prog)
793
    CASE ("EMPTY","G09","GAUSSIAN")
794
     Prog='GAUSSIAN'
795
     If (ProgExe=="EMPTY") ProgExe="g09"
796
     if (CalcName=="EMPTY") CalcName="Path"
797
     if (ISuffix=="EMPTY")  ISuffix="com"
798
     if (OSuffix=="EMPTY")  OSuffix="log"
799
     UnitProg="au"
800
     ConvE=au2kcal
801
    CASE ("G03")
802
     Prog='GAUSSIAN'
803
     If (ProgExe=="EMPTY") ProgExe="g03"
804
     if (CalcName=="EMPTY") CalcName="Path"
805
     if (ISuffix=="EMPTY")  ISuffix="com"
806
     if (OSuffix=="EMPTY")  OSuffix="log"
807
     UnitProg="au"
808
     ConvE=au2kcal
809
    CASE ("EXT")
810
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
811
     if (CalcName=="EMPTY") CalcName="Ext"
812
     if (ISuffix=="EMPTY")  ISuffix="in"
813
     if (OSuffix=="EMPTY")  OSuffix="out"
814
     UnitProg="au"
815
     ConvE=au2kcal
816
    CASE ('MOPAC')
817
     If (ProgExe=="EMPTY") ProgExe="mopac"
818
     if (CalcName=="EMPTY") CalcName="Path"
819
     if (ISuffix=="EMPTY")  ISuffix="mop"
820
     if (OSuffix=="EMPTY")  OSuffix="out"
821
     UnitProg="au"
822
     ConvE=au2kcal
823
    CASE ("SIESTA")
824
     If (ProgExe=="EMPTY") ProgExe="siesta"
825
     if (CalcName=="EMPTY") CalcName="Path"
826
     if (ISuffix=="EMPTY")  ISuffix="fdf"
827
     if (OSuffix=="EMPTY")  OSuffix="out"
828
     UnitProg="eV"
829
     ConvE=eV2kcal
830
    CASE ("VASP")
831
     If (ProgExe=="EMPTY") ProgExe="VASP"
832
     UnitProg="eV"
833
     ConvE=eV2kcal
834
     FPBC=.TRUE.
835
     IPer=3
836
     Kabeg=-1
837
     kaEnd=1
838
     kbBeg=-1
839
     kbEnd=1
840
     kcBeg=-1
841
     kcEnd=1
842
    CASE ("TM","TURBOMOLE")
843
     Prog="TURBOMOLE"
844
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
845
     UnitProg="au"
846
     ConvE=au2kcal
847
     FPBC=.FALSE.
848
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
849
     ProgExe=""
850
     UnitProg="au"
851
     ConvE=au2kcal
852
    CASE DEFAULT
853
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
854
     STOP
855
  END SELECT 
856

    
857
  if (Input.EQ.'EMPTY') THEN
858
     Select Case (Prog)
859
       CASE ("VASP")
860
          Input=Prog
861
       CASE ("CHAMFRE")
862
          Input="VASP"
863
       CASE DEFAULT
864
          Input='XYZ'
865
     END SELECT
866
    WRITE(*,*) "Input has been set to the default: ",INPUT
867
 END IF
868

    
869
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
870

    
871
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
872
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
873
     STOP
874
  END IF
875

    
876
  if (OptGeom.GE.1) THEN
877
     FPrintGeom=.FALSE.
878
     If (GeomFile/='EMPTY') THEN
879
        FPrintGeom=.TRUE.
880
     END IF
881
     IF (IGeomRef.LE.0) THEN
882
        IGeomRef=OptGeom
883
     ELSE
884
        IF (IGeomRef.NE.OptGeom) THEN
885
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
886
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
887

    
888
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
889
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
890
           STOP
891
        END IF
892
     END IF
893
  END IF
894

    
895
  ALLOCATE(Cart(Nat))
896
  Cart=0
897

    
898
  ALLOCATE(Frozen(Nat))
899
  Frozen=0
900

    
901
  IF (FCart) THEN
902
! We rewind the file to be sure to browse all of it to read the Cart namelist
903
     REWIND(IOIN)
904
     List=0
905
     READ(IOIN,cartlist)
906
     IF (COORD.NE.'CART') THEN
907
        Cart=List(1:Nat)
908
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
909
     ELSE
910
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
911
        WRITE(*,*) "I will discard the cartlist"
912
        Cart=0
913
     END IF
914
  END IF
915

    
916
  if (FFrozen) THEN
917

    
918
     REWIND(IOIN)
919
     List=0
920
     READ(IOIN,Frozenlist)
921
     Frozen=List(1:Nat)
922
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
923
  END IF
924

    
925
  If (FCart) Call Expand(Cart,NCart,Nat)
926
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
927

    
928
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
929
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
930

    
931
  if (AnaGeom) THEN
932
     REWIND(IOIN)
933
     READ(IOIN,AnaList,IOSTAT=IoS)
934
     IF (IOS==0) THEN ! we have a AnaList namelist 
935
        Call ReadAnaList
936
     ELSE
937
! No AnaList namelist, we will print only the energy
938
        FormAna='(1X,F8.3,1X,1X,F7.2,1X,F15.6)'
939
     END IF
940
     IF (AnaFile=="EMPTY") THEN
941
        AnaFile=TRIM(PathName) // '.dat'
942
     END IF
943
     OPEN(IODat,File=AnaFile)
944
     IF (GplotFile=="EMPTY") THEN
945
        GplotFile=TRIM(PathName) // '.gplot'
946
     END IF
947
     
948
     OPEN(IODat,File=AnaFile)
949
     If (OptGeom<=0) THEN
950
        OPEN(IOGplot,File=GPlotFile)
951
        WRITE(IoGplot,'(A)') "#!/usr/bin/gnuplot -persist"
952
        WRITE(IoGplot,'(A)') " set pointsize 2"
953
        WRITE(IoGplot,'(A)') " xcol=1"
954
        WRITE(IoGplot,'(A,I5)') " Ecol=",NbVar+3
955
     END IF
956
  END IF
957

    
958
        
959

    
960
  if (NMaxPtPath.EQ.-1) then
961
     NMaxPtPath=min(50*NGeomF,2000)
962
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
963
  end if
964

    
965
  !  NMaxPtPath=10000
966

    
967
  ! We ensure that FTan is in [0.:1.]
968
  FTan=min(abs(FTan),1.d0)
969

    
970
! PFL 2011 Mar 14 ->
971
! Added some printing for analyses with Anapath
972
  if (debug) THEN
973
     WRITE(IOOUT,path)
974
  ELSE
975
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC
976
! and PATHNAME 
977
     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat
978
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI
979
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF
980
     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC
981
     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)
982
  END IF
983

    
984
  FTan=1.-FTan
985

    
986
  !Thresholds...
987
  if (SThresh.LE.0) SThresh=0.01d0
988
  If (GThresh.LE.0) GThresh=SThresh/4.
989
  SThreshRms=Sthresh*2./3.
990
  GThreshRms=Gthresh*2./3.
991

    
992
  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
993
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
994
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
995
  ALLOCATE(MassAt(NAt),OrderInv(Nat))
996

    
997
  AtName=""
998
  MassAt=1.
999

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

    
1004
  REWIND(IOIN)
1005
! We search the '&path' namelist
1006
  READ(IOIN,'(A)',Iostat=Ios) Line
1007
  Fini=.FALSE.
1008
  DO WHILE ((Ios==0).AND.InString(Line,"&PATH")==0)
1009
     if (debug) WRITE(*,*) "Lookgin for Path:",TRIM(Line)
1010
     READ(IOIN,'(A)',Iostat=Ios) Line
1011
  END DO
1012
  if (debug) WRITE(*,*) "Path found:",TRIM(LINE)
1013
  ! We are on the &path line, we move to the end
1014
  Call NoComment(Line)
1015
  DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
1016
     if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
1017
     READ(IOIN,'(A)',Iostat=Ios) Line
1018
     Call NoComment(Line)
1019
     Call noString(Line)
1020
  END DO
1021
! We are at then end of &Path / namelist
1022
! We scan for other Namelists
1023
  READ(IOIN,'(A)',Iostat=Ios) Line
1024
  if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
1025
  DO WHILE  ((IoS==0).AND.(Index(Line,'&')/=0))
1026
     Call NoComment(Line)
1027
     if (debug) WRITe(*,*) "Another Namelist ?",TRIM(Line)
1028
     Idx=InString(Line,'Analist')
1029
     DO WHILE ((Ios==0).AND.InString(Line,"/")==0)
1030
        if (debug) WRITE(*,*) "Lookgin for /:",TRIM(Line)
1031
        READ(IOIN,'(A)',Iostat=Ios) Line
1032
        Call NoComment(Line)
1033
     END DO
1034
     If (Idx/=0) THEN
1035
! we have just read &Analist namelist, we have to skip the variables description
1036
        DO I=1,Nb
1037
           READ(IOIN,'(A)',Iostat=Ios) Line
1038
           if (debug) WRITe(*,*) "We skip Analist",TRIM(LINE)
1039
        END DO
1040
     END IF
1041
     READ(IOIN,'(A)',Iostat=Ios) Line
1042
     if (debug) WRITe(*,*) "We skip the line:",TRIM(LINE)
1043
  END DO
1044
  BACKSPACE(IOIN)
1045

    
1046
  ! We read the initial geometries
1047
  Call Read_Geom(input)
1048

    
1049
  ! We convert atome names into atom mass number
1050
  DO I=1,NAt
1051
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
1052
     Atome(I)=ConvertNumAt(AtName(I))
1053
  END DO
1054

    
1055
 If (AnaGeom) THEN
1056
    If (NbVar>0) THEN
1057
! We print the description of the varialbes in AnaFile
1058
       Call PrintAnaList(IoDat)
1059
       If (OptGeom<=0) THEN
1060
          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
1061
       ELSE
1062
          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
1063
       END IF
1064
    ELSE
1065
       If (OptGeom<=0) THEN
1066
          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
1067
       ELSE
1068
          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"       
1069
       END IF
1070
    END IF
1071
 END IF
1072

    
1073
  ! PFL 23 Sep 2008 ->
1074
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
1075
  ! between the first geometry and the others.
1076
  IF (NFroz.GT.0) THEN
1077
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
1078
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
1079
     ListAtFroz=.FALSE.
1080

    
1081
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
1082
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
1083
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
1084

    
1085
     SELECT CASE (NFroz)
1086
        ! We have Frozen atoms
1087
        ! PFL 28 Oct 2008 ->
1088
        ! It might happen that the geometries are translated/rotated.
1089
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
1090
        ! three frozen atoms.
1091

    
1092
     CASE (3:)
1093
        DO I=1,NFroz
1094
           ListAtFroz(Frozen(I))=.TRUE.
1095
        END DO
1096
     CASE (2)
1097
        IAt1=Frozen(1)
1098
        IAt2=Frozen(2)
1099
        ListAtFroz(Iat1)=.TRUE.
1100
        ListAtFroz(Iat2)=.TRUE.
1101
        x2(1)=x1(Iat2)-x1(Iat1)
1102
        y2(1)=y1(Iat2)-y1(Iat1)
1103
        z2(1)=z1(Iat2)-z1(Iat1)
1104
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1105
        Dist=-1.
1106
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1107
        DO I=1,Nat
1108
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
1109
           x2(2)=x1(Iat2)-x1(I)
1110
           y2(2)=y1(Iat2)-y1(I)
1111
           z2(2)=z1(Iat2)-z1(I)
1112
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
1113
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
1114
           if (abs(ca)<0.9) THEN
1115
              IF (Norm2>Dist) THEN
1116
                 Iat3=I
1117
                 Dist=Norm2
1118
              END IF
1119
           END IF
1120
        END DO
1121
        ListAtFroz(Iat3)=.TRUE.
1122
        if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3
1123
     CASE (1)
1124
        IAt1=Frozen(1)
1125
        ListAtFroz(Iat1)=.TRUE.
1126
        Dist=-1.
1127
        ! We scan all atoms to find one that is further from At1
1128
        DO I=1,Nat
1129
           if (I.EQ.Iat1) Cycle
1130
           x2(1)=x1(Iat1)-x1(I)
1131
           y2(1)=y1(Iat1)-y1(I)
1132
           z2(1)=z1(Iat1)-z1(I)
1133
           Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1134
           IF (Norm1>Dist) THEN
1135
              Iat2=I
1136
              Dist=Norm1
1137
           END IF
1138
        END DO
1139
        ListAtFroz(Iat2)=.TRUE.
1140
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
1141
        x2(1)=x1(Iat2)-x1(Iat1)
1142
        y2(1)=y1(Iat2)-y1(Iat1)
1143
        z2(1)=z1(Iat2)-z1(Iat1)
1144
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
1145
        Dist=-1.
1146
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
1147
        Iat3=1
1148
        DO WHILE (ListAtFroz(Iat3).AND.(Iat3<=Nat))
1149
           Iat3=Iat3+1
1150
        END DO
1151
        DO I=1,Nat
1152
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
1153
           x2(2)=x1(Iat2)-x1(I)
1154
           y2(2)=y1(Iat2)-y1(I)
1155
           z2(2)=z1(Iat2)-z1(I)
1156
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
1157
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
1158
           if (abs(ca)<0.9) THEN
1159
              IF (Norm2>Dist) THEN
1160
                 Iat3=I
1161
                 Dist=Norm2
1162
              END IF
1163
           END IF
1164
        END DO
1165
        ListAtFroz(Iat3)=.TRUE.
1166
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3
1167
     END SELECT
1168

    
1169
     DO I=2,NGeomI
1170
        ! First, we align the Ith geometry on the first one, using only
1171
        ! the positions of the "frozen" atoms. (see above: it might be that
1172
        ! ListAtFroz contains non frozen atoms used to align the geometries
1173
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
1174
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
1175
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
1176
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
1177

    
1178
        FChkFrozen=.FALSE.
1179
        DO J=1,NFroz
1180
           IAt=Frozen(J)
1181
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
1182
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
1183
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
1184
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
1185
              FChkFrozen=.TRUE.
1186
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
1187
              WRITE(*,*) "from geometry 1 to geometry ",I
1188
           END IF
1189
        END DO
1190
     END DO
1191
     IF (FChkFrozen) THEN
1192
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
1193
        STOP
1194
     END IF
1195
  END IF
1196

    
1197

    
1198
  N3at=3*Nat
1199
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
1200

    
1201

    
1202
  Call ReadInput
1203

    
1204
  if (FPBC)  THEN
1205
     Allocate(XGeomRefPBC(Nat),YGeomRefPBC(Nat),ZGeomRefPBC(Nat))
1206
     XGeomRefPBC(1:Nat)=XyzGeomI(1,1,1:Nat)
1207
     YGeomRefPBC(1:Nat)=XyzGeomI(1,2,1:Nat)
1208
     ZGeomRefPBC(1:Nat)=XyzGeomI(1,3,1:Nat)
1209

    
1210
     Call CheckPeriodicBound
1211
  END IF
1212

    
1213
  IF (COORD=="MIXED") Call TestCart(AutoCart)
1214

    
1215
  SELECT CASE(NFroz)
1216
  CASE (2)
1217
     IntFroz=1
1218
  CASE (3)
1219
     IntFroz=3
1220
  CASE (4:)
1221
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
1222
  END SELECT
1223

    
1224
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
1225
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
1226
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
1227

    
1228
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
1229
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
1230
     COORD="ZMAT"
1231
  END IF
1232

    
1233
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
1234

    
1235
  SELECT CASE(COORD)
1236
  CASE ('ZMAT')
1237
     NCoord=Nfree
1238
     ALLOCATE(dzdc(3,nat,3,nat))
1239
  CASE ('MIXED')
1240
     SELECT CASE (NCart+NFroz)
1241
     CASE (1) 
1242
        NCoord=N3At-3
1243
     CASE (2)
1244
        NCoord=N3At-1
1245
     CASE (3:)
1246
        NCoord=N3At
1247
     END SELECT
1248
     ALLOCATE(dzdc(3,nat,3,nat))
1249
  CASE ('BAKER')
1250
     Symmetry_elimination=0
1251
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
1252
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
1253
     ALLOCATE(BTransInv(NCoord,3*Nat))
1254
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
1255
  CASE ('HYBRID')
1256
     NCoord=N3at
1257
  CASE ('CART')
1258
     NCoord=N3at
1259
  END SELECT
1260

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

    
1263
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
1264
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
1265
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
1266
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
1267
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
1268
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
1269
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
1270
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
1271

    
1272
  ZeroVec=0.d0
1273
  DO I =1, NGeomF
1274
     IntTangent(I,:)=0.d0
1275
  END DO
1276

    
1277
  if (debug) THEN
1278
     Print *, 'L886, IntTangent(1,:)='
1279
     Print *, IntTangent(1,:)
1280
  END IF
1281

    
1282
  if (.NOT.OptReac) Energies(1)=Ereac
1283
  if (.NOT.OptProd) Energies(NGeomF)=EProd
1284
  MaxStep=SMax
1285

    
1286
  ! We initialize the mass array
1287
  if (MW) THEN
1288
     WRITE(*,*) 'Working in MW coordinates'
1289
     DO I=1,Nat
1290
        massAt(I)=Mass(Atome(I))
1291
     END DO
1292
  ELSE
1293
     DO I=1,Nat
1294
        MassAt(I)=1.d0
1295
     END DO
1296
  END IF
1297

    
1298
  WRITE(*,*) "Prog=",TRIM(PROG)
1299

    
1300
  ALLOCATE(FrozAtoms(Nat))
1301

    
1302
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
1303
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
1304
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
1305
  IF (NFroz.EQ.0) THEN
1306
     FrozAtoms=.TRUE.
1307
  ELSE
1308
     I=1
1309
     NFr=0
1310
     FrozAtoms=.False.
1311
     DO WHILE (Frozen(I).NE.0)
1312
        IF (Frozen(I).LT.0) THEN
1313
           DO J=Frozen(I-1),abs(Frozen(I))
1314
              IF (.NOT.FrozAtoms(J)) THEN
1315
                 NFr=NFr+1
1316
                 FrozAtoms(J)=.TRUE.
1317
              END IF
1318
           END DO
1319
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1320
           FrozAtoms(Frozen(I))=.TRUE.
1321
           NFr=NFr+1
1322
        END IF
1323
        I=I+1
1324
     END DO
1325
     IF (NFr.NE.NFroz) THEN
1326
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1327
        STOP
1328
     END IF
1329
  END IF
1330

    
1331
!  if (debug) THEN
1332
     !Print *, 'L968, IntTangent(1,:)='
1333
     !Print *, IntTangent(1,:)
1334
!  END IF
1335

    
1336
  ! We have read everything we needed... time to work :-)
1337
  IOpt=0
1338
  FirstTimePathCreate = .TRUE. ! Initialized.
1339
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1340
  FirstTimePathCreate = .FALSE.
1341

    
1342
  if (debug) THEN
1343
     Print *, 'L980, IntTangent(1,:)='
1344
     Print *, IntTangent(1,:)
1345
  END IF
1346

    
1347
  ! PFL 30 Mar 2008 ->
1348
  ! The following is now allocated in Calc_Baker. This is more logical to me
1349
  !   IF (COORD .EQ. "BAKER") Then
1350
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1351
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1352
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1353
  !   END IF
1354
  ! <- PFL 30 mar 2008
1355

    
1356
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1357
  !
1358
  !
1359
  ! OptGeom is the index of the geometry which is to be optimized.
1360
  IF (Optgeom.GT.0) THEN
1361
     Flag_Opt_Geom=.TRUE.
1362
     SELECT CASE(Coord)
1363
     CASE ('CART','HYBRID')
1364
!!! CAUTION : PFL 29.AUG.2008 ->
1365
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1366
        ! This should be made  consistent !!!!!!!!!!!!!!!
1367
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1368
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1369
!!! <- CAUTION : PFL 29.AUG.2008
1370
     CASE ('ZMAT','MIXED')
1371
        !GeomTmp=IntCoordF(OptGeom,:)
1372
        GeomTmp=IntCoordI(OptGeom,:)
1373
     CASE ('BAKER')
1374
        !GeomTmp=IntCoordF(OptGeom,:)
1375
        GeomTmp=IntCoordI(OptGeom,:)
1376
     CASE DEFAULT
1377
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1378
        STOP
1379
     END SELECT
1380

    
1381
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1382
     Flag_Opt_Geom = .TRUE.
1383
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1384

    
1385
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1386
     if (debug) WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1387
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1388
     STOP
1389
  END IF ! IF (Optgeom.GT.0) THEN
1390

    
1391
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392
  ! End of GeomOpt
1393
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1394

    
1395
  IF (PathOnly) THEN
1396
     Call Header("PathOnly=.T. , Stop here")
1397
     Call Write_path(-1)
1398
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1399
     STOP
1400
  END IF
1401

    
1402
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1403

    
1404
  DEALLOCATE(XyzGeomI,IntCoordI)
1405
  NGeomI=NGeomF
1406
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1407

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

    
1413
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1414
     DO IGeom=1,NGeomF
1415
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1416
        DO I=1,3*Nat-6
1417
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1418
        END DO
1419
     END DO
1420
  END IF
1421

    
1422
  ! To calculate B^T^-1 for all extrapolated final geometries:
1423
  ! PFL 18 Jan 2008 ->
1424
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1425
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1426
  ! <-- PFL 18 Jan 2008
1427

    
1428
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1429
     DO IGeom=1,NGeomF
1430
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
1431
        DO I=1,3*Nat-6
1432
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1433
        END DO
1434
     END DO
1435
  END IF
1436

    
1437
  IOpt=0
1438
  Call EgradPath(IOpt)
1439

    
1440
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1441
  Call Write_path(IOpt)
1442
! Shall we analyze the geometries ?
1443
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1444

    
1445
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1446

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

    
1451
  ! By default, Hess=Id
1452
  Hess=0.
1453
  IF(Coord .EQ. "ZMAT") Then
1454
     ! We use the fact that we know that approximate good hessian values
1455
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1456
     DO IGeom=1,NGeomF
1457
        Hess(IGeom,1,1)=0.5d0
1458
        Hess(IGeom,2,2)=0.5d0
1459
        Hess(IGeom,3,3)=0.1d0
1460
        DO J=1,Nat-3
1461
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1462
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1463
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1464
        END DO
1465
        IF (HInv) THEN
1466
           DO I=1,NCoord
1467
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1468
           END DO
1469
        END IF
1470
     END DO
1471
  ELSE
1472
     DO I=1,NCoord
1473
        DO J=1,NGeomF
1474
           Hess(J,I,I)=1.d0
1475
        END DO
1476
     END DO
1477
  END IF
1478

    
1479
  IF (COORD .EQ. "BAKER") THEN
1480
     ! UMat(NPrim,NCoord)
1481
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1482
     Hprim=0.d0
1483
     ScanCoord=>Coordinate
1484
     I=0
1485
     DO WHILE (Associated(ScanCoord%next))
1486
        I=I+1
1487
        SELECT CASE (ScanCoord%Type)
1488
        CASE ('BOND')
1489
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1490
           Hprim(I,I) = 0.5d0
1491
           !END DO
1492
        CASE ('ANGLE')
1493
           Hprim(I,I) = 0.2d0
1494
        CASE ('DIHEDRAL')
1495
           Hprim(I,I) = 0.1d0
1496
        END SELECT
1497
        ScanCoord => ScanCoord%next
1498
     END DO
1499

    
1500
     ! HprimU = H^prim U:
1501
     HprimU=0.d0
1502
     DO I=1, 3*Nat-6
1503
        DO J=1,NPrim
1504
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1505
        END DO
1506
     END DO
1507

    
1508
     ! Hess = U^T HprimU = U^T H^prim U:
1509
     Hess=0.d0
1510
     DO I=1, 3*Nat-6
1511
        DO J=1,NPrim
1512
           ! UMat^T is needed below.
1513
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1514
        END DO
1515
     END DO
1516
     DO K=2,NGeomF
1517
        Hess(K,:,:)=Hess(1,:,:)
1518
     END DO
1519
     DEALLOCATE(Hprim,HprimU)
1520
  end if !  ends IF COORD=BAKER
1521
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1522
  ALLOCATE(HessTmp(NCoord,NCoord))
1523

    
1524
  IF (IniHUp) THEN
1525
     IF (FCalcHess) THEN
1526
        ! The numerical calculation of the Hessian has been switched on
1527
        DO IGeom=1,NGeomF
1528
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1529
              GeomTmp=IntCoordF(IGeom,:)
1530
           ELSE
1531
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1532
           END IF
1533
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1534
        END DO
1535
     ELSE
1536
        ! We generate a forward hessian and a backward one and we mix them.
1537
        ! First, the forward way:
1538
        DO IGeom=2,NGeomF-1
1539
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1540
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1541
           ELSE
1542
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1543
           END IF
1544
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1545
           HessTmp=Hess(IGeom-1,:,:)
1546
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1547
           Hess(IGeom,:,:)=HessTmp
1548
        END DO
1549

    
1550
        ! Now, the backward way:
1551
        HessTmp=Hess(NGeomF,:,:)
1552
        DO IGeom=NGeomF-1,2,-1
1553
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1554
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1555
           ELSE
1556
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1557
           END IF
1558
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1559
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1560
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1561
           ca=cos(alpha)
1562
           sa=sin(alpha)
1563
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1564
        END DO
1565
     END IF ! matches IF (FCalcHess)
1566
  END IF ! matches IF (IniHUp) THEN
1567

    
1568
  ! For debug purposes, we diagonalize the Hessian matrices
1569
  IF (debug) THEN
1570
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1571
     DO I=1,NGeomF
1572
        WRITE(*,*) "DBG Main - Point ",I
1573
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1574
        DO J=1,NCoord
1575
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1576
        END DO
1577
     END DO
1578
  END IF ! matches IF (debug) THEN
1579

    
1580
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1581
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1582
  !
1583
  ! Main LOOP : Optimization of the Path
1584
  !
1585
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1586
  if (debug)   Print *, 'NGeomF=', NGeomF
1587
  allocation_flag = .TRUE.
1588

    
1589
  Fini=.FALSE.
1590

    
1591
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1592
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1593
     IOpt=IOpt+1
1594

    
1595
     SELECT CASE (COORD)
1596
     CASE ('ZMAT','MIXED','BAKER')
1597
        GeomOld_all=IntCoordF
1598
     CASE ('CART','HYBRID')
1599
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1600
     CASE DEFAULT
1601
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1602
        STOP
1603
     END SELECT
1604

    
1605
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1606
          valid("printtangent")) THEN
1607
        WRITE(*,*) "Visualization of tangents"
1608
        Idx=min(12,NCoord)
1609
        DO I=1,NGeomF
1610
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1611
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1612
           WRITE(*,*) 
1613
        END DO
1614
        WRITE(*,*) "END of tangents"
1615
     END IF
1616

    
1617
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1618
          valid("printgrad")) THEN
1619
        WRITE(*,*) "Visualization of gradients"
1620
        Idx=min(12,NCoord)
1621
        DO I=1,NGeomF
1622
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1623
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1624
           WRITe(*,*) 
1625
        END DO
1626
        WRITE(*,*) "END of gradients"
1627
     END IF
1628

    
1629
     Fini=.TRUE.
1630
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1631
        IGeom=1
1632
        if (debug) WRITE(*,*) '**************************************'
1633
        if (debug) WRITE(*,*) 'Optimizing reactant'
1634
        if (debug) WRITE(*,*) '**************************************'
1635
        SELECT CASE (COORD)
1636
        CASE ('ZMAT','MIXED','BAKER')
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,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1641
              Print *, GradTmp
1642
           CASE ('GDIIS')
1643
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1644
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1645
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1646
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1647
              Print *, GradTmp
1648
           CASE ('GEDIIS')
1649
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1650
              ! Energies are updated in EgradPath.f90
1651
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1652
                   NCoord,allocation_flag,ZeroVec)
1653
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1654
              !GradTmp = GeomTmp - IntCoordF(1,:)
1655
              !Print *, 'Old Geometry:'
1656
              !Print *, IntCoordF(1,:)
1657
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1658
              Print *, GradTmp
1659
              !Print *, 'GeomTmp:'
1660
              !Print *, GeomTmp
1661
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1662
           CASE DEFAULT
1663
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1664
              STOP
1665
           END SELECT
1666
        CASE ('HYBRID','CART')
1667
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1668
        CASE DEFAULT
1669
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1670
           STOP
1671
        END SELECT
1672

    
1673
        IF (debug) THEN
1674
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1675
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1676
        END IF
1677

    
1678
        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
1679
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1680
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
1681
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1682
        OptReac=(NormStep.GT.SThresh)
1683
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
1684
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
1685
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1686
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
1687
        !Print *, 'Grad(1,:):'
1688
        !Print *, Grad(1,:)
1689
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac
1690

    
1691
        GradTmp=GradTmp*FactStep
1692

    
1693
        ! we take care that frozen atoms do not move.
1694
        IF (NFroz.NE.0) THEN
1695
           SELECT CASE (COORD)
1696
           CASE ('ZMAT','MIXED')
1697
              GradTmp(1:IntFroz)=0.d0
1698
           CASE ('CART','HYBRID')
1699
              DO I=1,Nat
1700
                 IF (FrozAtoms(I)) THEN
1701
                    Iat=Order(I)
1702
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1703
                 END IF
1704
              END DO
1705
           CASE('BAKER')
1706
              GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..)
1707
           CASE DEFAULT
1708
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP'
1709
              STOP
1710
           END SELECT
1711
        END IF ! matches IF (NFroz.NE.0) THEN
1712

    
1713
        IF (debug) THEN
1714
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1715
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1716
        END IF
1717

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

    
1730
        SELECT CASE (COORD)
1731
        CASE ('ZMAT','MIXED','BAKER')
1732
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1733
        CASE ('HYBRID','CART')
1734
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1735
        CASE DEFAULT
1736
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1737
           STOP
1738
        END SELECT
1739

    
1740
        IF (debug) THEN
1741
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1742
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1743
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1744
           ELSE
1745
              DO Iat=1,Nat
1746

    
1747
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1748
                 ! It was with indzmat. IndZmat cannot appear here...
1749
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1750
                      XyzGeomF(IGeom,1:3,Iat)
1751
              END DO
1752
           END IF
1753
        END IF ! matches IF (debug) THEN
1754

    
1755
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1756
     ELSE ! Optreac
1757
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1758
     END IF ! matches IF (OptReac) THEN
1759

    
1760
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1761
        IGeom=NGeomF
1762
        if (debug) WRITE(*,*) '******************'
1763
        if (debug) WRITE(*,*) 'Optimizing product'
1764
        if (debug) WRITE(*,*) '******************'
1765
        SELECT CASE (COORD)
1766
        CASE ('ZMAT','MIXED','BAKER')
1767
           Print *, 'Optimizing product'
1768
           SELECT CASE (Step_method)
1769
           CASE ('RFO')
1770
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1771
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1772
              Print *, GradTmp
1773
           CASE ('GDIIS')
1774
              ! GradTmp becomes "step" below and has INTENT(OUT):
1775
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1776
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1777
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1778
              Print *, GradTmp
1779
           CASE ('GEDIIS')
1780
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1781
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1782
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1783
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1784
              Print *, GradTmp
1785
           CASE DEFAULT
1786
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1787
              STOP
1788
           END SELECT
1789
        CASE ('HYBRID','CART')
1790
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1791
        CASE DEFAULT
1792
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1793
           STOP
1794
        END SELECT
1795

    
1796
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1797
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1798
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1799
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1800
        OptProd=.NOT.(NormStep.LE.SThresh)
1801
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1802
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1803
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1804
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1805

    
1806
        GradTmp=GradTmp*FactStep
1807

    
1808
        ! we take care that frozen atoms do not move
1809
        IF (NFroz.NE.0) THEN
1810
           SELECT CASE (COORD)
1811
           CASE ('ZMAT','MIXED','BAKER')
1812
              GradTmp(1:IntFroz)=0.d0
1813
           CASE ('CART','HYBRID')
1814
              DO I=1,Nat
1815
                 IF (FrozAtoms(I)) THEN
1816
                    Iat=Order(I)
1817
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1818
                 END IF
1819
              END DO
1820
           CASE DEFAULT
1821
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1822
              STOP
1823
           END SELECT
1824
        END IF ! matches IF (NFroz.NE.0) THEN
1825

    
1826
        IF (debug) THEN
1827
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
1828
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1829
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1830
           ELSE
1831
              DO Iat=1,Nat
1832
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1833
                      XyzGeomF(IGeom,1:3,Iat)
1834
              END DO
1835
           END IF
1836
        END IF
1837

    
1838
        SELECT CASE (COORD)
1839
        CASE ('ZMAT','MIXED','BAKER')
1840
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1841
        CASE ('HYBRID','CART')
1842
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1843
        CASE DEFAULT
1844
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1845
           STOP
1846
        END SELECT
1847

    
1848
        IF (debug) THEN
1849
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1850
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1851
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1852
           ELSE
1853
              DO Iat=1,Nat
1854
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1855
                      XyzGeomF(IGeom,1:3,Iat)
1856
              END DO
1857
           END IF
1858
        END IF
1859
     ELSE ! Optprod
1860
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1861
     END IF !matches IF (OptProd) THEN 
1862

    
1863
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1864

    
1865
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1866
     !
1867
     !  Optimizing other geometries
1868
     !
1869
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1870

    
1871

    
1872

    
1873
     DO IGeom=2,NGeomF-1 ! matches at L1556
1874
        if (debug) WRITE(*,*) '****************************'
1875
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1876
        if (debug) WRITE(*,*) '****************************'
1877

    
1878
        SELECT CASE (COORD)
1879
        CASE ('ZMAT','BAKER','MIXED')
1880
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
1881
           ! PFL 6 Apr 2008 ->
1882
           ! Special case, if FTan=0. we do not consider Tangent at all.
1883
           ! To DO: add the full treatment of FTan
1884
           if (FTan==0.) GradTmp2=ZeroVec
1885
           ! <- PFL 6 Apr 2008
1886
           if (debug) THEN
1887
              Print *, 'L1500, IntTangent(',IGeom,',:)='
1888
              Print *, IntTangent(IGeom,:)
1889
           END IF
1890
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
1891
           SELECT CASE (Step_method)
1892
           CASE ('RFO')
1893
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1894
           CASE ('GDIIS')
1895
              !The following has no effect at IOpt==1
1896
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
1897
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&        
1898
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1899
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
1900
              Print *, GradTmp
1901
           CASE ('GEDIIS')
1902
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1903
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
1904
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1905
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
1906
              Print *, GradTmp
1907
           CASE DEFAULT
1908
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
1909
              STOP
1910
           END SELECT
1911
        CASE ('HYBRID','CART')
1912
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
1913
           ! so we change it
1914
           DO I=1,Nat
1915
              DO J=1,3
1916
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
1917
              END DO
1918
           END DO
1919
           !           GradTmp2=XyzTangent(IGeom,:)
1920
           ! PFL 6 Apr 2008 ->
1921
           ! Special case, if FTan=0. we do not consider Tangent at all.
1922
           ! To DO: add the full treatment of FTan
1923
           if (FTan==0.) GradTmp2=ZeroVec
1924
           ! <- PFL 6 Apr 2008
1925

    
1926
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1927
        CASE DEFAULT
1928
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
1929
           STOP
1930
        END SELECT
1931

    
1932
        ! we take care that frozen atoms do not move
1933
        IF (NFroz.NE.0) THEN
1934
           SELECT CASE (COORD)
1935
           CASE ('ZMAT','MIXED','BAKER')
1936
              IF (debug) THEN 
1937
                 WRITE(*,*) 'Step computed. Coord=',Coord
1938
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
1939
              END IF
1940
              GradTmp(1:IntFroz)=0.d0
1941
           CASE ('CART','HYBRID')
1942
              if (debug) THEN
1943
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
1944
                 DO I=1,Nat
1945
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
1946
                 END DO
1947
              END IF
1948
              DO I=1,Nat
1949
                 IF (FrozAtoms(I)) THEN
1950
                    if (debug) THEN
1951
                       write(*,*) 'Step Computed. Coord=',Coord
1952
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
1953
                    END IF
1954
                    Iat=Order(I)
1955
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1956
                 END IF
1957
              END DO
1958

    
1959
                 if (debug) THEN
1960
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
1961
                    DO I=1,Nat
1962
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
1963
                    END DO
1964
                 END IF
1965

    
1966
           CASE DEFAULT
1967
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
1968
              STOP
1969
           END SELECT
1970
        END IF !matches IF (NFroz.NE.0) THEN
1971

    
1972
        IF (debug) THEN
1973
           SELECT CASE (COORD)
1974
           CASE ('ZMAT','MIXED','BAKER')
1975
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
1976
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1977
           CASE ('CART','HYBRID')
1978
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
1979
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1980
           CASE DEFAULT
1981
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
1982
              STOP
1983
           END SELECT
1984
        END IF ! matches if (debug) THEN
1985

    
1986
        ! we check that the maximum displacement is not greater than Smax
1987
        If (debug) WRITE(*,*) "Fini before test:",Fini
1988
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1989
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
1990
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1991
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
1992

    
1993
        GradTmp=GraDTmp*FactStep
1994

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

    
1998
        IF (debug) THEN
1999
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
2000
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2001
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
2002
           ELSE
2003
              DO Iat=1,Nat
2004
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
2005
                      XyzGeomF(IGeom,1:3,Iat)
2006
              END DO
2007
           END IF
2008
        END IF ! matches if (debug) THEN
2009

    
2010
        SELECT CASE (COORD)
2011
        CASE ('ZMAT')
2012
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2013
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2014
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2015
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2016
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2017
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
2018
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2019
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
2020
           Idx=4
2021
           DO Iat=4,Nat
2022
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2023
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2024
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2025
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2026
              Idx=Idx+3
2027
           END DO
2028

    
2029
        CASE ('BAKER')
2030
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2031
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2032
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2033
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
2034

    
2035
        CASE ('MIXED')
2036
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2037
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2038
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2039
           DO Iat=1,NCart
2040
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
2041
           END DO
2042
           Idx=3*NCart+1
2043
           IF (NCart.EQ.1) THEN
2044
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2045
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
2046
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2047
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
2048

    
2049
              Idx=7
2050
           END IF
2051
           IF (NCart.EQ.2) THEN
2052
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2053
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
2054
              Idx=9
2055
           END IF
2056

    
2057
           
2058
           DO Iat=max(NCart+1,4),Nat
2059
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2060
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2061
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2062
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2063
              Idx=Idx+3
2064
           END DO
2065
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2066

    
2067
        CASE ('HYBRID','CART')
2068
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
2069
        CASE DEFAULT
2070
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
2071
           STOP
2072
        END SELECT
2073

    
2074
        IF (debug) THEN
2075
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
2076
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2077
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
2078
           ELSE
2079
              DO Iat=1,Nat
2080
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
2081
                      XyzGeomF(IGeom,1:3,Iat)
2082
              END DO
2083
           END IF
2084
        END IF ! matches if (debug) THEN
2085

    
2086
        ! We project out the tangential part of the gradient to know if we are done
2087
        GradTmp=Grad(IGeom,:)
2088
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
2089
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
2090
        SELECT CASE (COORD)
2091
        CASE ('ZMAT','MIXED','BAKER')
2092
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2093
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
2094
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2095
           Ca=S/(sqrt(Norm)*NormGrad)
2096
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2097
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2098
        CASE ('CART','HYBRID')
2099
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2100
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
2101
           Ca=S/(sqrt(Norm)*NormGrad)
2102
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2103
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2104
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2105
        CASE DEFAULT
2106
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
2107
           STOP
2108
        END SELECT
2109
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
2110
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
2111
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
2112
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
2113

    
2114
     END DO ! matches DO IGeom=2,NGeomF-1
2115
     ! we save the old gradients
2116
     GradOld=Grad
2117
     EPathOld=Energies
2118

    
2119

    
2120

    
2121
     ! Shall we re-parameterize the path ?
2122
     ! PFL 2007/Apr/10 ->
2123
     ! We call PathCreate in all cases, it will take care of the 
2124
     ! reparameterization, as well as calculating the tangents.
2125
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
2126
     XyzGeomI=XyzGeomF
2127
     IntCoordI=IntCoordF
2128

    
2129
     Call PathCreate(IOpt)
2130
     !     END IF
2131
     ! <- PFL 2007/Apr/10
2132

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

    
2136
        ! We have the new path, we calculate its energy and gradient
2137

    
2138
        Call EgradPath(IOpt)
2139
        !IF(IOPT .EQ. 10) Then
2140
        !         Print *, 'Stopping at my checkpoint.'
2141
        !           STOP !This is my temporary checkpoint.
2142
        !ENDIF
2143

    
2144
        ! PFL 31 Mar 2008 ->
2145
        ! Taken from v3.99 but we added a flag...
2146
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
2147
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2148
        ! else it is decreased by 0.8
2149

    
2150
        if (DynMaxStep) THEN
2151
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
2152
           if (OptReac) THEN
2153
              If (Energies(1)<EPathOld(1)) Then
2154
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
2155
              ELSE
2156
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
2157
              END IF
2158
           END IF
2159

    
2160
           if (OptProd) THEN
2161
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
2162
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
2163
              ELSE
2164
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
2165
              END IF
2166
           END IF
2167

    
2168
           DO IGeom=2,NGeomF-1
2169
              If (Energies(IGeom)<EPathOld(IGeom)) Then
2170
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
2171
              ELSE
2172
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
2173
              END IF
2174
           END DO
2175
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2176
        END IF ! test on DynMaxStep
2177
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2178
        ! <- PFL 31 MAr 2008
2179
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
2180
        !  It should get updated for other cases too.
2181

    
2182
        DO IGeom=1,NGeomF
2183
           SELECT CASE (COORD)
2184
           CASE('ZMAT')
2185
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2186
              GeomTmp=IntCoordF(IGeom,:)
2187
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2188
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2189
                   OrderInv(indzmat(2,2)),Geomtmp(1)
2190
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2191
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
2192
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
2193
              Idx=4
2194
              DO Iat=4,Nat
2195
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2196
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
2197
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
2198
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
2199
                 Idx=Idx+3
2200
              END DO
2201
           CASE('BAKER')
2202
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2203
              !GeomTmp=IntCoordF(IGeom,:)
2204
           CASE('CART','HYBRID','MIXED')
2205
              WRITE(IOOUT,'(1X,I5)') Nat
2206
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2207
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2208
              DO I=1,Nat
2209
                 Iat=I
2210
                 If (renum) Iat=OrderInv(I)
2211
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
2212
              END DO
2213
           CASE DEFAULT
2214
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
2215
              STOP
2216
           END SELECT
2217
        END DO ! matches DO IGeom=1,NGeomF
2218

    
2219
        Call Write_path(IOpt)
2220
! Shall we analyze the geometries ?
2221
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
2222

    
2223
        ! We update the Hessian matrices
2224
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
2225
        ! First using the displacement from the old path
2226
        IGeom0=2
2227
        IGeomF=NGeomF-1
2228
        IF (OptReac) IGeom0=1
2229
        IF (OptProd) IGeomF=NGeomF
2230

    
2231
        IF (FCalcHess) THEN
2232
           DO IGeom=IGeom0,IGeomF
2233
              SELECT CASE (COORD)
2234
              CASE ('ZMAT','MIXED','BAKER')
2235
                 GeomTmp2=IntCoordF(IGeom,:)
2236
              CASE ('CART','HYBRID')
2237
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2238
              CASE DEFAULT
2239
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2240
                 STOP
2241
              END SELECT
2242
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
2243
           END DO
2244
        ELSE
2245
           DO IGeom=IGeom0,IGeomF
2246
              GeomTmp=GeomOld_all(IGeom,:)
2247
              SELECT CASE (COORD)
2248
              CASE ('ZMAT','MIXED','BAKER')
2249
                 GeomTmp2=IntCoordF(IGeom,:)
2250
              CASE ('CART','HYBRID')
2251
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2252
              CASE DEFAULT
2253
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2254
                 STOP
2255
              END SELECT
2256
              GeomTmp=GeomTmp2-GeomTmp
2257
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
2258
              HessTmp=Hess(IGeom,:,:)
2259
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2260
              Hess(IGeom,:,:)=HessTmp
2261
           END DO ! matches DO IGeom=IGeom0,IGeomF
2262

    
2263
           ! Second using the neighbour points
2264
           IF (HupNeighbour) THEN
2265
              ! Only one point for end points.
2266
              IF (OptReac) THEN
2267
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2268
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
2269
                 ELSE
2270
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
2271
                 END IF
2272
                 GradTmp=Grad(1,:)-Grad(2,:)
2273
                 HessTmp=Hess(1,:,:)
2274
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2275
                 Hess(1,:,:)=HessTmp
2276
              END IF
2277

    
2278
              IF (OptProd) THEN
2279
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2280
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
2281
                 ELSE
2282
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
2283
                 END IF
2284
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
2285
                 HessTmp=Hess(NGeomF,:,:)
2286
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2287
                 Hess(NGeomF,:,:)=HessTmp
2288
              END IF
2289

    
2290
              ! Two points for the rest of the path
2291
              DO IGeom=2,NGeomF-1
2292
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
2293
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
2294
                 ELSE
2295
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
2296
                 END IF
2297
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
2298
                 HessTmp=Hess(IGeom,:,:)
2299
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2300

    
2301
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2302
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
2303
                 ELSE
2304
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
2305
                 END IF
2306
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
2307

    
2308
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2309
                 Hess(IGeom,:,:)=HessTmp
2310
              END DO ! matches DO IGeom=2,NGeomF-1
2311
           END IF !matches IF (HupNeighbour) THEN
2312
        END IF ! matches IF (FCalcHess)
2313
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
2314
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
2315

    
2316
!   IF (PROG=="GAUSSIAN") THEN
2317
!      DEALLOCATE(Gauss_paste)
2318
!      DEALLOCATE(Gauss_root)
2319
!      DEALLOCATE(Gauss_comment)
2320
!      DEALLOCATE(Gauss_end)
2321
!   END IF
2322

    
2323
  DEALLOCATE(XyzGeomF, IntCoordF)
2324
  DEALLOCATE(XyzGeomI, IntCoordI)
2325
  DEALLOCATE(XyzTangent,IntTangent)
2326
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2327
  DEALLOCATE(GeomOld_all)
2328

    
2329
  if (PROG=="GAUSSIAN") THEN
2330
     ! We de-allocate the Gauss_XX lists
2331
     ! transverse the list and deallocate each node
2332
     previous => Gauss_root ! make current point to head of list
2333
     DO
2334
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2335
        current => previous%next ! make list point to next node of head
2336
        DEALLOCATE(previous) ! deallocate current head node
2337
        previous => current  ! make current point to new head
2338
     END DO
2339

    
2340
     previous => Gauss_comment ! make current point to head of list
2341
     DO
2342
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2343
        current => previous%next ! make list point to next node of head
2344
        DEALLOCATE(previous) ! deallocate current head node
2345
        previous => current  ! make current point to new head
2346
     END DO
2347

    
2348
     previous => Gauss_end ! make current point to head of list
2349
     DO
2350
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2351
        current => previous%next ! make list point to next node of head
2352
        DEALLOCATE(previous) ! deallocate current head node
2353
        previous => current  ! make current point to new head
2354
     END DO
2355

    
2356

    
2357
  END IF
2358

    
2359
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2360
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2361
  IF (COORD.EQ.'BAKER') THEN
2362
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2363
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2364
  END IF
2365

    
2366
  If (FPBC)   DeAllocate(XGeomRefPBC,YGeomRefPBC,ZGeomRefPBC)
2367
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2368

    
2369
  Close(IOIN)
2370
  Close(IOOUT)
2371
  IF (AnaGeom) Close(IODAT)
2372
  IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot)
2373

    
2374
END PROGRAM PathOpt