Statistics
| Revision:

root / src / Path.f90 @ 5

History | View | Annotate | Download (77.3 kB)

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

    
220
PROGRAM PathOpt
221

    
222
  use VarTypes
223
  use Path_module
224
  use Io_module
225

    
226
  IMPLICIT NONE
227

    
228
  CHARACTER(132) :: FileIn, FileOut, Line
229

    
230
  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF
231
  INTEGER(KINT) :: IOpt
232
  INTEGER(KINT) :: Idx
233
  INTEGER(KINT) :: N3at, NFree, Nfr
234

    
235
  INTEGER(KINT) :: List(2*MaxFroz)
236

    
237
  REAL(KREAL) :: E,FactStep !E=Energy
238
  REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad
239
  REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E
240
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord
241
  REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord)
242
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord
243
  REAL(KREAL), ALLOCATABLE :: HessTmp(:,:)
244
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
245
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
246
  LOGICAL, SAVE :: allocation_flag
247

    
248
  ! Flag to see if frozen atoms are frozen (see below)
249
  LOGICAL :: FChkFrozen
250

    
251
  ! Energies for all points along the path at the previous iteration
252
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF, where it is deallocated: Prakash
253
  ! Maximum step allowed for this geometry
254
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF, where it is deallocated: Prakash 
255

    
256
! these are used to read temporary coordinates
257
  LOGICAL :: FFrozen,FCart
258

    
259
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
260
  !
261
  ! To analyse the frozen atoms
262
  !
263
  REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat
264
  REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat
265
  REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist
266
  LOGICAL, ALLOCATABLE :: ListAtFroz(:)
267
  INTEGER(KINT) :: Iat1, Iat2, Iat3
268

    
269
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
270

    
271
  LOGICAL :: Debug, Fini,Flag_Opt_Geom
272

    
273
!  INTEGER(KINT), EXTERNAL :: Iargc
274
  INTEGER(KINT), EXTERNAL ::  ConvertNumAt
275

    
276
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
277
  !
278
  ! For Debugging purposes
279
  !
280
  LOGICAL :: FCalcHess
281
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282

    
283

    
284
  INTERFACE
285
     function valid(string) result (isValid)
286
       CHARACTER(*), intent(in) :: string
287
       logical                  :: isValid
288
     END function VALID
289

    
290
     SUBROUTINE Read_Geom(Input)
291
       CHARACTER(32), INTENT(IN) :: input
292
     END SUBROUTINE READ_GEOM
293

    
294

    
295
     subroutine hupdate_all (n, dx, dgrad, HessOld)
296

    
297
       IMPLICIT NONE
298

    
299
       INTEGER, PARAMETER :: KINT=KIND(1)
300
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
301

    
302
       INTEGER(KINT), INTENT(IN) :: n
303
       real(KREAL) :: dx(n), dgrad(n), HessOld(n*n)
304
     END subroutine hupdate_all
305

    
306
     SUBROUTINE Hinvup_BFGS_new(Nfree,ds,dgrad,Hess)
307

    
308
       IMPLICIT NONE
309
       
310
       INTEGER, PARAMETER :: KINT=KIND(1)
311
       INTEGER, PARAMETER :: KREAL=KIND(1.D0)
312
       
313
       INTEGER(KINT), INTENT(IN) :: NFREE
314
       REAL(KREAL), INTENT(INOUT) :: Hess(Nfree,Nfree)
315
       REAL(KREAL), INTENT(IN) :: ds(Nfree)
316
       REAL(KREAL), INTENT(IN) :: dGrad(Nfree)
317
       
318
     END subroutine Hinvup_BFGS_new
319

    
320

    
321
  END INTERFACE
322

    
323
  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
324
       MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, &
325
       SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, &
326
       input,prog, ProgExe,RunMode, AtTypes,poscar,PathName,  &
327
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
328
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
329
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
330
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile
331

    
332
  Namelist/cartlist/list
333
  Namelist/frozenlist/list
334

    
335

    
336
  Flag_Opt_Geom = .FALSE. ! Initialized.
337

    
338
!!!!!!!!!!!!!!!!!!!!!!!!!
339
  ! 
340
  ! We print the Version info in file
341
  !
342
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
343

    
344
  ! We read the files names
345
!  SELECT CASE (iargc())
346
  SELECT CASE (command_argument_count())
347
  CASE (2)
348
     call getarg(1,FileIn)
349
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
350
     call getarg(2,FileOut)
351
     OPEN(IOOUT,FILE=FileOut,STATUS='UNKNOWN')
352
  CASE (1)
353
     call getarg(1,FileIn)
354
     IF ((FileIn=="-help").OR.(FileIn=="--help")) THEN
355
        WRITE(*,*) "Path mini-help"
356
        WRITE(*,*) "--------------"
357
        WRITE(*,*) ""
358
        WRITE(*,*) "Use: Path Input_file Output_file"
359
        WRITE(*,*) "Input_file starts with a Namelist called path"
360
        WRITE(*,*) ""
361
        WRITE(*,*) "Compulsory variables are:"
362
        WRITE(*,*) "NGeomi: Number of geometries defining the Initial path"
363
        WRITE(*,*) "NGeomf: Number of geometries defining the Final path"
364
        WRITE(*,*) "Nat   : Number of atoms"
365
        WRITE(*,*) ""
366
        WRITE(*,*) ""
367
        WRITE(*,*) "Other options"
368
        WRITE(*,*) "Input: string that indicates the type of the input geometries."
369
        WRITE(*,*) "Accepted values are: Cart (or Xmol or Xyz) or Vasp"
370
        WRITE(*,*) "Prog: string that indicates the program that will be used for energy and gradient calculations."
371
        WRITE(*,*) "      Accepted values are: Gaussian, Vasp, Mopac, Test, Siesta or Ext"
372
        WRITE(*,*) "      In case of a Gaussian or Mopac calculations, input must be set to Cart."
373
        WRITE(*,*) "      One example of a gaussian/mopac input should be added at the end of the"
374
        WRITE(*,*) "      input file.See example file Test_G03.path or Test_Mopac.path"
375
        WRITE(*,*) "      In the case of a VASP calculation, if input is set to Cart, then"
376
        WRITE(*,*) "    the preamble of a VASP calculation should be added at the end of"
377
        WRITE(*,*) "    the input file. See example file Test_VASP_cart.path"
378
        WRITE(*,*) "    In the case of a VASP calculation, one should also give value of the "
379
        WRITE(*,*) "    Runmode variable"
380
        WRITE(*,*) "Runmode: This indicates wether one should use VASP routine to calculate the energy"
381
        WRITE(*,*) "       and gradient of the whole path or not. If one wants to use VASP,"
382
        WRITE(*,*) "       Runmode must be set to PARA."
383
        WRITE(*,*) "PathName: Prefix used to save the path. Default is Path"
384
        WRITE(*,*) "Poscar: string that will be used as the prefix for the different "
385
        WRITE(*,*) "        POSCAR files in a VASP calculations. Usefull only if PathOnly=.TRUE.,"
386
        WRITE(*,*) "        not used for internal calculations."
387
        WRITE(*,*) " CalcName: Prefix for the files used for the energy and gradient calculations"
388
        WRITE(*,*) " ISuffix: Suffix for the input file."
389
        WRITE(*,*) " OSuffix: suffix for the output file."
390
        WRITE(*,*) "IGeomRef: Index of the geometry used to construct the internal coordinates."
391
        WRITE(*,*) "      Valid only for Coord=Zmat, Hybrid or Mixed"
392
        WRITE(*,*) "Fact: REAL used to define if two atoms are linked."
393
        WRITE(*,*) "      If d(A,B)<=fact*(rcov(A)+rcov(B)), then A and B are considered Linked."
394
        WRITE(*,*) "debugFile: Name of the file that indicates which subroutine should print debug info."
395
        WRITE(*,*) "Coord: System of coordinates to use. Possible choices are:"
396
        WRITE(*,*) "       - CART (or Xyz): works in cartesian"
397
        WRITE(*,*) "       - Zmat: works in internal coordinates (Zmat)"
398
        WRITE(*,*) "       - Mixed: frozen atoms, as well as atoms defined by the "
399
        WRITE(*,*) "       'cart' array(see below) are describe in CARTESIAN, whereas the"
400
        WRITE(*,*) "       others are described in Zmat" 
401
        WRITE(*,*) "       - Baker: use of Baker coordinates, also called delocalized internal coordinates"
402
        WRITE(*,*) "       - Hybrid: geometries are described in zmat, but the gradient are used in cartesian"
403
        WRITE(*,*) "Step_method: method to compute the step for optimizing a geometry; choices are:" 
404
        WRITE(*,*) "       - RFO: Rational function optimization"
405
        WRITE(*,*) "       - GDIIS: Geometry optimization using direct inversion in the iterative supspace"
406
!        WRITE(*,*) " HesUpd: Deprecated. Use Hupdate."
407
        WRITE(*,*) " HUpdate: method to update the hessian. By default, it is Murtagh-Sargent"
408
        WRITE(*,*) "      Except for geometry optimization where it is BFGS."
409
        WRITE(*,*) " GeomFile: Name of the file to print the geometries and their energy"
410
        WRITE(*,*) "      during a geometry optimization. If this variable is not given"
411
        WRITE(*,*) "      then nothing is printed."
412
        WRITE(*,*) "MaxCyc: maximum number of iterations for the path optimization"
413
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
414
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
415
        WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
416
        WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
417
        WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
418
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path."
419
        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
420
        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
421
        WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
422
        WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, "
423
        WRITE(*,*) "        it might happen that an atom moves out of the unit cell."
424
        WRITE(*,*) "        Path detects this by comparing the displacement to boxtol:"
425
        WRITE(*,*) "        if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5"
426
        WRITE(*,*) ""
427
        WRITE(*,*) "Arrays:"
428
        WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements."
429
        WRITE(*,*) "      You can modify it using, rcov(6)=0.8."
430
        WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements."
431
        WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations."
432
        WRITe(*,*) "If not given, Path will read the POTCAR file."
433
        WRITE(*,*) ""
434
        WRITE(*,*) "Flags:"
435
        WRITE(*,*) "MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
436
        WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time."
437
        WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
438
        WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
439
        WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para"
440
        WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para"
441
        WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
442
        WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
443
        WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
444
        WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path."
445
        WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
446
        WRITE(*,*) "         If True, a &frozenlist namelist containing the list of frozen atoms must be given."
447
        WRITE(*,*) "          If VASP is used, and frozen is not given"
448
        WRITE(*,*) " here, the program will use the F flags of the input geometry"
449
        WRITE(*,*) "FCart:  True if one wants to describe some atoms using cartesian coordinates. "
450
        WRITE(*,*) "        *** Only used in 'mixed' calculations. ***"
451
        WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
452
        WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
453
        WRITE(*,*) ""
454
        WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
455
        WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
456
        WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
457
        WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
458
        WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now"
459
        WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files."
460
        WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP."
461
        WRITE(*,*) ""
462
        STOP
463
     END IF
464

    
465
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
466
     IOOUT=6
467
  CASE (0)
468
     IOIN=5
469
     IOOUT=6
470
  END SELECT
471

    
472

    
473

    
474
  ! We initiliaze variables
475
  Pi=dacos(-1.0d0)
476
  Nat=1
477
  Fact=1.1
478
  IGeomRef=-1
479
  NGeomI=1
480
  NGeomF=8
481
  IReparam=5
482
  IReparamT=-1
483
  ISpline=5
484
  debug=.False.
485
  MW=.TRUE.
486
  bohr=.False.
487
  Coord='ZMAT'
488
  Step_method='RFO'
489
  DebugFile='Path.valid'
490
  PathName="Path"
491
  MaxCyc=50
492
  SMax=0.1D0
493
  SThresh=-1.
494
  GThresh=-1.
495
  FTan=0.
496
  Hinv=.TRUE.
497
  IniHUp=.FALSE.
498
  HupNeighbour=.FALSE.
499
  HesUpd="EMPTY"
500
  HUpdate="EMPTY"
501
  FFrozen=.FALSE.
502
  FCart=.FALSE.
503
  List=0
504
  renum=.TRUE.
505
  OptReac=.FALSE.
506
  OptProd=.FALSE.
507
  CalcEReac=.FALSE.
508
  CalcEProd=.FALSE.
509
  EReac=0.d0
510
  EProd=0.d0
511
  OptGeom=-1
512
  GeomFile="EMPTY"
513
  PathOnly=.FALSE.
514
  Prog='EMPTY'
515
  ProgExe='EMPTY'
516
  Runmode='SERIAL'
517
  Input='EMPTY'
518
  Poscar="POSCAR"
519
  AutoCart=.FALSE.
520
  VMD=.TRUE.
521
  WriteVASP=.FALSE.
522
  DynMaxStep=.FALSE.
523
  CalcName="EMPTY"
524
  ISuffix="EMPTY"
525
  OSuffix="EMPTY"
526
  NGintMax=10
527
  Align=.TRUE.
528

    
529
  ! Number of point used to compute the length of the path,
530
  ! and to reparameterize the path
531
  NMaxPtPath=-1 
532
  FrozTol=1e-4
533

    
534
  ! Read the namelist
535
  read (IOIN,path)
536

    
537
  debug=valid("Main").or.valid("Path")
538

    
539
  FCalcHess=valid("CalcHess")
540
  PathName=AdjustL(PathName)
541
  Coord=AdjustL(Coord)
542
  CALL UpCase(coord)
543

    
544
  Step_method=AdjustL(Step_method)
545
  CALL UpCase(Step_method)
546

    
547
  Prog=AdjustL(Prog)
548
  CALL UpCase(prog)
549

    
550
  Input=AdjustL(Input)
551
  CALL UpCase(input)
552

    
553
  Runmode=AdjustL(Runmode)
554
  Call UpCase(Runmode)
555
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"
556

    
557
  if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN
558
     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"
559
     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"
560
     CalcEReac=.FALSE.
561
     CalcEProd=.FALSE.
562
  END IF
563

    
564
  If (IReparamT<=0) IReparamT=IReparam
565

    
566
! We take care of HesUpd flag
567
! this is needed because some users might still use it
568
  if (HesUpd.NE."EMPTY") THEN
569
! We are pityless:
570
     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"
571
     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"
572
     STOP
573
  END IF
574

    
575
  IF (HUpdate=='EMPTY') THEN
576
     IF (OptGeom>=1) THEN
577
        HupDate="BFGS"
578
     ELSE
579
        HUpdate="MS"
580
     END IF
581
  END IF
582

    
583
  If (COORD.EQ.'XYZ') THEN 
584
     COORD='CART'
585
  END IF
586

    
587
  IF (COORD.EQ.'CART') THEN
588
     Renum=.FALSE.
589
     IF (OptGeom>0) THEN
590
        IGeomRef=OptGeom
591
     ELSE
592
        IGeomRef=1
593
     END IF
594
  END IF
595

    
596
  
597
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
598
  if (Prog.EQ.'2D') Prog="TEST2D"
599
  if (Prog.EQ.'TEST3') Prog="LEPS"
600

    
601

    
602
  Select Case (Prog)
603
    CASE ("EMPTY","G09","GAUSSIAN")
604
     Prog='GAUSSIAN'
605
     If (ProgExe=="EMPTY") ProgExe="g09"
606
     if (CalcName=="EMPTY") CalcName="Path"
607
     if (ISuffix=="EMPTY")  ISuffix="com"
608
     if (OSuffix=="EMPTY")  OSuffix="log"
609
    CASE ("G03")
610
     Prog='GAUSSIAN'
611
     If (ProgExe=="EMPTY") ProgExe="g03"
612
     if (CalcName=="EMPTY") CalcName="Path"
613
     if (ISuffix=="EMPTY")  ISuffix="com"
614
     if (OSuffix=="EMPTY")  OSuffix="log"
615
    CASE ("EXT")
616
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
617
     if (CalcName=="EMPTY") CalcName="Ext"
618
     if (ISuffix=="EMPTY")  ISuffix="in"
619
     if (OSuffix=="EMPTY")  OSuffix="out"
620
    CASE ('MOPAC')
621
     If (ProgExe=="EMPTY") ProgExe="mopac"
622
     if (CalcName=="EMPTY") CalcName="Path"
623
     if (ISuffix=="EMPTY")  ISuffix="mop"
624
     if (OSuffix=="EMPTY")  OSuffix="out"
625
    CASE ("SIESTA")
626
     If (ProgExe=="EMPTY") ProgExe="siesta"
627
     if (CalcName=="EMPTY") CalcName="Path"
628
     if (ISuffix=="EMPTY")  ISuffix="fdf"
629
     if (OSuffix=="EMPTY")  OSuffix="out"
630
    CASE ("VASP")
631
     If (ProgExe=="EMPTY") ProgExe="VASP"
632
    CASE ("TM","TURBOMOLE")
633
     Prog="TURBOMOLE"
634
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
635
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
636
     ProgExe=""
637
    CASE DEFAULT
638
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
639
     STOP
640
  END SELECT 
641

    
642
  if (Input.EQ.'EMPTY') THEN
643
     Select Case (Prog)
644
       CASE ("VASP")
645
          Input=Prog
646
       CASE ("CHAMFRE")
647
          Input="VASP"
648
       CASE DEFAULT
649
          Input='XYZ'
650
     END SELECT
651
    WRITE(*,*) "Input has been set to the default: ",INPUT
652
 END IF
653

    
654

    
655

    
656
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
657

    
658
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
659
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
660
     STOP
661
  END IF
662

    
663
  if (OptGeom.GE.1) THEN
664
     FPrintGeom=.FALSE.
665
     If (GeomFile/='EMPTY') THEN
666
        FPrintGeom=.TRUE.
667
     END IF
668
     IF (IGeomRef.LE.0) THEN
669
        IGeomRef=OptGeom
670
     ELSE
671
        IF (IGeomRef.NE.OptGeom) THEN
672
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
673
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
674

    
675
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
676
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
677
           STOP
678
        END IF
679
     END IF
680
  END IF
681

    
682
  ALLOCATE(Cart(Nat))
683
  Cart=0
684

    
685
  ALLOCATE(Frozen(Nat))
686
  Frozen=0
687

    
688
  IF (FCart) THEN
689
! We rewind the file to be sure to browse all of it to read the Cart namelist
690
     REWIND(IOIN)
691
     List=0
692
     READ(IOIN,cartlist)
693
     IF (COORD.NE.'CART') THEN
694
        Cart=List(1:Nat)
695
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
696
     ELSE
697
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
698
        WRITE(*,*) "I will discard the cartlist"
699
        Cart=0
700
     END IF
701
  END IF
702

    
703
  if (FFrozen) THEN
704

    
705
     REWIND(IOIN)
706
     List=0
707
     READ(IOIN,Frozenlist)
708
     Frozen=List(1:Nat)
709
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
710
  END IF
711

    
712
  If (FCart) Call Expand(Cart,NCart,Nat)
713
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
714

    
715
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
716
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
717

    
718
  if (NMaxPtPath.EQ.-1) then
719
     NMaxPtPath=min(50*NGeomF,2000)
720
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
721
  end if
722

    
723
  !  NMaxPtPath=10000
724

    
725
  ! We ensure that FTan is in [0.:1.]
726
  FTan=min(abs(FTan),1.d0)
727

    
728
! PFL 2011 Mar 14 ->
729
! Added some printing for analyses with Anapath
730
  if (debug) THEN
731
     WRITE(IOOUT,path)
732
  ELSE
733
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC
734
! and PATHNAME 
735
     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat
736
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI
737
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF
738
     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC
739
     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)
740
  END IF
741

    
742
  FTan=1.-FTan
743

    
744
  !Thresholds...
745
  if (SThresh.LE.0) SThresh=0.01d0
746
  If (GThresh.LE.0) GThresh=SThresh/4.
747
  SThreshRms=Sthresh*2./3.
748
  GThreshRms=Gthresh*2./3.
749

    
750
  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
751
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
752
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
753
  ALLOCATE(MassAt(NAt),OrderInv(Nat))
754

    
755
  ! We read the initial geometries
756
  Call Read_Geom(input)
757

    
758
  ! We convert atome names into atom mass number
759
  DO I=1,NAt
760
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
761
     Atome(I)=ConvertNumAt(AtName(I))
762
  END DO
763

    
764
!*********** HERE *************
765
!* To DO: 
766
!* deplacer le test sur les frozen dans une sous routine
767
!* deplacer la lecture de sinput dans une sous routine
768
!* creer lecture input siesta
769
!* creer egrad siesta
770
!*
771
!*************************************
772

    
773
  ! PFL 23 Sep 2008 ->
774
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
775
  ! between the first geometry and the others.
776
  IF (NFroz.GT.0) THEN
777
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
778
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
779
     ListAtFroz=.FALSE.
780

    
781
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
782
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
783
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
784

    
785
     SELECT CASE (NFroz)
786
        ! We have Frozen atoms
787
        ! PFL 28 Oct 2008 ->
788
        ! It might happen that the geometries are translated/rotated.
789
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
790
        ! three frozen atoms.
791

    
792
     CASE (3:)
793
        DO I=1,NFroz
794
           ListAtFroz(Frozen(I))=.TRUE.
795
        END DO
796
     CASE (2)
797
        IAt1=Frozen(1)
798
        IAt2=Frozen(2)
799
        ListAtFroz(Iat1)=.TRUE.
800
        ListAtFroz(Iat2)=.TRUE.
801
        x2(1)=x1(Iat2)-x1(Iat1)
802
        y2(1)=y1(Iat2)-y1(Iat1)
803
        z2(1)=z1(Iat2)-z1(Iat1)
804
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
805
        Dist=-1.
806
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
807
        DO I=1,Nat
808
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
809
           x2(2)=x1(Iat2)-x1(I)
810
           y2(2)=y1(Iat2)-y1(I)
811
           z2(2)=z1(Iat2)-z1(I)
812
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
813
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
814
           if (abs(ca)<0.9) THEN
815
              IF (Norm2>Dist) THEN
816
                 Iat3=I
817
                 Dist=Norm2
818
              END IF
819
           END IF
820
        END DO
821
        ListAtFroz(Iat3)=.TRUE.
822
        if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3
823
     CASE (1)
824
        IAt1=Frozen(1)
825
        ListAtFroz(Iat1)=.TRUE.
826
        Dist=-1.
827
        ! We scan all atoms to find one that is further from At1
828
        DO I=1,Nat
829
           if (I.EQ.Iat1) Cycle
830
           x2(1)=x1(Iat1)-x1(I)
831
           y2(1)=y1(Iat1)-y1(I)
832
           z2(1)=z1(Iat1)-z1(I)
833
           Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
834
           IF (Norm1>Dist) THEN
835
              Iat2=I
836
              Dist=Norm1
837
           END IF
838
        END DO
839
        IAt2=Frozen(2)
840
        ListAtFroz(Iat2)=.TRUE.
841
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
842
        x2(1)=x1(Iat2)-x1(Iat1)
843
        y2(1)=y1(Iat2)-y1(Iat1)
844
        z2(1)=z1(Iat2)-z1(Iat1)
845
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
846
        Dist=-1.
847
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
848
        DO I=1,Nat
849
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
850
           x2(2)=x1(Iat2)-x1(I)
851
           y2(2)=y1(Iat2)-y1(I)
852
           z2(2)=z1(Iat2)-z1(I)
853
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
854
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
855
           if (abs(ca)<0.9) THEN
856
              IF (Norm2>Dist) THEN
857
                 Iat3=I
858
                 Dist=Norm2
859
              END IF
860
           END IF
861
        END DO
862
        ListAtFroz(Iat3)=.TRUE.
863
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3
864
     END SELECT
865

    
866
     DO I=2,NGeomI
867
        ! First, we align the Ith geometry on the first one, using only
868
        ! the positions of the "frozen" atoms. (see above: it might be that
869
        ! ListAtFroz contains non frozen atoms usd to align the geometries
870
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
871
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
872
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
873
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
874

    
875

    
876
        FChkFrozen=.FALSE.
877
        DO J=1,NFroz
878
           IAt=Frozen(J)
879
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
880
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
881
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
882
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
883
              FChkFrozen=.TRUE.
884
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
885
              WRITE(*,*) "from geometry 1 to geometry ",I
886
           END IF
887
        END DO
888
     END DO
889
     IF (FChkFrozen) THEN
890
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
891
        STOP
892
     END IF
893
  END IF
894

    
895

    
896
  N3at=3*Nat
897
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
898

    
899

    
900
  Call ReadInput
901

    
902
  IF (COORD=="MIXED") Call TestCart(AutoCart)
903

    
904
  SELECT CASE(NFroz)
905
  CASE (2)
906
     IntFroz=1
907
  CASE (3)
908
     IntFroz=3
909
  CASE (4:)
910
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
911
  END SELECT
912

    
913
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
914
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
915
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
916

    
917
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
918
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
919
     COORD="ZMAT"
920
  END IF
921

    
922
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
923

    
924
  SELECT CASE(COORD)
925
  CASE ('ZMAT')
926
     NCoord=Nfree
927
     ALLOCATE(dzdc(3,nat,3,nat))
928
  CASE ('MIXED')
929
     SELECT CASE (NCart+NFroz)
930
     CASE (1) 
931
        NCoord=N3At-3
932
     CASE (2)
933
        NCoord=N3At-1
934
     CASE (3:)
935
        NCoord=N3At
936
     END SELECT
937
     ALLOCATE(dzdc(3,nat,3,nat))
938
  CASE ('BAKER')
939
     Symmetry_elimination=0
940
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
941
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
942
     ALLOCATE(BTransInv(NCoord,3*Nat))
943
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
944
  CASE ('HYBRID')
945
     NCoord=N3at
946
  CASE ('CART')
947
     NCoord=N3at
948
  END SELECT
949

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

    
952
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
953
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
954
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
955
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
956
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
957
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
958
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
959
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
960

    
961
  ZeroVec=0.d0
962
  DO I =1, NGeomF
963
     IntTangent(I,:)=0.d0
964
  END DO
965

    
966
  if (debug) THEN
967
     Print *, 'L886, IntTangent(1,:)='
968
     Print *, IntTangent(1,:)
969
  END IF
970

    
971
  if (.NOT.OptReac) Energies(1)=Ereac
972
  if (.NOT.OptProd) Energies(NGeomF)=EProd
973
  MaxStep=SMax
974

    
975
  ! We initialize the mass array
976
  if (MW) THEN
977
     WRITE(*,*) 'Working in MW coordinates'
978
     DO I=1,Nat
979
        massAt(I)=Mass(Atome(I))
980
     END DO
981
  ELSE
982
     DO I=1,Nat
983
        MassAt(I)=1.d0
984
     END DO
985
  END IF
986

    
987
  WRITE(*,*) "Prog=",TRIM(PROG)
988

    
989
  ALLOCATE(FrozAtoms(Nat))
990

    
991
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
992
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
993
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
994
  IF (NFroz.EQ.0) THEN
995
     FrozAtoms=.TRUE.
996
  ELSE
997
     I=1
998
     NFr=0
999
     FrozAtoms=.False.
1000
     DO WHILE (Frozen(I).NE.0)
1001
        IF (Frozen(I).LT.0) THEN
1002
           DO J=Frozen(I-1),abs(Frozen(I))
1003
              IF (.NOT.FrozAtoms(J)) THEN
1004
                 NFr=NFr+1
1005
                 FrozAtoms(J)=.TRUE.
1006
              END IF
1007
           END DO
1008
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1009
           FrozAtoms(Frozen(I))=.TRUE.
1010
           NFr=NFr+1
1011
        END IF
1012
        I=I+1
1013
     END DO
1014
     IF (NFr.NE.NFroz) THEN
1015
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1016
        STOP
1017
     END IF
1018
  END IF
1019

    
1020
!  if (debug) THEN
1021
     !Print *, 'L968, IntTangent(1,:)='
1022
     !Print *, IntTangent(1,:)
1023
!  END IF
1024

    
1025
  ! We have read everything we needed... time to work :-)
1026
  IOpt=0
1027
  FirstTimePathCreate = .TRUE. ! Initialized.
1028
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1029
  FirstTimePathCreate = .FALSE.
1030

    
1031
  if (debug) THEN
1032
     Print *, 'L980, IntTangent(1,:)='
1033
     Print *, IntTangent(1,:)
1034
  END IF
1035

    
1036
  ! PFL 30 Mar 2008 ->
1037
  ! The following is now allocated in Calc_Baker. This is more logical to me
1038
  !   IF (COORD .EQ. "BAKER") Then
1039
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1040
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1041
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1042
  !   END IF
1043
  ! <- PFL 30 mar 2008
1044

    
1045
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1046
  !
1047
  ! For Debugging purposes
1048
  !
1049
  ! OptGeom is the index of the geometry which is to be optimized.
1050
  IF (Optgeom.GT.0) THEN
1051
     Flag_Opt_Geom=.TRUE.
1052
     SELECT CASE(Coord)
1053
     CASE ('CART','HYBRID')
1054
!!! CAUTION : PFL 29.AUG.2008 ->
1055
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1056
        ! This should be made  consistent !!!!!!!!!!!!!!!
1057
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1058
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1059
!!! <- CAUTION : PFL 29.AUG.2008
1060
     CASE ('ZMAT','MIXED')
1061
        !GeomTmp=IntCoordF(OptGeom,:)
1062
        GeomTmp=IntCoordI(OptGeom,:)
1063
     CASE ('BAKER')
1064
        !GeomTmp=IntCoordF(OptGeom,:)
1065
        GeomTmp=IntCoordI(OptGeom,:)
1066
     CASE DEFAULT
1067
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1068
        STOP
1069
     END SELECT
1070

    
1071
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1072
     Flag_Opt_Geom = .TRUE.
1073
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1074

    
1075
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1076
     WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1077
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1078
     STOP
1079
  END IF ! IF (Optgeom.GT.0) THEN
1080

    
1081
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1082
  ! End of GeomOpt
1083
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1084

    
1085
  IF (PathOnly) THEN
1086
     WRITE(*,*) "PathOnly=.T. , Stop here"
1087
     Call Write_path(-1)
1088
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1089
     STOP
1090
  END IF
1091

    
1092
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1093

    
1094
  DEALLOCATE(XyzGeomI,IntCoordI)
1095
  NGeomI=NGeomF
1096
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1097

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

    
1103
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1104
     DO IGeom=1,NGeomF
1105
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1106
        DO I=1,3*Nat-6
1107
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1108
        END DO
1109
     END DO
1110
  END IF
1111

    
1112
  ! To calculate B^T^-1 for all extrapolated final geometries:
1113
  ! PFL 18 Jan 2008 ->
1114
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1115
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1116
  ! <-- PFL 18 Jan 2008
1117

    
1118
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1119
     DO IGeom=1,NGeomF
1120
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
1121
        DO I=1,3*Nat-6
1122
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1123
        END DO
1124
     END DO
1125
  END IF
1126

    
1127
  IOpt=0
1128
  Call EgradPath(IOpt)
1129

    
1130
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1131
  Call Write_path(IOpt)
1132
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1133

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

    
1138
  ! By default, Hess=Id
1139
  Hess=0.
1140
  IF(Coord .EQ. "ZMAT") Then
1141
     ! We use the fact that we know that approximate good hessian values
1142
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1143
     DO IGeom=1,NGeomF
1144
        Hess(IGeom,1,1)=0.5d0
1145
        Hess(IGeom,2,2)=0.5d0
1146
        Hess(IGeom,3,3)=0.1d0
1147
        DO J=1,Nat-3
1148
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1149
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1150
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1151
        END DO
1152
        IF (HInv) THEN
1153
           DO I=1,NCoord
1154
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1155
           END DO
1156
        END IF
1157
     END DO
1158
  ELSE
1159
     DO I=1,NCoord
1160
        DO J=1,NGeomF
1161
           Hess(J,I,I)=1.d0
1162
        END DO
1163
     END DO
1164
  END IF
1165

    
1166
  IF (COORD .EQ. "BAKER") THEN
1167
     ! UMat(NPrim,NCoord)
1168
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1169
     Hprim=0.d0
1170
     ScanCoord=>Coordinate
1171
     I=0
1172
     DO WHILE (Associated(ScanCoord%next))
1173
        I=I+1
1174
        SELECT CASE (ScanCoord%Type)
1175
        CASE ('BOND')
1176
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1177
           Hprim(I,I) = 0.5d0
1178
           !END DO
1179
        CASE ('ANGLE')
1180
           Hprim(I,I) = 0.2d0
1181
        CASE ('DIHEDRAL')
1182
           Hprim(I,I) = 0.1d0
1183
        END SELECT
1184
        ScanCoord => ScanCoord%next
1185
     END DO
1186

    
1187
     ! HprimU = H^prim U:
1188
     HprimU=0.d0
1189
     DO I=1, 3*Nat-6
1190
        DO J=1,NPrim
1191
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1192
        END DO
1193
     END DO
1194

    
1195
     ! Hess = U^T HprimU = U^T H^prim U:
1196
     Hess=0.d0
1197
     DO I=1, 3*Nat-6
1198
        DO J=1,NPrim
1199
           ! UMat^T is needed below.
1200
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1201
        END DO
1202
     END DO
1203
     DO K=2,NGeomF
1204
        Hess(K,:,:)=Hess(1,:,:)
1205
     END DO
1206
     DEALLOCATE(Hprim,HprimU)
1207
  end if !  ends IF COORD=BAKER
1208
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1209
  ALLOCATE(HessTmp(NCoord,NCoord))
1210

    
1211
  IF (IniHUp) THEN
1212
     IF (FCalcHess) THEN
1213
        ! The numerical calculation of the Hessian has been switched on
1214
        DO IGeom=1,NGeomF
1215
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1216
              GeomTmp=IntCoordF(IGeom,:)
1217
           ELSE
1218
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1219
           END IF
1220
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1221
        END DO
1222
     ELSE
1223
        ! We generate a forward hessian and a backward one and we mix them.
1224
        ! First, the forward way:
1225
        DO IGeom=2,NGeomF-1
1226
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1227
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1228
           ELSE
1229
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1230
           END IF
1231
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1232
           HessTmp=Hess(IGeom-1,:,:)
1233
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1234
           Hess(IGeom,:,:)=HessTmp
1235
        END DO
1236

    
1237
        ! Now, the backward way:
1238
        HessTmp=Hess(NGeomF,:,:)
1239
        DO IGeom=NGeomF-1,2,-1
1240
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1241
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1242
           ELSE
1243
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1244
           END IF
1245
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1246
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1247
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1248
           ca=cos(alpha)
1249
           sa=sin(alpha)
1250
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1251
        END DO
1252
     END IF ! matches IF (FCalcHess)
1253
  END IF ! matches IF (IniHUp) THEN
1254

    
1255
  ! For debug purposes, we diagonalize the Hessian matrices
1256
  IF (debug) THEN
1257
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1258
     DO I=1,NGeomF
1259
        WRITE(*,*) "DBG Main - Point ",I
1260
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1261
        DO J=1,NCoord
1262
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1263
        END DO
1264
     END DO
1265
  END IF ! matches IF (debug) THEN
1266

    
1267
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1268
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1269
  !
1270
  ! Main LOOP : Optimization of the Path
1271
  !
1272
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1273
  if (debug)   Print *, 'NGeomF=', NGeomF
1274
  allocation_flag = .TRUE.
1275

    
1276
  Fini=.FALSE.
1277

    
1278
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1279
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1280
     IOpt=IOpt+1
1281

    
1282
     SELECT CASE (COORD)
1283
     CASE ('ZMAT','MIXED','BAKER')
1284
        GeomOld_all=IntCoordF
1285
     CASE ('CART','HYBRID')
1286
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1287
     CASE DEFAULT
1288
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1289
        STOP
1290
     END SELECT
1291

    
1292
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1293
          valid("printtangent")) THEN
1294
        WRITE(*,*) "Visualization of tangents"
1295
        Idx=min(12,NCoord)
1296
        DO I=1,NGeomF
1297
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1298
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1299
           WRITE(*,*) 
1300
        END DO
1301
        WRITE(*,*) "END of tangents"
1302
     END IF
1303

    
1304
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1305
          valid("printgrad")) THEN
1306
        WRITE(*,*) "Visualization of gradients"
1307
        Idx=min(12,NCoord)
1308
        DO I=1,NGeomF
1309
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1310
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1311
           WRITe(*,*) 
1312
        END DO
1313
        WRITE(*,*) "END of gradients"
1314
     END IF
1315

    
1316
     Fini=.TRUE.
1317
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1318
        IGeom=1
1319
        if (debug) WRITE(*,*) '**************************************'
1320
        if (debug) WRITE(*,*) 'Optimizing reactant'
1321
        if (debug) WRITE(*,*) '**************************************'
1322
        SELECT CASE (COORD)
1323
        CASE ('ZMAT','MIXED','BAKER')
1324
           SELECT CASE (Step_method)
1325
           CASE ('RFO')
1326
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
1327
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1328
              Print *, GradTmp
1329
           CASE ('GDIIS')
1330
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1331
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1332
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1333
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1334
              Print *, GradTmp
1335
           CASE ('GEDIIS')
1336
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1337
              ! Energies are updated in EgradPath.f90
1338
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1339
                   NCoord,allocation_flag,ZeroVec)
1340
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1341
              !GradTmp = GeomTmp - IntCoordF(1,:)
1342
              !Print *, 'Old Geometry:'
1343
              !Print *, IntCoordF(1,:)
1344
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1345
              Print *, GradTmp
1346
              !Print *, 'GeomTmp:'
1347
              !Print *, GeomTmp
1348
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1349
           CASE DEFAULT
1350
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1351
              STOP
1352
           END SELECT
1353
        CASE ('HYBRID','CART')
1354
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1355
        CASE DEFAULT
1356
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1357
           STOP
1358
        END SELECT
1359

    
1360
        IF (debug) THEN
1361
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1362
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1363
        END IF
1364

    
1365
        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
1366
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1367
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
1368
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1369
        OptReac=(NormStep.GT.SThresh)
1370
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
1371
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
1372
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1373
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
1374
        !Print *, 'Grad(1,:):'
1375
        !Print *, Grad(1,:)
1376
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac
1377

    
1378
        GradTmp=GradTmp*FactStep
1379

    
1380
        ! we take care that frozen atoms do not move.
1381
        IF (NFroz.NE.0) THEN
1382
           SELECT CASE (COORD)
1383
           CASE ('ZMAT','MIXED')
1384
              GradTmp(1:IntFroz)=0.d0
1385
           CASE ('CART','HYBRID')
1386
              DO I=1,Nat
1387
                 IF (FrozAtoms(I)) THEN
1388
                    Iat=Order(I)
1389
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1390
                 END IF
1391
              END DO
1392
           CASE('BAKER')
1393
              GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..)
1394
           CASE DEFAULT
1395
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP'
1396
              STOP
1397
           END SELECT
1398
        END IF ! matches IF (NFroz.NE.0) THEN
1399

    
1400
        IF (debug) THEN
1401
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1402
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1403
        END IF
1404

    
1405
        IF (debug) THEN
1406
           !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1'
1407
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1408
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1409
           ELSE
1410
              DO Iat=1,Nat
1411
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1412
                      XyzGeomF(IGeom,1:3,Iat)
1413
              END DO
1414
           END IF
1415
        END IF ! matches IF (debug) THEN
1416

    
1417
        SELECT CASE (COORD)
1418
        CASE ('ZMAT','MIXED','BAKER')
1419
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1420
        CASE ('HYBRID','CART')
1421
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1422
        CASE DEFAULT
1423
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1424
           STOP
1425
        END SELECT
1426

    
1427
        IF (debug) THEN
1428
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1429
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1430
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1431
           ELSE
1432
              DO Iat=1,Nat
1433

    
1434
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1435
                 ! It was with indzmat. IndZmat cannot appear here...
1436
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1437
                      XyzGeomF(IGeom,1:3,Iat)
1438
              END DO
1439
           END IF
1440
        END IF ! matches IF (debug) THEN
1441

    
1442
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1443
     ELSE ! Optreac
1444
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1445
     END IF ! matches IF (OptReac) THEN
1446

    
1447
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1448
        IGeom=NGeomF
1449
        if (debug) WRITE(*,*) '******************'
1450
        if (debug) WRITE(*,*) 'Optimizing product'
1451
        if (debug) WRITE(*,*) '******************'
1452
        SELECT CASE (COORD)
1453
        CASE ('ZMAT','MIXED','BAKER')
1454
           Print *, 'Optimizing product'
1455
           SELECT CASE (Step_method)
1456
           CASE ('RFO')
1457
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1458
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1459
              Print *, GradTmp
1460
           CASE ('GDIIS')
1461
              ! GradTmp becomes "step" below and has INTENT(OUT):
1462
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1463
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1464
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1465
              Print *, GradTmp
1466
           CASE ('GEDIIS')
1467
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1468
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1469
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1470
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1471
              Print *, GradTmp
1472
           CASE DEFAULT
1473
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1474
              STOP
1475
           END SELECT
1476
        CASE ('HYBRID','CART')
1477
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1478
        CASE DEFAULT
1479
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1480
           STOP
1481
        END SELECT
1482

    
1483
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1484
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1485
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1486
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1487
        OptProd=.NOT.(NormStep.LE.SThresh)
1488
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1489
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1490
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1491
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1492

    
1493
        GradTmp=GradTmp*FactStep
1494

    
1495
        ! we take care that frozen atoms do not move
1496
        IF (NFroz.NE.0) THEN
1497
           SELECT CASE (COORD)
1498
           CASE ('ZMAT','MIXED','BAKER')
1499
              GradTmp(1:IntFroz)=0.d0
1500
           CASE ('CART','HYBRID')
1501
              DO I=1,Nat
1502
                 IF (FrozAtoms(I)) THEN
1503
                    Iat=Order(I)
1504
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1505
                 END IF
1506
              END DO
1507
           CASE DEFAULT
1508
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1509
              STOP
1510
           END SELECT
1511
        END IF ! matches IF (NFroz.NE.0) THEN
1512

    
1513
        IF (debug) THEN
1514
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
1515
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1516
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1517
           ELSE
1518
              DO Iat=1,Nat
1519
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1520
                      XyzGeomF(IGeom,1:3,Iat)
1521
              END DO
1522
           END IF
1523
        END IF
1524

    
1525
        SELECT CASE (COORD)
1526
        CASE ('ZMAT','MIXED','BAKER')
1527
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1528
        CASE ('HYBRID','CART')
1529
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1530
        CASE DEFAULT
1531
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1532
           STOP
1533
        END SELECT
1534

    
1535
        IF (debug) THEN
1536
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1537
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1538
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1539
           ELSE
1540
              DO Iat=1,Nat
1541
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1542
                      XyzGeomF(IGeom,1:3,Iat)
1543
              END DO
1544
           END IF
1545
        END IF
1546
     ELSE ! Optprod
1547
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1548
     END IF !matches IF (OptProd) THEN 
1549

    
1550
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1551

    
1552
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1553
     !
1554
     !  Optimizing other geometries
1555
     !
1556
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1557

    
1558

    
1559

    
1560
     DO IGeom=2,NGeomF-1 ! matches at L1556
1561
        if (debug) WRITE(*,*) '****************************'
1562
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1563
        if (debug) WRITE(*,*) '****************************'
1564

    
1565
        SELECT CASE (COORD)
1566
        CASE ('ZMAT','BAKER','MIXED')
1567
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
1568
           ! PFL 6 Apr 2008 ->
1569
           ! Special case, if FTan=0. we do not consider Tangent at all.
1570
           ! To DO: add the full treatment of FTan
1571
           if (FTan==0.) GradTmp2=ZeroVec
1572
           ! <- PFL 6 Apr 2008
1573
           if (debug) THEN
1574
              Print *, 'L1500, IntTangent(',IGeom,',:)='
1575
              Print *, IntTangent(IGeom,:)
1576
           END IF
1577
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
1578
           SELECT CASE (Step_method)
1579
           CASE ('RFO')
1580
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1581
           CASE ('GDIIS')
1582
              !The following has no effect at IOpt==1
1583
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
1584
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&        
1585
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1586
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
1587
              Print *, GradTmp
1588
           CASE ('GEDIIS')
1589
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1590
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
1591
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1592
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
1593
              Print *, GradTmp
1594
           CASE DEFAULT
1595
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
1596
              STOP
1597
           END SELECT
1598
        CASE ('HYBRID','CART')
1599
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
1600
           ! so we change it
1601
           DO I=1,Nat
1602
              DO J=1,3
1603
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
1604
              END DO
1605
           END DO
1606
           !           GradTmp2=XyzTangent(IGeom,:)
1607
           ! PFL 6 Apr 2008 ->
1608
           ! Special case, if FTan=0. we do not consider Tangent at all.
1609
           ! To DO: add the full treatment of FTan
1610
           if (FTan==0.) GradTmp2=ZeroVec
1611
           ! <- PFL 6 Apr 2008
1612

    
1613
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1614
        CASE DEFAULT
1615
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
1616
           STOP
1617
        END SELECT
1618

    
1619
        ! we take care that frozen atoms do not move
1620
        IF (NFroz.NE.0) THEN
1621
           SELECT CASE (COORD)
1622
           CASE ('ZMAT','MIXED','BAKER')
1623
              IF (debug) THEN 
1624
                 WRITE(*,*) 'Step computed. Coord=',Coord
1625
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
1626
              END IF
1627
              GradTmp(1:IntFroz)=0.d0
1628
           CASE ('CART','HYBRID')
1629
              if (debug) THEN
1630
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
1631
                 DO I=1,Nat
1632
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
1633
                 END DO
1634
              END IF
1635
              DO I=1,Nat
1636
                 IF (FrozAtoms(I)) THEN
1637
                    if (debug) THEN
1638
                       write(*,*) 'Step Computed. Coord=',Coord
1639
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
1640
                    END IF
1641
                    Iat=Order(I)
1642
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1643
                 END IF
1644
              END DO
1645

    
1646
                 if (debug) THEN
1647
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
1648
                    DO I=1,Nat
1649
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
1650
                    END DO
1651
                 END IF
1652

    
1653
           CASE DEFAULT
1654
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
1655
              STOP
1656
           END SELECT
1657
        END IF !matches IF (NFroz.NE.0) THEN
1658

    
1659
        IF (debug) THEN
1660
           SELECT CASE (COORD)
1661
           CASE ('ZMAT','MIXED','BAKER')
1662
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
1663
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1664
           CASE ('CART','HYBRID')
1665
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
1666
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1667
           CASE DEFAULT
1668
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
1669
              STOP
1670
           END SELECT
1671
        END IF ! matches if (debug) THEN
1672

    
1673
        ! we check that the maximum displacement is not greater than Smax
1674
        If (debug) WRITE(*,*) "Fini before test:",Fini
1675
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1676
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
1677
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1678
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
1679

    
1680
        GradTmp=GraDTmp*FactStep
1681

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

    
1685
        IF (debug) THEN
1686
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
1687
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1688
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1689
           ELSE
1690
              DO Iat=1,Nat
1691
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1692
                      XyzGeomF(IGeom,1:3,Iat)
1693
              END DO
1694
           END IF
1695
        END IF ! matches if (debug) THEN
1696

    
1697
        SELECT CASE (COORD)
1698
        CASE ('ZMAT')
1699
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1700
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1701
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1702
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1703
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1704
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
1705
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1706
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
1707
           Idx=4
1708
           DO Iat=4,Nat
1709
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1710
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
1711
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
1712
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
1713
              Idx=Idx+3
1714
           END DO
1715

    
1716
        CASE ('BAKER')
1717
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1718
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1719
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1720
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
1721

    
1722
        CASE ('MIXED')
1723
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1724
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1725
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1726
           DO Iat=1,NCart
1727
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
1728
           END DO
1729
           Idx=3*NCart+1
1730
           IF (NCart.EQ.1) THEN
1731
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1732
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
1733
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1734
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
1735

    
1736
              Idx=7
1737
           END IF
1738
           IF (NCart.EQ.2) THEN
1739
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1740
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
1741
              Idx=9
1742
           END IF
1743

    
1744
           
1745
           DO Iat=max(NCart+1,4),Nat
1746
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1747
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
1748
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
1749
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
1750
              Idx=Idx+3
1751
           END DO
1752
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1753

    
1754
        CASE ('HYBRID','CART')
1755
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1756
        CASE DEFAULT
1757
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
1758
           STOP
1759
        END SELECT
1760

    
1761
        IF (debug) THEN
1762
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
1763
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1764
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1765
           ELSE
1766
              DO Iat=1,Nat
1767
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1768
                      XyzGeomF(IGeom,1:3,Iat)
1769
              END DO
1770
           END IF
1771
        END IF ! matches if (debug) THEN
1772

    
1773
        ! We project out the tangential part of the gradient to know if we are done
1774
        GradTmp=Grad(IGeom,:)
1775
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
1776
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
1777
        SELECT CASE (COORD)
1778
        CASE ('ZMAT','MIXED','BAKER')
1779
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
1780
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
1781
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
1782
           Ca=S/(sqrt(Norm)*NormGrad)
1783
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
1784
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
1785
        CASE ('CART','HYBRID')
1786
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
1787
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
1788
           Ca=S/(sqrt(Norm)*NormGrad)
1789
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
1790
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
1791
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
1792
        CASE DEFAULT
1793
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
1794
           STOP
1795
        END SELECT
1796
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
1797
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
1798
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1799
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
1800

    
1801
     END DO ! matches DO IGeom=2,NGeomF-1
1802
     ! we save the old gradients
1803
     GradOld=Grad
1804
     EPathOld=Energies
1805

    
1806
     ! Shall we re-parameterize the path ?
1807
     ! PFL 2007/Apr/10 ->
1808
     ! We call PathCreate in all cases, it will take care of the 
1809
     ! reparameterization, as well as calculating the tangents.
1810
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
1811
     XyzGeomI=XyzGeomF
1812
     IntCoordI=IntCoordF
1813

    
1814
     Call PathCreate(IOpt)
1815
     !     END IF
1816
     ! <- PFL 2007/Apr/10
1817

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

    
1821
        ! We have the new path, we calculate its energy and gradient
1822

    
1823
        Call EgradPath(IOpt)
1824
        !IF(IOPT .EQ. 10) Then
1825
        !         Print *, 'Stopping at my checkpoint.'
1826
        !           STOP !This is my temporary checkpoint.
1827
        !ENDIF
1828

    
1829
        ! PFL 31 Mar 2008 ->
1830
        ! Taken from v3.99 but we added a flag...
1831
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
1832
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
1833
        ! else it is decreased by 0.8
1834

    
1835
        if (DynMaxStep) THEN
1836
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
1837
           if (OptReac) THEN
1838
              If (Energies(1)<EPathOld(1)) Then
1839
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
1840
              ELSE
1841
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
1842
              END IF
1843
           END IF
1844

    
1845
           if (OptProd) THEN
1846
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
1847
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
1848
              ELSE
1849
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
1850
              END IF
1851
           END IF
1852

    
1853
           DO IGeom=2,NGeomF-1
1854
              If (Energies(IGeom)<EPathOld(IGeom)) Then
1855
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
1856
              ELSE
1857
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
1858
              END IF
1859
           END DO
1860
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
1861
        END IF ! test on DynMaxStep
1862
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
1863
        ! <- PFL 31 MAr 2008
1864
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
1865
        !  It should get updated for other cases too.
1866

    
1867
        DO IGeom=1,NGeomF
1868
           SELECT CASE (COORD)
1869
           CASE('ZMAT')
1870
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1871
              GeomTmp=IntCoordF(IGeom,:)
1872
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1873
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1874
                   OrderInv(indzmat(2,2)),Geomtmp(1)
1875
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1876
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
1877
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
1878
              Idx=4
1879
              DO Iat=4,Nat
1880
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1881
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
1882
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
1883
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
1884
                 Idx=Idx+3
1885
              END DO
1886
           CASE('BAKER')
1887
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1888
              !GeomTmp=IntCoordF(IGeom,:)
1889
           CASE('CART','HYBRID','MIXED')
1890
              WRITE(IOOUT,'(1X,I5)') Nat
1891
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1892
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1893
              DO I=1,Nat
1894
                 Iat=I
1895
                 If (renum) Iat=OrderInv(I)
1896
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
1897
              END DO
1898
           CASE DEFAULT
1899
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
1900
              STOP
1901
           END SELECT
1902
        END DO ! matches DO IGeom=1,NGeomF
1903

    
1904
        Call Write_path(IOpt)
1905

    
1906
        ! We update the Hessian matrices
1907
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
1908
        ! First using the displacement from the old path
1909
        IGeom0=2
1910
        IGeomF=NGeomF-1
1911
        IF (OptReac) IGeom0=1
1912
        IF (OptProd) IGeomF=NGeomF
1913

    
1914
        IF (FCalcHess) THEN
1915
           DO IGeom=IGeom0,IGeomF
1916
              SELECT CASE (COORD)
1917
              CASE ('ZMAT','MIXED','BAKER')
1918
                 GeomTmp2=IntCoordF(IGeom,:)
1919
              CASE ('CART','HYBRID')
1920
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1921
              CASE DEFAULT
1922
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
1923
                 STOP
1924
              END SELECT
1925
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
1926
           END DO
1927
        ELSE
1928
           DO IGeom=IGeom0,IGeomF
1929
              GeomTmp=GeomOld_all(IGeom,:)
1930
              SELECT CASE (COORD)
1931
              CASE ('ZMAT','MIXED','BAKER')
1932
                 GeomTmp2=IntCoordF(IGeom,:)
1933
              CASE ('CART','HYBRID')
1934
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1935
              CASE DEFAULT
1936
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
1937
                 STOP
1938
              END SELECT
1939
              GeomTmp=GeomTmp2-GeomTmp
1940
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
1941
              HessTmp=Hess(IGeom,:,:)
1942
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1943
              Hess(IGeom,:,:)=HessTmp
1944
           END DO ! matches DO IGeom=IGeom0,IGeomF
1945

    
1946
           ! Second using the neighbour points
1947
           IF (HupNeighbour) THEN
1948
              ! Only one point for end points.
1949
              IF (OptReac) THEN
1950
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1951
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
1952
                 ELSE
1953
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
1954
                 END IF
1955
                 GradTmp=Grad(1,:)-Grad(2,:)
1956
                 HessTmp=Hess(1,:,:)
1957
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1958
                 Hess(1,:,:)=HessTmp
1959
              END IF
1960

    
1961
              IF (OptProd) THEN
1962
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1963
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
1964
                 ELSE
1965
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
1966
                 END IF
1967
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
1968
                 HessTmp=Hess(NGeomF,:,:)
1969
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
1970
                 Hess(NGeomF,:,:)=HessTmp
1971
              END IF
1972

    
1973
              ! Two points for the rest of the path
1974
              DO IGeom=2,NGeomF-1
1975
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
1976
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1977
                 ELSE
1978
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1979
                 END IF
1980
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1981
                 HessTmp=Hess(IGeom,:,:)
1982
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
1983

    
1984
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1985
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1986
                 ELSE
1987
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1988
                 END IF
1989
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
1990

    
1991
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
1992
                 Hess(IGeom,:,:)=HessTmp
1993
              END DO ! matches DO IGeom=2,NGeomF-1
1994
           END IF !matches IF (HupNeighbour) THEN
1995
        END IF ! matches IF (FCalcHess)
1996
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
1997
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
1998

    
1999
!   IF (PROG=="GAUSSIAN") THEN
2000
!      DEALLOCATE(Gauss_paste)
2001
!      DEALLOCATE(Gauss_root)
2002
!      DEALLOCATE(Gauss_comment)
2003
!      DEALLOCATE(Gauss_end)
2004
!   END IF
2005

    
2006
  DEALLOCATE(XyzGeomF, IntCoordF)
2007
  DEALLOCATE(XyzGeomI, IntCoordI)
2008
  DEALLOCATE(XyzTangent,IntTangent)
2009
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2010
  DEALLOCATE(GeomOld_all)
2011

    
2012
  if (PROG=="GAUSSIAN") THEN
2013
     ! We de-allocate the Gauss_XX lists
2014
     ! transverse the list and deallocate each node
2015
     previous => Gauss_root ! make current point to head of list
2016
     DO
2017
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2018
        current => previous%next ! make list point to next node of head
2019
        DEALLOCATE(previous) ! deallocate current head node
2020
        previous => current  ! make current point to new head
2021
     END DO
2022

    
2023
     previous => Gauss_comment ! make current point to head of list
2024
     DO
2025
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2026
        current => previous%next ! make list point to next node of head
2027
        DEALLOCATE(previous) ! deallocate current head node
2028
        previous => current  ! make current point to new head
2029
     END DO
2030

    
2031
     previous => Gauss_end ! make current point to head of list
2032
     DO
2033
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2034
        current => previous%next ! make list point to next node of head
2035
        DEALLOCATE(previous) ! deallocate current head node
2036
        previous => current  ! make current point to new head
2037
     END DO
2038

    
2039

    
2040
  END IF
2041

    
2042
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2043
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2044
  IF (COORD.EQ.'BAKER') THEN
2045
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2046
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2047
  END IF
2048

    
2049
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2050

    
2051
END PROGRAM PathOpt