Statistiques
| Révision :

root / src / Path.f90 @ 9

Historique | Voir | Annoter | Télécharger (89,46 ko)

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

    
232
PROGRAM PathOpt
233

    
234
  use VarTypes
235
  use Path_module
236
  use Io_module
237

    
238
  IMPLICIT NONE
239

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

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

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

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

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

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

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

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

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

    
283
  LOGICAL :: Debug, Fini,Flag_Opt_Geom
284

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

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

    
296

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

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

    
307

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

    
310
       IMPLICIT NONE
311

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

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

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

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

    
333

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

    
336
      Use VarTypes
337

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

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

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

    
353
      Use VarTypes
354
      Use io_module
355

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

    
361
    END SUBROUTINE die
362

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

    
365
      Use VarTypes
366
      Use io_module
367

    
368
      implicit none
369

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

    
374
    END SUBROUTINE Warning
375

    
376

    
377
  END INTERFACE
378

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

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

    
392

    
393
  Flag_Opt_Geom = .FALSE. ! Initialized.
394

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

    
401
  ! We read the files  names
402
!  SELECT CASE (iargc())
403
  SELECT CASE (command_argument_count())
404
  CASE (2)
405
     call getarg(1,FileIn)
406
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
407
     call getarg(2,FileOut)
408
     OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN')
409
  CASE (1)
410
     call getarg(1,FileIn)
411
     IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
412
        WRITE(*,*) "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
    CASE ("TM","TURBOMOLE")
835
     Prog="TURBOMOLE"
836
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
837
     UnitProg="au"
838
     ConvE=au2kcal
839
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
840
     ProgExe=""
841
     UnitProg="au"
842
     ConvE=au2kcal
843
    CASE DEFAULT
844
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
845
     STOP
846
  END SELECT 
847

    
848
  if (Input.EQ.'EMPTY') THEN
849
     Select Case (Prog)
850
       CASE ("VASP")
851
          Input=Prog
852
       CASE ("CHAMFRE")
853
          Input="VASP"
854
       CASE DEFAULT
855
          Input='XYZ'
856
     END SELECT
857
    WRITE(*,*) "Input has been set to the default: ",INPUT
858
 END IF
859

    
860
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
861

    
862
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
863
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
864
     STOP
865
  END IF
866

    
867
  if (OptGeom.GE.1) THEN
868
     FPrintGeom=.FALSE.
869
     If (GeomFile/='EMPTY') THEN
870
        FPrintGeom=.TRUE.
871
     END IF
872
     IF (IGeomRef.LE.0) THEN
873
        IGeomRef=OptGeom
874
     ELSE
875
        IF (IGeomRef.NE.OptGeom) THEN
876
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
877
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
878

    
879
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
880
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
881
           STOP
882
        END IF
883
     END IF
884
  END IF
885

    
886
  ALLOCATE(Cart(Nat))
887
  Cart=0
888

    
889
  ALLOCATE(Frozen(Nat))
890
  Frozen=0
891

    
892
  IF (FCart) THEN
893
! We rewind the file to be sure to browse all of it to read the Cart namelist
894
     REWIND(IOIN)
895
     List=0
896
     READ(IOIN,cartlist)
897
     IF (COORD.NE.'CART') THEN
898
        Cart=List(1:Nat)
899
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
900
     ELSE
901
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
902
        WRITE(*,*) "I will discard the cartlist"
903
        Cart=0
904
     END IF
905
  END IF
906

    
907
  if (FFrozen) THEN
908

    
909
     REWIND(IOIN)
910
     List=0
911
     READ(IOIN,Frozenlist)
912
     Frozen=List(1:Nat)
913
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
914
  END IF
915

    
916
  If (FCart) Call Expand(Cart,NCart,Nat)
917
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
918

    
919
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
920
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
921

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

    
949
        
950

    
951
  if (NMaxPtPath.EQ.-1) then
952
     NMaxPtPath=min(50*NGeomF,2000)
953
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
954
  end if
955

    
956
  !  NMaxPtPath=10000
957

    
958
  ! We ensure that FTan is in [0.:1.]
959
  FTan=min(abs(FTan),1.d0)
960

    
961
! PFL 2011 Mar 14 ->
962
! Added some printing for analyses with Anapath
963
  if (debug) THEN
964
     WRITE(IOOUT,path)
965
  ELSE
966
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC
967
! and PATHNAME 
968
     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat
969
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI
970
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF
971
     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC
972
     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)
973
  END IF
974

    
975
  FTan=1.-FTan
976

    
977
  !Thresholds...
978
  if (SThresh.LE.0) SThresh=0.01d0
979
  If (GThresh.LE.0) GThresh=SThresh/4.
980
  SThreshRms=Sthresh*2./3.
981
  GThreshRms=Gthresh*2./3.
982

    
983
  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
984
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
985
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
986
  ALLOCATE(MassAt(NAt),OrderInv(Nat))
987

    
988
  AtName=""
989
  MassAt=1.
990

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

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

    
1037
  ! We read the initial geometries
1038
  Call Read_Geom(input)
1039

    
1040
  ! We convert atome names into atom mass number
1041
  DO I=1,NAt
1042
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
1043
     Atome(I)=ConvertNumAt(AtName(I))
1044
  END DO
1045

    
1046
 If (AnaGeom) THEN
1047
    If (NbVar>0) THEN
1048
! We print the description of the varialbes in AnaFile
1049
       Call PrintAnaList(IoDat)
1050
       If (OptGeom<=0) THEN
1051
          WRITE(IoDat,'(A)') "# s Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
1052
       ELSE
1053
          WRITE(IoDat,'(A)') "# Iter Var1 ... VarN Erel(kcal/mol) E (" // TRIM(UnitProg) //")"
1054
       END IF
1055
    ELSE
1056
       If (OptGeom<=0) THEN
1057
          WRITE(IoDat,'(A)') "#  s    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"
1058
       ELSE
1059
          WRITE(IoDat,'(A)') "#  Iter    Erel(kcal/mol)   E (" // TRIM(UnitProg) //")"       
1060
       END IF
1061
    END IF
1062
 END IF
1063

    
1064
  ! PFL 23 Sep 2008 ->
1065
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
1066
  ! between the first geometry and the others.
1067
  IF (NFroz.GT.0) THEN
1068
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
1069
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
1070
     ListAtFroz=.FALSE.
1071

    
1072
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
1073
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
1074
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
1075

    
1076
     SELECT CASE (NFroz)
1077
        ! We have Frozen atoms
1078
        ! PFL 28 Oct 2008 ->
1079
        ! It might happen that the geometries are translated/rotated.
1080
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
1081
        ! three frozen atoms.
1082

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

    
1160
     DO I=2,NGeomI
1161
        ! First, we align the Ith geometry on the first one, using only
1162
        ! the positions of the "frozen" atoms. (see above: it might be that
1163
        ! ListAtFroz contains non frozen atoms used to align the geometries
1164
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
1165
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
1166
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
1167
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
1168

    
1169
        FChkFrozen=.FALSE.
1170
        DO J=1,NFroz
1171
           IAt=Frozen(J)
1172
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
1173
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
1174
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
1175
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
1176
              FChkFrozen=.TRUE.
1177
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
1178
              WRITE(*,*) "from geometry 1 to geometry ",I
1179
           END IF
1180
        END DO
1181
     END DO
1182
     IF (FChkFrozen) THEN
1183
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
1184
        STOP
1185
     END IF
1186
  END IF
1187

    
1188

    
1189
  N3at=3*Nat
1190
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
1191

    
1192

    
1193
  Call ReadInput
1194

    
1195
  IF (COORD=="MIXED") Call TestCart(AutoCart)
1196

    
1197
  SELECT CASE(NFroz)
1198
  CASE (2)
1199
     IntFroz=1
1200
  CASE (3)
1201
     IntFroz=3
1202
  CASE (4:)
1203
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
1204
  END SELECT
1205

    
1206
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
1207
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
1208
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
1209

    
1210
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
1211
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
1212
     COORD="ZMAT"
1213
  END IF
1214

    
1215
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
1216

    
1217
  SELECT CASE(COORD)
1218
  CASE ('ZMAT')
1219
     NCoord=Nfree
1220
     ALLOCATE(dzdc(3,nat,3,nat))
1221
  CASE ('MIXED')
1222
     SELECT CASE (NCart+NFroz)
1223
     CASE (1) 
1224
        NCoord=N3At-3
1225
     CASE (2)
1226
        NCoord=N3At-1
1227
     CASE (3:)
1228
        NCoord=N3At
1229
     END SELECT
1230
     ALLOCATE(dzdc(3,nat,3,nat))
1231
  CASE ('BAKER')
1232
     Symmetry_elimination=0
1233
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
1234
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
1235
     ALLOCATE(BTransInv(NCoord,3*Nat))
1236
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
1237
  CASE ('HYBRID')
1238
     NCoord=N3at
1239
  CASE ('CART')
1240
     NCoord=N3at
1241
  END SELECT
1242

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

    
1245
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
1246
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
1247
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
1248
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
1249
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
1250
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
1251
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
1252
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
1253

    
1254
  ZeroVec=0.d0
1255
  DO I =1, NGeomF
1256
     IntTangent(I,:)=0.d0
1257
  END DO
1258

    
1259
  if (debug) THEN
1260
     Print *, 'L886, IntTangent(1,:)='
1261
     Print *, IntTangent(1,:)
1262
  END IF
1263

    
1264
  if (.NOT.OptReac) Energies(1)=Ereac
1265
  if (.NOT.OptProd) Energies(NGeomF)=EProd
1266
  MaxStep=SMax
1267

    
1268
  ! We initialize the mass array
1269
  if (MW) THEN
1270
     WRITE(*,*) 'Working in MW coordinates'
1271
     DO I=1,Nat
1272
        massAt(I)=Mass(Atome(I))
1273
     END DO
1274
  ELSE
1275
     DO I=1,Nat
1276
        MassAt(I)=1.d0
1277
     END DO
1278
  END IF
1279

    
1280
  WRITE(*,*) "Prog=",TRIM(PROG)
1281

    
1282
  ALLOCATE(FrozAtoms(Nat))
1283

    
1284
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
1285
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
1286
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
1287
  IF (NFroz.EQ.0) THEN
1288
     FrozAtoms=.TRUE.
1289
  ELSE
1290
     I=1
1291
     NFr=0
1292
     FrozAtoms=.False.
1293
     DO WHILE (Frozen(I).NE.0)
1294
        IF (Frozen(I).LT.0) THEN
1295
           DO J=Frozen(I-1),abs(Frozen(I))
1296
              IF (.NOT.FrozAtoms(J)) THEN
1297
                 NFr=NFr+1
1298
                 FrozAtoms(J)=.TRUE.
1299
              END IF
1300
           END DO
1301
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1302
           FrozAtoms(Frozen(I))=.TRUE.
1303
           NFr=NFr+1
1304
        END IF
1305
        I=I+1
1306
     END DO
1307
     IF (NFr.NE.NFroz) THEN
1308
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1309
        STOP
1310
     END IF
1311
  END IF
1312

    
1313
!  if (debug) THEN
1314
     !Print *, 'L968, IntTangent(1,:)='
1315
     !Print *, IntTangent(1,:)
1316
!  END IF
1317

    
1318
  ! We have read everything we needed... time to work :-)
1319
  IOpt=0
1320
  FirstTimePathCreate = .TRUE. ! Initialized.
1321
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1322
  FirstTimePathCreate = .FALSE.
1323

    
1324
  if (debug) THEN
1325
     Print *, 'L980, IntTangent(1,:)='
1326
     Print *, IntTangent(1,:)
1327
  END IF
1328

    
1329
  ! PFL 30 Mar 2008 ->
1330
  ! The following is now allocated in Calc_Baker. This is more logical to me
1331
  !   IF (COORD .EQ. "BAKER") Then
1332
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1333
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1334
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1335
  !   END IF
1336
  ! <- PFL 30 mar 2008
1337

    
1338
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1339
  !
1340
  !
1341
  ! OptGeom is the index of the geometry which is to be optimized.
1342
  IF (Optgeom.GT.0) THEN
1343
     Flag_Opt_Geom=.TRUE.
1344
     SELECT CASE(Coord)
1345
     CASE ('CART','HYBRID')
1346
!!! CAUTION : PFL 29.AUG.2008 ->
1347
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1348
        ! This should be made  consistent !!!!!!!!!!!!!!!
1349
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1350
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1351
!!! <- CAUTION : PFL 29.AUG.2008
1352
     CASE ('ZMAT','MIXED')
1353
        !GeomTmp=IntCoordF(OptGeom,:)
1354
        GeomTmp=IntCoordI(OptGeom,:)
1355
     CASE ('BAKER')
1356
        !GeomTmp=IntCoordF(OptGeom,:)
1357
        GeomTmp=IntCoordI(OptGeom,:)
1358
     CASE DEFAULT
1359
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1360
        STOP
1361
     END SELECT
1362

    
1363
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1364
     Flag_Opt_Geom = .TRUE.
1365
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1366

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

    
1373
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1374
  ! End of GeomOpt
1375
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1376

    
1377
  IF (PathOnly) THEN
1378
     Call Header("PathOnly=.T. , Stop here")
1379
     Call Write_path(-1)
1380
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1381
     STOP
1382
  END IF
1383

    
1384
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1385

    
1386
  DEALLOCATE(XyzGeomI,IntCoordI)
1387
  NGeomI=NGeomF
1388
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1389

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

    
1395
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1396
     DO IGeom=1,NGeomF
1397
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1398
        DO I=1,3*Nat-6
1399
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1400
        END DO
1401
     END DO
1402
  END IF
1403

    
1404
  ! To calculate B^T^-1 for all extrapolated final geometries:
1405
  ! PFL 18 Jan 2008 ->
1406
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1407
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1408
  ! <-- PFL 18 Jan 2008
1409

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

    
1419
  IOpt=0
1420
  Call EgradPath(IOpt)
1421

    
1422
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1423
  Call Write_path(IOpt)
1424
! Shall we analyze the geometries ?
1425
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1426

    
1427
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1428

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

    
1433
  ! By default, Hess=Id
1434
  Hess=0.
1435
  IF(Coord .EQ. "ZMAT") Then
1436
     ! We use the fact that we know that approximate good hessian values
1437
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1438
     DO IGeom=1,NGeomF
1439
        Hess(IGeom,1,1)=0.5d0
1440
        Hess(IGeom,2,2)=0.5d0
1441
        Hess(IGeom,3,3)=0.1d0
1442
        DO J=1,Nat-3
1443
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1444
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1445
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1446
        END DO
1447
        IF (HInv) THEN
1448
           DO I=1,NCoord
1449
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1450
           END DO
1451
        END IF
1452
     END DO
1453
  ELSE
1454
     DO I=1,NCoord
1455
        DO J=1,NGeomF
1456
           Hess(J,I,I)=1.d0
1457
        END DO
1458
     END DO
1459
  END IF
1460

    
1461
  IF (COORD .EQ. "BAKER") THEN
1462
     ! UMat(NPrim,NCoord)
1463
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1464
     Hprim=0.d0
1465
     ScanCoord=>Coordinate
1466
     I=0
1467
     DO WHILE (Associated(ScanCoord%next))
1468
        I=I+1
1469
        SELECT CASE (ScanCoord%Type)
1470
        CASE ('BOND')
1471
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1472
           Hprim(I,I) = 0.5d0
1473
           !END DO
1474
        CASE ('ANGLE')
1475
           Hprim(I,I) = 0.2d0
1476
        CASE ('DIHEDRAL')
1477
           Hprim(I,I) = 0.1d0
1478
        END SELECT
1479
        ScanCoord => ScanCoord%next
1480
     END DO
1481

    
1482
     ! HprimU = H^prim U:
1483
     HprimU=0.d0
1484
     DO I=1, 3*Nat-6
1485
        DO J=1,NPrim
1486
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1487
        END DO
1488
     END DO
1489

    
1490
     ! Hess = U^T HprimU = U^T H^prim U:
1491
     Hess=0.d0
1492
     DO I=1, 3*Nat-6
1493
        DO J=1,NPrim
1494
           ! UMat^T is needed below.
1495
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1496
        END DO
1497
     END DO
1498
     DO K=2,NGeomF
1499
        Hess(K,:,:)=Hess(1,:,:)
1500
     END DO
1501
     DEALLOCATE(Hprim,HprimU)
1502
  end if !  ends IF COORD=BAKER
1503
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1504
  ALLOCATE(HessTmp(NCoord,NCoord))
1505

    
1506
  IF (IniHUp) THEN
1507
     IF (FCalcHess) THEN
1508
        ! The numerical calculation of the Hessian has been switched on
1509
        DO IGeom=1,NGeomF
1510
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1511
              GeomTmp=IntCoordF(IGeom,:)
1512
           ELSE
1513
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1514
           END IF
1515
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1516
        END DO
1517
     ELSE
1518
        ! We generate a forward hessian and a backward one and we mix them.
1519
        ! First, the forward way:
1520
        DO IGeom=2,NGeomF-1
1521
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1522
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1523
           ELSE
1524
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1525
           END IF
1526
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1527
           HessTmp=Hess(IGeom-1,:,:)
1528
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1529
           Hess(IGeom,:,:)=HessTmp
1530
        END DO
1531

    
1532
        ! Now, the backward way:
1533
        HessTmp=Hess(NGeomF,:,:)
1534
        DO IGeom=NGeomF-1,2,-1
1535
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1536
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1537
           ELSE
1538
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1539
           END IF
1540
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1541
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1542
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1543
           ca=cos(alpha)
1544
           sa=sin(alpha)
1545
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1546
        END DO
1547
     END IF ! matches IF (FCalcHess)
1548
  END IF ! matches IF (IniHUp) THEN
1549

    
1550
  ! For debug purposes, we diagonalize the Hessian matrices
1551
  IF (debug) THEN
1552
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1553
     DO I=1,NGeomF
1554
        WRITE(*,*) "DBG Main - Point ",I
1555
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1556
        DO J=1,NCoord
1557
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1558
        END DO
1559
     END DO
1560
  END IF ! matches IF (debug) THEN
1561

    
1562
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1563
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1564
  !
1565
  ! Main LOOP : Optimization of the Path
1566
  !
1567
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1568
  if (debug)   Print *, 'NGeomF=', NGeomF
1569
  allocation_flag = .TRUE.
1570

    
1571
  Fini=.FALSE.
1572

    
1573
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1574
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1575
     IOpt=IOpt+1
1576

    
1577
     SELECT CASE (COORD)
1578
     CASE ('ZMAT','MIXED','BAKER')
1579
        GeomOld_all=IntCoordF
1580
     CASE ('CART','HYBRID')
1581
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1582
     CASE DEFAULT
1583
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1584
        STOP
1585
     END SELECT
1586

    
1587
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1588
          valid("printtangent")) THEN
1589
        WRITE(*,*) "Visualization of tangents"
1590
        Idx=min(12,NCoord)
1591
        DO I=1,NGeomF
1592
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1593
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1594
           WRITE(*,*) 
1595
        END DO
1596
        WRITE(*,*) "END of tangents"
1597
     END IF
1598

    
1599
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1600
          valid("printgrad")) THEN
1601
        WRITE(*,*) "Visualization of gradients"
1602
        Idx=min(12,NCoord)
1603
        DO I=1,NGeomF
1604
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1605
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1606
           WRITe(*,*) 
1607
        END DO
1608
        WRITE(*,*) "END of gradients"
1609
     END IF
1610

    
1611
     Fini=.TRUE.
1612
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1613
        IGeom=1
1614
        if (debug) WRITE(*,*) '**************************************'
1615
        if (debug) WRITE(*,*) 'Optimizing reactant'
1616
        if (debug) WRITE(*,*) '**************************************'
1617
        SELECT CASE (COORD)
1618
        CASE ('ZMAT','MIXED','BAKER')
1619
           SELECT CASE (Step_method)
1620
           CASE ('RFO')
1621
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
1622
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1623
              Print *, GradTmp
1624
           CASE ('GDIIS')
1625
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1626
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1627
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1628
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1629
              Print *, GradTmp
1630
           CASE ('GEDIIS')
1631
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1632
              ! Energies are updated in EgradPath.f90
1633
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1634
                   NCoord,allocation_flag,ZeroVec)
1635
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1636
              !GradTmp = GeomTmp - IntCoordF(1,:)
1637
              !Print *, 'Old Geometry:'
1638
              !Print *, IntCoordF(1,:)
1639
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1640
              Print *, GradTmp
1641
              !Print *, 'GeomTmp:'
1642
              !Print *, GeomTmp
1643
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1644
           CASE DEFAULT
1645
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1646
              STOP
1647
           END SELECT
1648
        CASE ('HYBRID','CART')
1649
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1650
        CASE DEFAULT
1651
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1652
           STOP
1653
        END SELECT
1654

    
1655
        IF (debug) THEN
1656
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1657
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1658
        END IF
1659

    
1660
        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
1661
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1662
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
1663
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1664
        OptReac=(NormStep.GT.SThresh)
1665
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
1666
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
1667
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1668
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
1669
        !Print *, 'Grad(1,:):'
1670
        !Print *, Grad(1,:)
1671
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac
1672

    
1673
        GradTmp=GradTmp*FactStep
1674

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

    
1695
        IF (debug) THEN
1696
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1697
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1698
        END IF
1699

    
1700
        IF (debug) THEN
1701
           !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1'
1702
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1703
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1704
           ELSE
1705
              DO Iat=1,Nat
1706
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1707
                      XyzGeomF(IGeom,1:3,Iat)
1708
              END DO
1709
           END IF
1710
        END IF ! matches IF (debug) THEN
1711

    
1712
        SELECT CASE (COORD)
1713
        CASE ('ZMAT','MIXED','BAKER')
1714
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1715
        CASE ('HYBRID','CART')
1716
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1717
        CASE DEFAULT
1718
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1719
           STOP
1720
        END SELECT
1721

    
1722
        IF (debug) THEN
1723
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1724
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1725
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1726
           ELSE
1727
              DO Iat=1,Nat
1728

    
1729
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1730
                 ! It was with indzmat. IndZmat cannot appear here...
1731
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1732
                      XyzGeomF(IGeom,1:3,Iat)
1733
              END DO
1734
           END IF
1735
        END IF ! matches IF (debug) THEN
1736

    
1737
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1738
     ELSE ! Optreac
1739
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1740
     END IF ! matches IF (OptReac) THEN
1741

    
1742
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1743
        IGeom=NGeomF
1744
        if (debug) WRITE(*,*) '******************'
1745
        if (debug) WRITE(*,*) 'Optimizing product'
1746
        if (debug) WRITE(*,*) '******************'
1747
        SELECT CASE (COORD)
1748
        CASE ('ZMAT','MIXED','BAKER')
1749
           Print *, 'Optimizing product'
1750
           SELECT CASE (Step_method)
1751
           CASE ('RFO')
1752
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1753
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1754
              Print *, GradTmp
1755
           CASE ('GDIIS')
1756
              ! GradTmp becomes "step" below and has INTENT(OUT):
1757
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1758
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1759
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1760
              Print *, GradTmp
1761
           CASE ('GEDIIS')
1762
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1763
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1764
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1765
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1766
              Print *, GradTmp
1767
           CASE DEFAULT
1768
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1769
              STOP
1770
           END SELECT
1771
        CASE ('HYBRID','CART')
1772
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1773
        CASE DEFAULT
1774
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1775
           STOP
1776
        END SELECT
1777

    
1778
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1779
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1780
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1781
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1782
        OptProd=.NOT.(NormStep.LE.SThresh)
1783
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1784
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1785
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1786
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1787

    
1788
        GradTmp=GradTmp*FactStep
1789

    
1790
        ! we take care that frozen atoms do not move
1791
        IF (NFroz.NE.0) THEN
1792
           SELECT CASE (COORD)
1793
           CASE ('ZMAT','MIXED','BAKER')
1794
              GradTmp(1:IntFroz)=0.d0
1795
           CASE ('CART','HYBRID')
1796
              DO I=1,Nat
1797
                 IF (FrozAtoms(I)) THEN
1798
                    Iat=Order(I)
1799
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1800
                 END IF
1801
              END DO
1802
           CASE DEFAULT
1803
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1804
              STOP
1805
           END SELECT
1806
        END IF ! matches IF (NFroz.NE.0) THEN
1807

    
1808
        IF (debug) THEN
1809
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
1810
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1811
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1812
           ELSE
1813
              DO Iat=1,Nat
1814
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1815
                      XyzGeomF(IGeom,1:3,Iat)
1816
              END DO
1817
           END IF
1818
        END IF
1819

    
1820
        SELECT CASE (COORD)
1821
        CASE ('ZMAT','MIXED','BAKER')
1822
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1823
        CASE ('HYBRID','CART')
1824
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1825
        CASE DEFAULT
1826
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1827
           STOP
1828
        END SELECT
1829

    
1830
        IF (debug) THEN
1831
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1832
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1833
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1834
           ELSE
1835
              DO Iat=1,Nat
1836
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1837
                      XyzGeomF(IGeom,1:3,Iat)
1838
              END DO
1839
           END IF
1840
        END IF
1841
     ELSE ! Optprod
1842
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1843
     END IF !matches IF (OptProd) THEN 
1844

    
1845
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1846

    
1847
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1848
     !
1849
     !  Optimizing other geometries
1850
     !
1851
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1852

    
1853

    
1854

    
1855
     DO IGeom=2,NGeomF-1 ! matches at L1556
1856
        if (debug) WRITE(*,*) '****************************'
1857
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1858
        if (debug) WRITE(*,*) '****************************'
1859

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

    
1908
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1909
        CASE DEFAULT
1910
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
1911
           STOP
1912
        END SELECT
1913

    
1914
        ! we take care that frozen atoms do not move
1915
        IF (NFroz.NE.0) THEN
1916
           SELECT CASE (COORD)
1917
           CASE ('ZMAT','MIXED','BAKER')
1918
              IF (debug) THEN 
1919
                 WRITE(*,*) 'Step computed. Coord=',Coord
1920
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
1921
              END IF
1922
              GradTmp(1:IntFroz)=0.d0
1923
           CASE ('CART','HYBRID')
1924
              if (debug) THEN
1925
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
1926
                 DO I=1,Nat
1927
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
1928
                 END DO
1929
              END IF
1930
              DO I=1,Nat
1931
                 IF (FrozAtoms(I)) THEN
1932
                    if (debug) THEN
1933
                       write(*,*) 'Step Computed. Coord=',Coord
1934
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
1935
                    END IF
1936
                    Iat=Order(I)
1937
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1938
                 END IF
1939
              END DO
1940

    
1941
                 if (debug) THEN
1942
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
1943
                    DO I=1,Nat
1944
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
1945
                    END DO
1946
                 END IF
1947

    
1948
           CASE DEFAULT
1949
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
1950
              STOP
1951
           END SELECT
1952
        END IF !matches IF (NFroz.NE.0) THEN
1953

    
1954
        IF (debug) THEN
1955
           SELECT CASE (COORD)
1956
           CASE ('ZMAT','MIXED','BAKER')
1957
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
1958
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1959
           CASE ('CART','HYBRID')
1960
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
1961
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1962
           CASE DEFAULT
1963
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
1964
              STOP
1965
           END SELECT
1966
        END IF ! matches if (debug) THEN
1967

    
1968
        ! we check that the maximum displacement is not greater than Smax
1969
        If (debug) WRITE(*,*) "Fini before test:",Fini
1970
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1971
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
1972
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1973
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
1974

    
1975
        GradTmp=GraDTmp*FactStep
1976

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

    
1980
        IF (debug) THEN
1981
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
1982
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1983
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1984
           ELSE
1985
              DO Iat=1,Nat
1986
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1987
                      XyzGeomF(IGeom,1:3,Iat)
1988
              END DO
1989
           END IF
1990
        END IF ! matches if (debug) THEN
1991

    
1992
        SELECT CASE (COORD)
1993
        CASE ('ZMAT')
1994
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1995
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1996
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1997
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1998
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1999
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
2000
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2001
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
2002
           Idx=4
2003
           DO Iat=4,Nat
2004
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2005
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2006
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2007
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2008
              Idx=Idx+3
2009
           END DO
2010

    
2011
        CASE ('BAKER')
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,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
2016

    
2017
        CASE ('MIXED')
2018
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2019
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2020
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2021
           DO Iat=1,NCart
2022
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
2023
           END DO
2024
           Idx=3*NCart+1
2025
           IF (NCart.EQ.1) THEN
2026
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2027
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
2028
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2029
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
2030

    
2031
              Idx=7
2032
           END IF
2033
           IF (NCart.EQ.2) THEN
2034
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2035
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
2036
              Idx=9
2037
           END IF
2038

    
2039
           
2040
           DO Iat=max(NCart+1,4),Nat
2041
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2042
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2043
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2044
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2045
              Idx=Idx+3
2046
           END DO
2047
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2048

    
2049
        CASE ('HYBRID','CART')
2050
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
2051
        CASE DEFAULT
2052
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
2053
           STOP
2054
        END SELECT
2055

    
2056
        IF (debug) THEN
2057
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
2058
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2059
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
2060
           ELSE
2061
              DO Iat=1,Nat
2062
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
2063
                      XyzGeomF(IGeom,1:3,Iat)
2064
              END DO
2065
           END IF
2066
        END IF ! matches if (debug) THEN
2067

    
2068
        ! We project out the tangential part of the gradient to know if we are done
2069
        GradTmp=Grad(IGeom,:)
2070
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
2071
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
2072
        SELECT CASE (COORD)
2073
        CASE ('ZMAT','MIXED','BAKER')
2074
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2075
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
2076
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2077
           Ca=S/(sqrt(Norm)*NormGrad)
2078
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2079
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2080
        CASE ('CART','HYBRID')
2081
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2082
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
2083
           Ca=S/(sqrt(Norm)*NormGrad)
2084
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2085
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2086
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2087
        CASE DEFAULT
2088
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
2089
           STOP
2090
        END SELECT
2091
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
2092
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
2093
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
2094
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
2095

    
2096
     END DO ! matches DO IGeom=2,NGeomF-1
2097
     ! we save the old gradients
2098
     GradOld=Grad
2099
     EPathOld=Energies
2100

    
2101

    
2102

    
2103
     ! Shall we re-parameterize the path ?
2104
     ! PFL 2007/Apr/10 ->
2105
     ! We call PathCreate in all cases, it will take care of the 
2106
     ! reparameterization, as well as calculating the tangents.
2107
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
2108
     XyzGeomI=XyzGeomF
2109
     IntCoordI=IntCoordF
2110

    
2111
     Call PathCreate(IOpt)
2112
     !     END IF
2113
     ! <- PFL 2007/Apr/10
2114

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

    
2118
        ! We have the new path, we calculate its energy and gradient
2119

    
2120
        Call EgradPath(IOpt)
2121
        !IF(IOPT .EQ. 10) Then
2122
        !         Print *, 'Stopping at my checkpoint.'
2123
        !           STOP !This is my temporary checkpoint.
2124
        !ENDIF
2125

    
2126
        ! PFL 31 Mar 2008 ->
2127
        ! Taken from v3.99 but we added a flag...
2128
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
2129
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2130
        ! else it is decreased by 0.8
2131

    
2132
        if (DynMaxStep) THEN
2133
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
2134
           if (OptReac) THEN
2135
              If (Energies(1)<EPathOld(1)) Then
2136
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
2137
              ELSE
2138
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
2139
              END IF
2140
           END IF
2141

    
2142
           if (OptProd) THEN
2143
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
2144
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
2145
              ELSE
2146
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
2147
              END IF
2148
           END IF
2149

    
2150
           DO IGeom=2,NGeomF-1
2151
              If (Energies(IGeom)<EPathOld(IGeom)) Then
2152
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
2153
              ELSE
2154
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
2155
              END IF
2156
           END DO
2157
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2158
        END IF ! test on DynMaxStep
2159
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2160
        ! <- PFL 31 MAr 2008
2161
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
2162
        !  It should get updated for other cases too.
2163

    
2164
        DO IGeom=1,NGeomF
2165
           SELECT CASE (COORD)
2166
           CASE('ZMAT')
2167
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2168
              GeomTmp=IntCoordF(IGeom,:)
2169
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2170
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2171
                   OrderInv(indzmat(2,2)),Geomtmp(1)
2172
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2173
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
2174
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
2175
              Idx=4
2176
              DO Iat=4,Nat
2177
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2178
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
2179
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
2180
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
2181
                 Idx=Idx+3
2182
              END DO
2183
           CASE('BAKER')
2184
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2185
              !GeomTmp=IntCoordF(IGeom,:)
2186
           CASE('CART','HYBRID','MIXED')
2187
              WRITE(IOOUT,'(1X,I5)') Nat
2188
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2189
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2190
              DO I=1,Nat
2191
                 Iat=I
2192
                 If (renum) Iat=OrderInv(I)
2193
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
2194
              END DO
2195
           CASE DEFAULT
2196
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
2197
              STOP
2198
           END SELECT
2199
        END DO ! matches DO IGeom=1,NGeomF
2200

    
2201
        Call Write_path(IOpt)
2202
! Shall we analyze the geometries ?
2203
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
2204

    
2205
        ! We update the Hessian matrices
2206
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
2207
        ! First using the displacement from the old path
2208
        IGeom0=2
2209
        IGeomF=NGeomF-1
2210
        IF (OptReac) IGeom0=1
2211
        IF (OptProd) IGeomF=NGeomF
2212

    
2213
        IF (FCalcHess) THEN
2214
           DO IGeom=IGeom0,IGeomF
2215
              SELECT CASE (COORD)
2216
              CASE ('ZMAT','MIXED','BAKER')
2217
                 GeomTmp2=IntCoordF(IGeom,:)
2218
              CASE ('CART','HYBRID')
2219
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2220
              CASE DEFAULT
2221
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2222
                 STOP
2223
              END SELECT
2224
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
2225
           END DO
2226
        ELSE
2227
           DO IGeom=IGeom0,IGeomF
2228
              GeomTmp=GeomOld_all(IGeom,:)
2229
              SELECT CASE (COORD)
2230
              CASE ('ZMAT','MIXED','BAKER')
2231
                 GeomTmp2=IntCoordF(IGeom,:)
2232
              CASE ('CART','HYBRID')
2233
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2234
              CASE DEFAULT
2235
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2236
                 STOP
2237
              END SELECT
2238
              GeomTmp=GeomTmp2-GeomTmp
2239
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
2240
              HessTmp=Hess(IGeom,:,:)
2241
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2242
              Hess(IGeom,:,:)=HessTmp
2243
           END DO ! matches DO IGeom=IGeom0,IGeomF
2244

    
2245
           ! Second using the neighbour points
2246
           IF (HupNeighbour) THEN
2247
              ! Only one point for end points.
2248
              IF (OptReac) THEN
2249
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2250
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
2251
                 ELSE
2252
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
2253
                 END IF
2254
                 GradTmp=Grad(1,:)-Grad(2,:)
2255
                 HessTmp=Hess(1,:,:)
2256
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2257
                 Hess(1,:,:)=HessTmp
2258
              END IF
2259

    
2260
              IF (OptProd) THEN
2261
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2262
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
2263
                 ELSE
2264
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
2265
                 END IF
2266
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
2267
                 HessTmp=Hess(NGeomF,:,:)
2268
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2269
                 Hess(NGeomF,:,:)=HessTmp
2270
              END IF
2271

    
2272
              ! Two points for the rest of the path
2273
              DO IGeom=2,NGeomF-1
2274
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
2275
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
2276
                 ELSE
2277
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
2278
                 END IF
2279
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
2280
                 HessTmp=Hess(IGeom,:,:)
2281
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2282

    
2283
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2284
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
2285
                 ELSE
2286
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
2287
                 END IF
2288
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
2289

    
2290
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2291
                 Hess(IGeom,:,:)=HessTmp
2292
              END DO ! matches DO IGeom=2,NGeomF-1
2293
           END IF !matches IF (HupNeighbour) THEN
2294
        END IF ! matches IF (FCalcHess)
2295
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
2296
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
2297

    
2298
!   IF (PROG=="GAUSSIAN") THEN
2299
!      DEALLOCATE(Gauss_paste)
2300
!      DEALLOCATE(Gauss_root)
2301
!      DEALLOCATE(Gauss_comment)
2302
!      DEALLOCATE(Gauss_end)
2303
!   END IF
2304

    
2305
  DEALLOCATE(XyzGeomF, IntCoordF)
2306
  DEALLOCATE(XyzGeomI, IntCoordI)
2307
  DEALLOCATE(XyzTangent,IntTangent)
2308
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2309
  DEALLOCATE(GeomOld_all)
2310

    
2311
  if (PROG=="GAUSSIAN") THEN
2312
     ! We de-allocate the Gauss_XX lists
2313
     ! transverse the list and deallocate each node
2314
     previous => Gauss_root ! make current point to head of list
2315
     DO
2316
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2317
        current => previous%next ! make list point to next node of head
2318
        DEALLOCATE(previous) ! deallocate current head node
2319
        previous => current  ! make current point to new head
2320
     END DO
2321

    
2322
     previous => Gauss_comment ! make current point to head of list
2323
     DO
2324
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2325
        current => previous%next ! make list point to next node of head
2326
        DEALLOCATE(previous) ! deallocate current head node
2327
        previous => current  ! make current point to new head
2328
     END DO
2329

    
2330
     previous => Gauss_end ! make current point to head of list
2331
     DO
2332
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2333
        current => previous%next ! make list point to next node of head
2334
        DEALLOCATE(previous) ! deallocate current head node
2335
        previous => current  ! make current point to new head
2336
     END DO
2337

    
2338

    
2339
  END IF
2340

    
2341
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2342
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2343
  IF (COORD.EQ.'BAKER') THEN
2344
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2345
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2346
  END IF
2347

    
2348
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2349

    
2350
  Close(IOIN)
2351
  Close(IOOUT)
2352
  IF (AnaGeom) Close(IODAT)
2353
  IF (AnaGeom.AND.(OptGeom<=0)) Close(IOGPlot)
2354

    
2355
END PROGRAM PathOpt