Statistiques
| Révision :

root / src / Path.f90 @ 7

Historique | Voir | Annoter | Télécharger (79,55 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
  END INTERFACE
335

    
336
  Namelist/path/fact,r_cov,nat,ngeomi,ngeomf,debugFile, &
337
       MW,mass,Ffrozen, renum, OptReac,OptProd,Coord,Step_method,MaxCyc, &
338
       SMax, SThresh, GThresh,IReparam, IReparamT,Hinv, ISpline, PathOnly,Fcart, &
339
       input,prog, ProgExe,RunMode, AtTypes,poscar,PathName,  &
340
       OptGeom,IniHup, HupNeighbour, hesUpd, HUpdate, &
341
       FTan, EReac, EProd, IGeomRef, Vmd, AutoCart, DynMaxStep, &
342
       NMaxPtPath, FrozTol, BoxTol,CalcName,ISuffix,OSuffix, WriteVasp, &
343
       NGintMax, Align, CalcEReac,CalcEProd, GeomFile,AnAGeom,AnaFile
344

    
345
  Namelist/cartlist/list
346
  Namelist/frozenlist/list
347
  Namelist/analist/nb
348

    
349

    
350
  Flag_Opt_Geom = .FALSE. ! Initialized.
351

    
352
!!!!!!!!!!!!!!!!!!!!!!!!!
353
  ! 
354
  ! We print the Version info in file
355
  !
356
  WRITE(*,'(1X,A)') "Opt'n Path v1.49 (c) PFL/PD 2007-2013"
357

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

    
494
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
495
     IOOUT=6
496
  CASE (0)
497
     IOIN=5
498
     IOOUT=6
499
  END SELECT
500

    
501

    
502
  ! We initiliaze variables
503
  Pi=dacos(-1.0d0)
504
  Nat=1
505
  Fact=1.1
506
  IGeomRef=-1
507
  NGeomI=1
508
  NGeomF=8
509
  IReparam=5
510
  IReparamT=-1
511
  ISpline=5
512
  debug=.False.
513
  MW=.TRUE.
514
  bohr=.False.
515
  Coord='ZMAT'
516
  Step_method='RFO'
517
  DebugFile='Path.valid'
518
  PathName="Path"
519
  MaxCyc=50
520
  SMax=0.1D0
521
  SThresh=-1.
522
  GThresh=-1.
523
  FTan=0.
524
  Hinv=.TRUE.
525
  IniHUp=.FALSE.
526
  HupNeighbour=.FALSE.
527
  HesUpd="EMPTY"
528
  HUpdate="EMPTY"
529
  FFrozen=.FALSE.
530
  FCart=.FALSE.
531
  List=0
532
  renum=.TRUE.
533
  OptReac=.FALSE.
534
  OptProd=.FALSE.
535
  CalcEReac=.FALSE.
536
  CalcEProd=.FALSE.
537
  EReac=0.d0
538
  EProd=0.d0 
539
  OptGeom=-1
540
  GeomFile="EMPTY"
541
  AnaGeom=.FALSE.
542
  AnaFile="EMPTY"
543
  Nb=0
544
  NbCom=0
545
  PathOnly=.FALSE.
546
  Prog='EMPTY'
547
  ProgExe='EMPTY'
548
  Runmode='SERIAL'
549
  Input='EMPTY'
550
  Poscar="POSCAR"
551
  AutoCart=.FALSE.
552
  VMD=.TRUE.
553
  WriteVASP=.FALSE.
554
  DynMaxStep=.FALSE.
555
  CalcName="EMPTY"
556
  ISuffix="EMPTY"
557
  OSuffix="EMPTY"
558
  NGintMax=10
559
  Align=.TRUE.
560

    
561
  ! Number of point used to compute the length of the path,
562
  ! and to reparameterize the path
563
  NMaxPtPath=-1 
564
  FrozTol=1e-4
565

    
566
  ! Read the namelist
567
  read (IOIN,path)
568

    
569
  debug=valid("Main").or.valid("Path")
570

    
571
  FCalcHess=valid("CalcHess")
572
  PathName=AdjustL(PathName)
573
  Coord=AdjustL(Coord)
574
  CALL UpCase(coord)
575

    
576
  Step_method=AdjustL(Step_method)
577
  CALL UpCase(Step_method)
578

    
579
  Prog=AdjustL(Prog)
580
  CALL UpCase(prog)
581

    
582
  Input=AdjustL(Input)
583
  CALL UpCase(input)
584

    
585
  Runmode=AdjustL(Runmode)
586
  Call UpCase(Runmode)
587
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"
588

    
589
  if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN
590
     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"
591
     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"
592
     CalcEReac=.FALSE.
593
     CalcEProd=.FALSE.
594
  END IF
595

    
596
  If (IReparamT<=0) IReparamT=IReparam
597

    
598
! We take care of HesUpd flag
599
! this is needed because some users might still use it
600
  if (HesUpd.NE."EMPTY") THEN
601
! We are pityless:
602
     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"
603
     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"
604
     STOP
605
  END IF
606

    
607
  IF (HUpdate=='EMPTY') THEN
608
     IF (OptGeom>=1) THEN
609
        HupDate="BFGS"
610
     ELSE
611
        HUpdate="MS"
612
     END IF
613
  END IF
614

    
615
  If ((COORD.EQ.'XYZ').OR.(COORD(1:4).EQ.'CART')) THEN 
616
     COORD='CART'
617
  END IF
618

    
619
  IF (COORD.EQ.'CART') THEN
620
     Renum=.FALSE.
621
     IF (OptGeom>0) THEN
622
        IGeomRef=OptGeom
623
     ELSE
624
        IGeomRef=1
625
     END IF
626
  END IF
627

    
628
  
629
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
630
  if (Prog.EQ.'2D') Prog="TEST2D"
631
  if (Prog.EQ.'TEST3') Prog="LEPS"
632

    
633

    
634
  Select Case (Prog)
635
    CASE ("EMPTY","G09","GAUSSIAN")
636
     Prog='GAUSSIAN'
637
     If (ProgExe=="EMPTY") ProgExe="g09"
638
     if (CalcName=="EMPTY") CalcName="Path"
639
     if (ISuffix=="EMPTY")  ISuffix="com"
640
     if (OSuffix=="EMPTY")  OSuffix="log"
641
    CASE ("G03")
642
     Prog='GAUSSIAN'
643
     If (ProgExe=="EMPTY") ProgExe="g03"
644
     if (CalcName=="EMPTY") CalcName="Path"
645
     if (ISuffix=="EMPTY")  ISuffix="com"
646
     if (OSuffix=="EMPTY")  OSuffix="log"
647
    CASE ("EXT")
648
     If (ProgExe=="EMPTY") ProgExe="./Ext.exe"
649
     if (CalcName=="EMPTY") CalcName="Ext"
650
     if (ISuffix=="EMPTY")  ISuffix="in"
651
     if (OSuffix=="EMPTY")  OSuffix="out"
652
    CASE ('MOPAC')
653
     If (ProgExe=="EMPTY") ProgExe="mopac"
654
     if (CalcName=="EMPTY") CalcName="Path"
655
     if (ISuffix=="EMPTY")  ISuffix="mop"
656
     if (OSuffix=="EMPTY")  OSuffix="out"
657
    CASE ("SIESTA")
658
     If (ProgExe=="EMPTY") ProgExe="siesta"
659
     if (CalcName=="EMPTY") CalcName="Path"
660
     if (ISuffix=="EMPTY")  ISuffix="fdf"
661
     if (OSuffix=="EMPTY")  OSuffix="out"
662
    CASE ("VASP")
663
     If (ProgExe=="EMPTY") ProgExe="VASP"
664
    CASE ("TM","TURBOMOLE")
665
     Prog="TURBOMOLE"
666
     If (ProgExe=="EMPTY") ProgExe="jobex -c 1 -keep"
667
    CASE ("TEST","CHAMFRE","TEST2D","LEPS")
668
     ProgExe=""
669
    CASE DEFAULT
670
     WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
671
     STOP
672
  END SELECT 
673

    
674
  if (Input.EQ.'EMPTY') THEN
675
     Select Case (Prog)
676
       CASE ("VASP")
677
          Input=Prog
678
       CASE ("CHAMFRE")
679
          Input="VASP"
680
       CASE DEFAULT
681
          Input='XYZ'
682
     END SELECT
683
    WRITE(*,*) "Input has been set to the default: ",INPUT
684
 END IF
685

    
686
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
687

    
688
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
689
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
690
     STOP
691
  END IF
692

    
693
  if (OptGeom.GE.1) THEN
694
     FPrintGeom=.FALSE.
695
     If (GeomFile/='EMPTY') THEN
696
        FPrintGeom=.TRUE.
697
     END IF
698
     IF (IGeomRef.LE.0) THEN
699
        IGeomRef=OptGeom
700
     ELSE
701
        IF (IGeomRef.NE.OptGeom) THEN
702
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
703
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
704

    
705
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
706
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
707
           STOP
708
        END IF
709
     END IF
710
  END IF
711

    
712
  ALLOCATE(Cart(Nat))
713
  Cart=0
714

    
715
  ALLOCATE(Frozen(Nat))
716
  Frozen=0
717

    
718
  IF (FCart) THEN
719
! We rewind the file to be sure to browse all of it to read the Cart namelist
720
     REWIND(IOIN)
721
     List=0
722
     READ(IOIN,cartlist)
723
     IF (COORD.NE.'CART') THEN
724
        Cart=List(1:Nat)
725
        if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
726
     ELSE
727
        WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
728
        WRITE(*,*) "I will discard the cartlist"
729
        Cart=0
730
     END IF
731
  END IF
732

    
733
  if (FFrozen) THEN
734

    
735
     REWIND(IOIN)
736
     List=0
737
     READ(IOIN,Frozenlist)
738
     Frozen=List(1:Nat)
739
     if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
740
  END IF
741

    
742
  If (FCart) Call Expand(Cart,NCart,Nat)
743
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
744

    
745
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
746
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
747

    
748
  if (AnaGeom) THEN
749
     REWIND(IOIN)
750
     READ(IOIN,AnaList,IOSTAT=IoS)
751
     IF (IOS==0) THEN ! we have a AnaList namelist 
752
        Call ReadAnaList
753
     END IF
754
     IF (AnaFile=="EMPTY") THEN
755
        AnaFile=TRIM(PathName) // '.dat'
756
     END IF
757
     OPEN(IODat,File=AnaFile)
758
  END IF
759

    
760
        
761

    
762
  if (NMaxPtPath.EQ.-1) then
763
     NMaxPtPath=min(50*NGeomF,2000)
764
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
765
  end if
766

    
767
  !  NMaxPtPath=10000
768

    
769
  ! We ensure that FTan is in [0.:1.]
770
  FTan=min(abs(FTan),1.d0)
771

    
772
! PFL 2011 Mar 14 ->
773
! Added some printing for analyses with Anapath
774
  if (debug) THEN
775
     WRITE(IOOUT,path)
776
  ELSE
777
! Even when debug is off we print NAT, NGEOMI, NGEOMF, MAXCYC
778
! and PATHNAME 
779
     WRITE(IOOUT,'(1X,A,I5)') "NAT = ", nat
780
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMI = ", NGeomI
781
     WRITE(IOOUT,'(1X,A,I5)') "NGEOMF = ", NGeomF
782
     WRITE(IOOUT,'(1X,A,I5)') "MAXCYC = ", MaxCYC
783
     WRITE(IOOUT,'(1X,A,1X,A)') "PATHNAME = ", Trim(PATHNAME)
784
  END IF
785

    
786
  FTan=1.-FTan
787

    
788
  !Thresholds...
789
  if (SThresh.LE.0) SThresh=0.01d0
790
  If (GThresh.LE.0) GThresh=SThresh/4.
791
  SThreshRms=Sthresh*2./3.
792
  GThreshRms=Gthresh*2./3.
793

    
794
  ! First batch of allocations. Arrays that depend only on NGeomI, Nat
795
  ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
796
  ALLOCATE(IndZmat(Nat,5),Order(Nat),Atome(Nat))
797
  ALLOCATE(MassAt(NAt),OrderInv(Nat))
798

    
799
  AtName=""
800
  MassAt=1.
801

    
802
  ! We read the initial geometries
803
  Call Read_Geom(input)
804

    
805
  ! We convert atome names into atom mass number
806
  DO I=1,NAt
807
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
808
     Atome(I)=ConvertNumAt(AtName(I))
809
  END DO
810

    
811
 If (AnaGeom) THEN
812
! We print the description of the varialbes in AnaFile
813
    Call PrintAnaList(IoDat)
814
 END IF
815

    
816
  ! PFL 23 Sep 2008 ->
817
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
818
  ! between the first geometry and the others.
819
  IF (NFroz.GT.0) THEN
820
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
821
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
822
     ListAtFroz=.FALSE.
823

    
824
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
825
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
826
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
827

    
828
     SELECT CASE (NFroz)
829
        ! We have Frozen atoms
830
        ! PFL 28 Oct 2008 ->
831
        ! It might happen that the geometries are translated/rotated.
832
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
833
        ! three frozen atoms.
834

    
835
     CASE (3:)
836
        DO I=1,NFroz
837
           ListAtFroz(Frozen(I))=.TRUE.
838
        END DO
839
     CASE (2)
840
        IAt1=Frozen(1)
841
        IAt2=Frozen(2)
842
        ListAtFroz(Iat1)=.TRUE.
843
        ListAtFroz(Iat2)=.TRUE.
844
        x2(1)=x1(Iat2)-x1(Iat1)
845
        y2(1)=y1(Iat2)-y1(Iat1)
846
        z2(1)=z1(Iat2)-z1(Iat1)
847
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
848
        Dist=-1.
849
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
850
        DO I=1,Nat
851
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
852
           x2(2)=x1(Iat2)-x1(I)
853
           y2(2)=y1(Iat2)-y1(I)
854
           z2(2)=z1(Iat2)-z1(I)
855
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
856
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
857
           if (abs(ca)<0.9) THEN
858
              IF (Norm2>Dist) THEN
859
                 Iat3=I
860
                 Dist=Norm2
861
              END IF
862
           END IF
863
        END DO
864
        ListAtFroz(Iat3)=.TRUE.
865
        if (debug) WRITE(*,*) "DBG Path, NFroz=2, third frozen=",Iat3
866
     CASE (1)
867
        IAt1=Frozen(1)
868
        ListAtFroz(Iat1)=.TRUE.
869
        Dist=-1.
870
        ! We scan all atoms to find one that is further from At1
871
        DO I=1,Nat
872
           if (I.EQ.Iat1) Cycle
873
           x2(1)=x1(Iat1)-x1(I)
874
           y2(1)=y1(Iat1)-y1(I)
875
           z2(1)=z1(Iat1)-z1(I)
876
           Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
877
           IF (Norm1>Dist) THEN
878
              Iat2=I
879
              Dist=Norm1
880
           END IF
881
        END DO
882
        IAt2=Frozen(2)
883
        ListAtFroz(Iat2)=.TRUE.
884
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, second frozen=",Iat2
885
        x2(1)=x1(Iat2)-x1(Iat1)
886
        y2(1)=y1(Iat2)-y1(Iat1)
887
        z2(1)=z1(Iat2)-z1(Iat1)
888
        Norm1=sqrt(x2(1)**2+y2(1)**2+z2(1)**2)
889
        Dist=-1.
890
        ! We scan all atoms to find one that is not aligned with Iat1-Iat2
891
        DO I=1,Nat
892
           if ((I.EQ.Iat1).OR.(I.EQ.Iat2)) Cycle
893
           x2(2)=x1(Iat2)-x1(I)
894
           y2(2)=y1(Iat2)-y1(I)
895
           z2(2)=z1(Iat2)-z1(I)
896
           Norm2=sqrt(x2(2)**2+y2(2)**2+z2(2)**2)
897
           ca=(x2(1)*x2(2)+y2(1)*y2(2)+z2(1)*Z2(2))/(Norm1*Norm2)
898
           if (abs(ca)<0.9) THEN
899
              IF (Norm2>Dist) THEN
900
                 Iat3=I
901
                 Dist=Norm2
902
              END IF
903
           END IF
904
        END DO
905
        ListAtFroz(Iat3)=.TRUE.
906
        if (debug) WRITE(*,*) "DBG Path, NFroz=1, third frozen=",Iat3
907
     END SELECT
908

    
909
     DO I=2,NGeomI
910
        ! First, we align the Ith geometry on the first one, using only
911
        ! the positions of the "frozen" atoms. (see above: it might be that
912
        ! ListAtFroz contains non frozen atoms usd to align the geometries
913
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
914
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
915
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
916
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
917

    
918

    
919
        FChkFrozen=.FALSE.
920
        DO J=1,NFroz
921
           IAt=Frozen(J)
922
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
923
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
924
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
925
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
926
              FChkFrozen=.TRUE.
927
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
928
              WRITE(*,*) "from geometry 1 to geometry ",I
929
           END IF
930
        END DO
931
     END DO
932
     IF (FChkFrozen) THEN
933
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
934
        STOP
935
     END IF
936
  END IF
937

    
938

    
939
  N3at=3*Nat
940
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
941

    
942

    
943
  Call ReadInput
944

    
945
  IF (COORD=="MIXED") Call TestCart(AutoCart)
946

    
947
  SELECT CASE(NFroz)
948
  CASE (2)
949
     IntFroz=1
950
  CASE (3)
951
     IntFroz=3
952
  CASE (4:)
953
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
954
  END SELECT
955

    
956
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
957
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
958
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
959

    
960
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
961
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
962
     COORD="ZMAT"
963
  END IF
964

    
965
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
966

    
967
  SELECT CASE(COORD)
968
  CASE ('ZMAT')
969
     NCoord=Nfree
970
     ALLOCATE(dzdc(3,nat,3,nat))
971
  CASE ('MIXED')
972
     SELECT CASE (NCart+NFroz)
973
     CASE (1) 
974
        NCoord=N3At-3
975
     CASE (2)
976
        NCoord=N3At-1
977
     CASE (3:)
978
        NCoord=N3At
979
     END SELECT
980
     ALLOCATE(dzdc(3,nat,3,nat))
981
  CASE ('BAKER')
982
     Symmetry_elimination=0
983
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
984
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
985
     ALLOCATE(BTransInv(NCoord,3*Nat))
986
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
987
  CASE ('HYBRID')
988
     NCoord=N3at
989
  CASE ('CART')
990
     NCoord=N3at
991
  END SELECT
992

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

    
995
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
996
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
997
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
998
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
999
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
1000
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
1001
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
1002
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
1003

    
1004
  ZeroVec=0.d0
1005
  DO I =1, NGeomF
1006
     IntTangent(I,:)=0.d0
1007
  END DO
1008

    
1009
  if (debug) THEN
1010
     Print *, 'L886, IntTangent(1,:)='
1011
     Print *, IntTangent(1,:)
1012
  END IF
1013

    
1014
  if (.NOT.OptReac) Energies(1)=Ereac
1015
  if (.NOT.OptProd) Energies(NGeomF)=EProd
1016
  MaxStep=SMax
1017

    
1018
  ! We initialize the mass array
1019
  if (MW) THEN
1020
     WRITE(*,*) 'Working in MW coordinates'
1021
     DO I=1,Nat
1022
        massAt(I)=Mass(Atome(I))
1023
     END DO
1024
  ELSE
1025
     DO I=1,Nat
1026
        MassAt(I)=1.d0
1027
     END DO
1028
  END IF
1029

    
1030
  WRITE(*,*) "Prog=",TRIM(PROG)
1031

    
1032
  ALLOCATE(FrozAtoms(Nat))
1033

    
1034
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
1035
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
1036
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
1037
  IF (NFroz.EQ.0) THEN
1038
     FrozAtoms=.TRUE.
1039
  ELSE
1040
     I=1
1041
     NFr=0
1042
     FrozAtoms=.False.
1043
     DO WHILE (Frozen(I).NE.0)
1044
        IF (Frozen(I).LT.0) THEN
1045
           DO J=Frozen(I-1),abs(Frozen(I))
1046
              IF (.NOT.FrozAtoms(J)) THEN
1047
                 NFr=NFr+1
1048
                 FrozAtoms(J)=.TRUE.
1049
              END IF
1050
           END DO
1051
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1052
           FrozAtoms(Frozen(I))=.TRUE.
1053
           NFr=NFr+1
1054
        END IF
1055
        I=I+1
1056
     END DO
1057
     IF (NFr.NE.NFroz) THEN
1058
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1059
        STOP
1060
     END IF
1061
  END IF
1062

    
1063
!  if (debug) THEN
1064
     !Print *, 'L968, IntTangent(1,:)='
1065
     !Print *, IntTangent(1,:)
1066
!  END IF
1067

    
1068
  ! We have read everything we needed... time to work :-)
1069
  IOpt=0
1070
  FirstTimePathCreate = .TRUE. ! Initialized.
1071
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1072
  FirstTimePathCreate = .FALSE.
1073

    
1074
  if (debug) THEN
1075
     Print *, 'L980, IntTangent(1,:)='
1076
     Print *, IntTangent(1,:)
1077
  END IF
1078

    
1079
  ! PFL 30 Mar 2008 ->
1080
  ! The following is now allocated in Calc_Baker. This is more logical to me
1081
  !   IF (COORD .EQ. "BAKER") Then
1082
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1083
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1084
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1085
  !   END IF
1086
  ! <- PFL 30 mar 2008
1087

    
1088
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1089
  !
1090
  ! For Debugging purposes
1091
  !
1092
  ! OptGeom is the index of the geometry which is to be optimized.
1093
  IF (Optgeom.GT.0) THEN
1094
     Flag_Opt_Geom=.TRUE.
1095
     SELECT CASE(Coord)
1096
     CASE ('CART','HYBRID')
1097
!!! CAUTION : PFL 29.AUG.2008 ->
1098
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1099
        ! This should be made  consistent !!!!!!!!!!!!!!!
1100
        GeomTmp=Reshape(reshape(XyzGeomI(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1101
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1102
!!! <- CAUTION : PFL 29.AUG.2008
1103
     CASE ('ZMAT','MIXED')
1104
        !GeomTmp=IntCoordF(OptGeom,:)
1105
        GeomTmp=IntCoordI(OptGeom,:)
1106
     CASE ('BAKER')
1107
        !GeomTmp=IntCoordF(OptGeom,:)
1108
        GeomTmp=IntCoordI(OptGeom,:)
1109
     CASE DEFAULT
1110
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1111
        STOP
1112
     END SELECT
1113

    
1114
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1115
     Flag_Opt_Geom = .TRUE.
1116
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1117

    
1118
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1119
     WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1120
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1121
     STOP
1122
  END IF ! IF (Optgeom.GT.0) THEN
1123

    
1124
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1125
  ! End of GeomOpt
1126
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1127

    
1128
  IF (PathOnly) THEN
1129
     WRITE(*,*) "PathOnly=.T. , Stop here"
1130
     Call Write_path(-1)
1131
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1132
     STOP
1133
  END IF
1134

    
1135
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1136

    
1137
  DEALLOCATE(XyzGeomI,IntCoordI)
1138
  NGeomI=NGeomF
1139
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1140

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

    
1146
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1147
     DO IGeom=1,NGeomF
1148
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1149
        DO I=1,3*Nat-6
1150
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1151
        END DO
1152
     END DO
1153
  END IF
1154

    
1155
  ! To calculate B^T^-1 for all extrapolated final geometries:
1156
  ! PFL 18 Jan 2008 ->
1157
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1158
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1159
  ! <-- PFL 18 Jan 2008
1160

    
1161
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1162
     DO IGeom=1,NGeomF
1163
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
1164
        DO I=1,3*Nat-6
1165
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1166
        END DO
1167
     END DO
1168
  END IF
1169

    
1170
  IOpt=0
1171
  Call EgradPath(IOpt)
1172

    
1173
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1174
  Call Write_path(IOpt)
1175
! Shall we analyze the geometries ?
1176
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1177

    
1178
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1179

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

    
1184
  ! By default, Hess=Id
1185
  Hess=0.
1186
  IF(Coord .EQ. "ZMAT") Then
1187
     ! We use the fact that we know that approximate good hessian values
1188
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1189
     DO IGeom=1,NGeomF
1190
        Hess(IGeom,1,1)=0.5d0
1191
        Hess(IGeom,2,2)=0.5d0
1192
        Hess(IGeom,3,3)=0.1d0
1193
        DO J=1,Nat-3
1194
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1195
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1196
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1197
        END DO
1198
        IF (HInv) THEN
1199
           DO I=1,NCoord
1200
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1201
           END DO
1202
        END IF
1203
     END DO
1204
  ELSE
1205
     DO I=1,NCoord
1206
        DO J=1,NGeomF
1207
           Hess(J,I,I)=1.d0
1208
        END DO
1209
     END DO
1210
  END IF
1211

    
1212
  IF (COORD .EQ. "BAKER") THEN
1213
     ! UMat(NPrim,NCoord)
1214
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1215
     Hprim=0.d0
1216
     ScanCoord=>Coordinate
1217
     I=0
1218
     DO WHILE (Associated(ScanCoord%next))
1219
        I=I+1
1220
        SELECT CASE (ScanCoord%Type)
1221
        CASE ('BOND')
1222
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1223
           Hprim(I,I) = 0.5d0
1224
           !END DO
1225
        CASE ('ANGLE')
1226
           Hprim(I,I) = 0.2d0
1227
        CASE ('DIHEDRAL')
1228
           Hprim(I,I) = 0.1d0
1229
        END SELECT
1230
        ScanCoord => ScanCoord%next
1231
     END DO
1232

    
1233
     ! HprimU = H^prim U:
1234
     HprimU=0.d0
1235
     DO I=1, 3*Nat-6
1236
        DO J=1,NPrim
1237
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1238
        END DO
1239
     END DO
1240

    
1241
     ! Hess = U^T HprimU = U^T H^prim U:
1242
     Hess=0.d0
1243
     DO I=1, 3*Nat-6
1244
        DO J=1,NPrim
1245
           ! UMat^T is needed below.
1246
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1247
        END DO
1248
     END DO
1249
     DO K=2,NGeomF
1250
        Hess(K,:,:)=Hess(1,:,:)
1251
     END DO
1252
     DEALLOCATE(Hprim,HprimU)
1253
  end if !  ends IF COORD=BAKER
1254
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1255
  ALLOCATE(HessTmp(NCoord,NCoord))
1256

    
1257
  IF (IniHUp) THEN
1258
     IF (FCalcHess) THEN
1259
        ! The numerical calculation of the Hessian has been switched on
1260
        DO IGeom=1,NGeomF
1261
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1262
              GeomTmp=IntCoordF(IGeom,:)
1263
           ELSE
1264
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1265
           END IF
1266
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1267
        END DO
1268
     ELSE
1269
        ! We generate a forward hessian and a backward one and we mix them.
1270
        ! First, the forward way:
1271
        DO IGeom=2,NGeomF-1
1272
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1273
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1274
           ELSE
1275
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1276
           END IF
1277
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1278
           HessTmp=Hess(IGeom-1,:,:)
1279
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1280
           Hess(IGeom,:,:)=HessTmp
1281
        END DO
1282

    
1283
        ! Now, the backward way:
1284
        HessTmp=Hess(NGeomF,:,:)
1285
        DO IGeom=NGeomF-1,2,-1
1286
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1287
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1288
           ELSE
1289
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1290
           END IF
1291
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1292
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1293
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1294
           ca=cos(alpha)
1295
           sa=sin(alpha)
1296
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1297
        END DO
1298
     END IF ! matches IF (FCalcHess)
1299
  END IF ! matches IF (IniHUp) THEN
1300

    
1301
  ! For debug purposes, we diagonalize the Hessian matrices
1302
  IF (debug) THEN
1303
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1304
     DO I=1,NGeomF
1305
        WRITE(*,*) "DBG Main - Point ",I
1306
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1307
        DO J=1,NCoord
1308
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1309
        END DO
1310
     END DO
1311
  END IF ! matches IF (debug) THEN
1312

    
1313
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1314
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1315
  !
1316
  ! Main LOOP : Optimization of the Path
1317
  !
1318
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1319
  if (debug)   Print *, 'NGeomF=', NGeomF
1320
  allocation_flag = .TRUE.
1321

    
1322
  Fini=.FALSE.
1323

    
1324
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1325
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1326
     IOpt=IOpt+1
1327

    
1328
     SELECT CASE (COORD)
1329
     CASE ('ZMAT','MIXED','BAKER')
1330
        GeomOld_all=IntCoordF
1331
     CASE ('CART','HYBRID')
1332
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1333
     CASE DEFAULT
1334
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1335
        STOP
1336
     END SELECT
1337

    
1338
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1339
          valid("printtangent")) THEN
1340
        WRITE(*,*) "Visualization of tangents"
1341
        Idx=min(12,NCoord)
1342
        DO I=1,NGeomF
1343
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1344
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1345
           WRITE(*,*) 
1346
        END DO
1347
        WRITE(*,*) "END of tangents"
1348
     END IF
1349

    
1350
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1351
          valid("printgrad")) THEN
1352
        WRITE(*,*) "Visualization of gradients"
1353
        Idx=min(12,NCoord)
1354
        DO I=1,NGeomF
1355
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1356
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1357
           WRITe(*,*) 
1358
        END DO
1359
        WRITE(*,*) "END of gradients"
1360
     END IF
1361

    
1362
     Fini=.TRUE.
1363
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1364
        IGeom=1
1365
        if (debug) WRITE(*,*) '**************************************'
1366
        if (debug) WRITE(*,*) 'Optimizing reactant'
1367
        if (debug) WRITE(*,*) '**************************************'
1368
        SELECT CASE (COORD)
1369
        CASE ('ZMAT','MIXED','BAKER')
1370
           SELECT CASE (Step_method)
1371
           CASE ('RFO')
1372
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
1373
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1374
              Print *, GradTmp
1375
           CASE ('GDIIS')
1376
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1377
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1378
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1379
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1380
              Print *, GradTmp
1381
           CASE ('GEDIIS')
1382
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1383
              ! Energies are updated in EgradPath.f90
1384
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1385
                   NCoord,allocation_flag,ZeroVec)
1386
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1387
              !GradTmp = GeomTmp - IntCoordF(1,:)
1388
              !Print *, 'Old Geometry:'
1389
              !Print *, IntCoordF(1,:)
1390
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1391
              Print *, GradTmp
1392
              !Print *, 'GeomTmp:'
1393
              !Print *, GeomTmp
1394
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1395
           CASE DEFAULT
1396
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1397
              STOP
1398
           END SELECT
1399
        CASE ('HYBRID','CART')
1400
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1401
        CASE DEFAULT
1402
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1403
           STOP
1404
        END SELECT
1405

    
1406
        IF (debug) THEN
1407
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1408
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1409
        END IF
1410

    
1411
        ! we have the Newton+RFO step, we will check that the maximum displacement is not greater than Smax
1412
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1413
        FactStep=min(1.d0,MaxStep(Igeom)/NormStep)
1414
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1415
        OptReac=(NormStep.GT.SThresh)
1416
        IF (debug) WRITE(*,*) "NormStep, SThresh, OptReac=",NormStep, SThresh, OptReac
1417
        NormGrad=sqrt(DOT_Product(reshape(Grad(1,:),(/Nfree/)),reshape(Grad(1,:),(/NFree/))))
1418
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1419
        OptReac=(OptReac.OR.(NormGrad.GT.GThresh))
1420
        !Print *, 'Grad(1,:):'
1421
        !Print *, Grad(1,:)
1422
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptReac=",NormGrad, GThresh, OptReac
1423

    
1424
        GradTmp=GradTmp*FactStep
1425

    
1426
        ! we take care that frozen atoms do not move.
1427
        IF (NFroz.NE.0) THEN
1428
           SELECT CASE (COORD)
1429
           CASE ('ZMAT','MIXED')
1430
              GradTmp(1:IntFroz)=0.d0
1431
           CASE ('CART','HYBRID')
1432
              DO I=1,Nat
1433
                 IF (FrozAtoms(I)) THEN
1434
                    Iat=Order(I)
1435
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1436
                 END IF
1437
              END DO
1438
           CASE('BAKER')
1439
              GradTmp(1:IntFroz)=0.d0 !GradTmp is "step" in Step_RFO_all(..)
1440
           CASE DEFAULT
1441
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1323.STOP'
1442
              STOP
1443
           END SELECT
1444
        END IF ! matches IF (NFroz.NE.0) THEN
1445

    
1446
        IF (debug) THEN
1447
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1448
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1449
        END IF
1450

    
1451
        IF (debug) THEN
1452
           !WRITE(*,*) "DBG Main, L1249, IOpt=",IOpt, ', IGeom=1'
1453
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1454
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1455
           ELSE
1456
              DO Iat=1,Nat
1457
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1458
                      XyzGeomF(IGeom,1:3,Iat)
1459
              END DO
1460
           END IF
1461
        END IF ! matches IF (debug) THEN
1462

    
1463
        SELECT CASE (COORD)
1464
        CASE ('ZMAT','MIXED','BAKER')
1465
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1466
        CASE ('HYBRID','CART')
1467
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1468
        CASE DEFAULT
1469
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1470
           STOP
1471
        END SELECT
1472

    
1473
        IF (debug) THEN
1474
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1475
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1476
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1477
           ELSE
1478
              DO Iat=1,Nat
1479

    
1480
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1481
                 ! It was with indzmat. IndZmat cannot appear here...
1482
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1483
                      XyzGeomF(IGeom,1:3,Iat)
1484
              END DO
1485
           END IF
1486
        END IF ! matches IF (debug) THEN
1487

    
1488
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1489
     ELSE ! Optreac
1490
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1491
     END IF ! matches IF (OptReac) THEN
1492

    
1493
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1494
        IGeom=NGeomF
1495
        if (debug) WRITE(*,*) '******************'
1496
        if (debug) WRITE(*,*) 'Optimizing product'
1497
        if (debug) WRITE(*,*) '******************'
1498
        SELECT CASE (COORD)
1499
        CASE ('ZMAT','MIXED','BAKER')
1500
           Print *, 'Optimizing product'
1501
           SELECT CASE (Step_method)
1502
           CASE ('RFO')
1503
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1504
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1505
              Print *, GradTmp
1506
           CASE ('GDIIS')
1507
              ! GradTmp becomes "step" below and has INTENT(OUT):
1508
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1509
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1510
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1511
              Print *, GradTmp
1512
           CASE ('GEDIIS')
1513
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1514
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1515
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1516
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1517
              Print *, GradTmp
1518
           CASE DEFAULT
1519
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1520
              STOP
1521
           END SELECT
1522
        CASE ('HYBRID','CART')
1523
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1524
        CASE DEFAULT
1525
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1526
           STOP
1527
        END SELECT
1528

    
1529
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1530
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1531
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1532
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1533
        OptProd=.NOT.(NormStep.LE.SThresh)
1534
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1535
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1536
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1537
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1538

    
1539
        GradTmp=GradTmp*FactStep
1540

    
1541
        ! we take care that frozen atoms do not move
1542
        IF (NFroz.NE.0) THEN
1543
           SELECT CASE (COORD)
1544
           CASE ('ZMAT','MIXED','BAKER')
1545
              GradTmp(1:IntFroz)=0.d0
1546
           CASE ('CART','HYBRID')
1547
              DO I=1,Nat
1548
                 IF (FrozAtoms(I)) THEN
1549
                    Iat=Order(I)
1550
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1551
                 END IF
1552
              END DO
1553
           CASE DEFAULT
1554
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1555
              STOP
1556
           END SELECT
1557
        END IF ! matches IF (NFroz.NE.0) THEN
1558

    
1559
        IF (debug) THEN
1560
           WRITE(*,*) "DBG Main, L1354, IGeom=, IOpt=",IGeom,IOpt
1561
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1562
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1563
           ELSE
1564
              DO Iat=1,Nat
1565
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1566
                      XyzGeomF(IGeom,1:3,Iat)
1567
              END DO
1568
           END IF
1569
        END IF
1570

    
1571
        SELECT CASE (COORD)
1572
        CASE ('ZMAT','MIXED','BAKER')
1573
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1574
        CASE ('HYBRID','CART')
1575
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1576
        CASE DEFAULT
1577
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1578
           STOP
1579
        END SELECT
1580

    
1581
        IF (debug) THEN
1582
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1583
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1584
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1585
           ELSE
1586
              DO Iat=1,Nat
1587
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1588
                      XyzGeomF(IGeom,1:3,Iat)
1589
              END DO
1590
           END IF
1591
        END IF
1592
     ELSE ! Optprod
1593
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1594
     END IF !matches IF (OptProd) THEN 
1595

    
1596
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1597

    
1598
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1599
     !
1600
     !  Optimizing other geometries
1601
     !
1602
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1603

    
1604

    
1605

    
1606
     DO IGeom=2,NGeomF-1 ! matches at L1556
1607
        if (debug) WRITE(*,*) '****************************'
1608
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1609
        if (debug) WRITE(*,*) '****************************'
1610

    
1611
        SELECT CASE (COORD)
1612
        CASE ('ZMAT','BAKER','MIXED')
1613
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
1614
           ! PFL 6 Apr 2008 ->
1615
           ! Special case, if FTan=0. we do not consider Tangent at all.
1616
           ! To DO: add the full treatment of FTan
1617
           if (FTan==0.) GradTmp2=ZeroVec
1618
           ! <- PFL 6 Apr 2008
1619
           if (debug) THEN
1620
              Print *, 'L1500, IntTangent(',IGeom,',:)='
1621
              Print *, IntTangent(IGeom,:)
1622
           END IF
1623
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
1624
           SELECT CASE (Step_method)
1625
           CASE ('RFO')
1626
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1627
           CASE ('GDIIS')
1628
              !The following has no effect at IOpt==1
1629
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
1630
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&        
1631
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1632
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
1633
              Print *, GradTmp
1634
           CASE ('GEDIIS')
1635
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1636
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
1637
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
1638
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
1639
              Print *, GradTmp
1640
           CASE DEFAULT
1641
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
1642
              STOP
1643
           END SELECT
1644
        CASE ('HYBRID','CART')
1645
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
1646
           ! so we change it
1647
           DO I=1,Nat
1648
              DO J=1,3
1649
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
1650
              END DO
1651
           END DO
1652
           !           GradTmp2=XyzTangent(IGeom,:)
1653
           ! PFL 6 Apr 2008 ->
1654
           ! Special case, if FTan=0. we do not consider Tangent at all.
1655
           ! To DO: add the full treatment of FTan
1656
           if (FTan==0.) GradTmp2=ZeroVec
1657
           ! <- PFL 6 Apr 2008
1658

    
1659
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
1660
        CASE DEFAULT
1661
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
1662
           STOP
1663
        END SELECT
1664

    
1665
        ! we take care that frozen atoms do not move
1666
        IF (NFroz.NE.0) THEN
1667
           SELECT CASE (COORD)
1668
           CASE ('ZMAT','MIXED','BAKER')
1669
              IF (debug) THEN 
1670
                 WRITE(*,*) 'Step computed. Coord=',Coord
1671
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
1672
              END IF
1673
              GradTmp(1:IntFroz)=0.d0
1674
           CASE ('CART','HYBRID')
1675
              if (debug) THEN
1676
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
1677
                 DO I=1,Nat
1678
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
1679
                 END DO
1680
              END IF
1681
              DO I=1,Nat
1682
                 IF (FrozAtoms(I)) THEN
1683
                    if (debug) THEN
1684
                       write(*,*) 'Step Computed. Coord=',Coord
1685
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
1686
                    END IF
1687
                    Iat=Order(I)
1688
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1689
                 END IF
1690
              END DO
1691

    
1692
                 if (debug) THEN
1693
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
1694
                    DO I=1,Nat
1695
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
1696
                    END DO
1697
                 END IF
1698

    
1699
           CASE DEFAULT
1700
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
1701
              STOP
1702
           END SELECT
1703
        END IF !matches IF (NFroz.NE.0) THEN
1704

    
1705
        IF (debug) THEN
1706
           SELECT CASE (COORD)
1707
           CASE ('ZMAT','MIXED','BAKER')
1708
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
1709
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1710
           CASE ('CART','HYBRID')
1711
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
1712
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
1713
           CASE DEFAULT
1714
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
1715
              STOP
1716
           END SELECT
1717
        END IF ! matches if (debug) THEN
1718

    
1719
        ! we check that the maximum displacement is not greater than Smax
1720
        If (debug) WRITE(*,*) "Fini before test:",Fini
1721
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1722
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
1723
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1724
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
1725

    
1726
        GradTmp=GraDTmp*FactStep
1727

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

    
1731
        IF (debug) THEN
1732
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
1733
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1734
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1735
           ELSE
1736
              DO Iat=1,Nat
1737
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1738
                      XyzGeomF(IGeom,1:3,Iat)
1739
              END DO
1740
           END IF
1741
        END IF ! matches if (debug) THEN
1742

    
1743
        SELECT CASE (COORD)
1744
        CASE ('ZMAT')
1745
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1746
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1747
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1748
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1749
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1750
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
1751
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1752
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
1753
           Idx=4
1754
           DO Iat=4,Nat
1755
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1756
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
1757
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
1758
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
1759
              Idx=Idx+3
1760
           END DO
1761

    
1762
        CASE ('BAKER')
1763
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1764
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1765
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1766
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
1767

    
1768
        CASE ('MIXED')
1769
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
1770
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1771
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1772
           DO Iat=1,NCart
1773
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
1774
           END DO
1775
           Idx=3*NCart+1
1776
           IF (NCart.EQ.1) THEN
1777
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1778
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
1779
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1780
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
1781

    
1782
              Idx=7
1783
           END IF
1784
           IF (NCart.EQ.2) THEN
1785
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1786
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
1787
              Idx=9
1788
           END IF
1789

    
1790
           
1791
           DO Iat=max(NCart+1,4),Nat
1792
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1793
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
1794
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
1795
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
1796
              Idx=Idx+3
1797
           END DO
1798
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
1799

    
1800
        CASE ('HYBRID','CART')
1801
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1802
        CASE DEFAULT
1803
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
1804
           STOP
1805
        END SELECT
1806

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

    
1819
        ! We project out the tangential part of the gradient to know if we are done
1820
        GradTmp=Grad(IGeom,:)
1821
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
1822
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
1823
        SELECT CASE (COORD)
1824
        CASE ('ZMAT','MIXED','BAKER')
1825
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
1826
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
1827
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
1828
           Ca=S/(sqrt(Norm)*NormGrad)
1829
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
1830
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
1831
        CASE ('CART','HYBRID')
1832
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
1833
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
1834
           Ca=S/(sqrt(Norm)*NormGrad)
1835
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
1836
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
1837
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
1838
        CASE DEFAULT
1839
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
1840
           STOP
1841
        END SELECT
1842
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
1843
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
1844
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1845
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
1846

    
1847
     END DO ! matches DO IGeom=2,NGeomF-1
1848
     ! we save the old gradients
1849
     GradOld=Grad
1850
     EPathOld=Energies
1851

    
1852

    
1853

    
1854
     ! Shall we re-parameterize the path ?
1855
     ! PFL 2007/Apr/10 ->
1856
     ! We call PathCreate in all cases, it will take care of the 
1857
     ! reparameterization, as well as calculating the tangents.
1858
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
1859
     XyzGeomI=XyzGeomF
1860
     IntCoordI=IntCoordF
1861

    
1862
     Call PathCreate(IOpt)
1863
     !     END IF
1864
     ! <- PFL 2007/Apr/10
1865

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

    
1869
        ! We have the new path, we calculate its energy and gradient
1870

    
1871
        Call EgradPath(IOpt)
1872
        !IF(IOPT .EQ. 10) Then
1873
        !         Print *, 'Stopping at my checkpoint.'
1874
        !           STOP !This is my temporary checkpoint.
1875
        !ENDIF
1876

    
1877
        ! PFL 31 Mar 2008 ->
1878
        ! Taken from v3.99 but we added a flag...
1879
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
1880
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
1881
        ! else it is decreased by 0.8
1882

    
1883
        if (DynMaxStep) THEN
1884
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
1885
           if (OptReac) THEN
1886
              If (Energies(1)<EPathOld(1)) Then
1887
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
1888
              ELSE
1889
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
1890
              END IF
1891
           END IF
1892

    
1893
           if (OptProd) THEN
1894
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
1895
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
1896
              ELSE
1897
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
1898
              END IF
1899
           END IF
1900

    
1901
           DO IGeom=2,NGeomF-1
1902
              If (Energies(IGeom)<EPathOld(IGeom)) Then
1903
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
1904
              ELSE
1905
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
1906
              END IF
1907
           END DO
1908
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
1909
        END IF ! test on DynMaxStep
1910
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
1911
        ! <- PFL 31 MAr 2008
1912
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
1913
        !  It should get updated for other cases too.
1914

    
1915
        DO IGeom=1,NGeomF
1916
           SELECT CASE (COORD)
1917
           CASE('ZMAT')
1918
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1919
              GeomTmp=IntCoordF(IGeom,:)
1920
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
1921
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
1922
                   OrderInv(indzmat(2,2)),Geomtmp(1)
1923
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
1924
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
1925
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
1926
              Idx=4
1927
              DO Iat=4,Nat
1928
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
1929
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
1930
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
1931
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
1932
                 Idx=Idx+3
1933
              END DO
1934
           CASE('BAKER')
1935
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1936
              !GeomTmp=IntCoordF(IGeom,:)
1937
           CASE('CART','HYBRID','MIXED')
1938
              WRITE(IOOUT,'(1X,I5)') Nat
1939
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
1940
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1941
              DO I=1,Nat
1942
                 Iat=I
1943
                 If (renum) Iat=OrderInv(I)
1944
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
1945
              END DO
1946
           CASE DEFAULT
1947
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
1948
              STOP
1949
           END SELECT
1950
        END DO ! matches DO IGeom=1,NGeomF
1951

    
1952
        Call Write_path(IOpt)
1953
! Shall we analyze the geometries ?
1954
     IF (AnaGeom) Call AnaPath(Iopt,IoDat)
1955

    
1956
        ! We update the Hessian matrices
1957
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
1958
        ! First using the displacement from the old path
1959
        IGeom0=2
1960
        IGeomF=NGeomF-1
1961
        IF (OptReac) IGeom0=1
1962
        IF (OptProd) IGeomF=NGeomF
1963

    
1964
        IF (FCalcHess) THEN
1965
           DO IGeom=IGeom0,IGeomF
1966
              SELECT CASE (COORD)
1967
              CASE ('ZMAT','MIXED','BAKER')
1968
                 GeomTmp2=IntCoordF(IGeom,:)
1969
              CASE ('CART','HYBRID')
1970
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1971
              CASE DEFAULT
1972
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
1973
                 STOP
1974
              END SELECT
1975
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
1976
           END DO
1977
        ELSE
1978
           DO IGeom=IGeom0,IGeomF
1979
              GeomTmp=GeomOld_all(IGeom,:)
1980
              SELECT CASE (COORD)
1981
              CASE ('ZMAT','MIXED','BAKER')
1982
                 GeomTmp2=IntCoordF(IGeom,:)
1983
              CASE ('CART','HYBRID')
1984
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1985
              CASE DEFAULT
1986
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
1987
                 STOP
1988
              END SELECT
1989
              GeomTmp=GeomTmp2-GeomTmp
1990
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
1991
              HessTmp=Hess(IGeom,:,:)
1992
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1993
              Hess(IGeom,:,:)=HessTmp
1994
           END DO ! matches DO IGeom=IGeom0,IGeomF
1995

    
1996
           ! Second using the neighbour points
1997
           IF (HupNeighbour) THEN
1998
              ! Only one point for end points.
1999
              IF (OptReac) THEN
2000
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2001
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
2002
                 ELSE
2003
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
2004
                 END IF
2005
                 GradTmp=Grad(1,:)-Grad(2,:)
2006
                 HessTmp=Hess(1,:,:)
2007
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2008
                 Hess(1,:,:)=HessTmp
2009
              END IF
2010

    
2011
              IF (OptProd) THEN
2012
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2013
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
2014
                 ELSE
2015
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
2016
                 END IF
2017
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
2018
                 HessTmp=Hess(NGeomF,:,:)
2019
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2020
                 Hess(NGeomF,:,:)=HessTmp
2021
              END IF
2022

    
2023
              ! Two points for the rest of the path
2024
              DO IGeom=2,NGeomF-1
2025
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
2026
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
2027
                 ELSE
2028
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
2029
                 END IF
2030
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
2031
                 HessTmp=Hess(IGeom,:,:)
2032
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2033

    
2034
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2035
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
2036
                 ELSE
2037
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
2038
                 END IF
2039
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
2040

    
2041
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2042
                 Hess(IGeom,:,:)=HessTmp
2043
              END DO ! matches DO IGeom=2,NGeomF-1
2044
           END IF !matches IF (HupNeighbour) THEN
2045
        END IF ! matches IF (FCalcHess)
2046
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
2047
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
2048

    
2049
!   IF (PROG=="GAUSSIAN") THEN
2050
!      DEALLOCATE(Gauss_paste)
2051
!      DEALLOCATE(Gauss_root)
2052
!      DEALLOCATE(Gauss_comment)
2053
!      DEALLOCATE(Gauss_end)
2054
!   END IF
2055

    
2056
  DEALLOCATE(XyzGeomF, IntCoordF)
2057
  DEALLOCATE(XyzGeomI, IntCoordI)
2058
  DEALLOCATE(XyzTangent,IntTangent)
2059
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2060
  DEALLOCATE(GeomOld_all)
2061

    
2062
  if (PROG=="GAUSSIAN") THEN
2063
     ! We de-allocate the Gauss_XX lists
2064
     ! transverse the list and deallocate each node
2065
     previous => Gauss_root ! make current point to head of list
2066
     DO
2067
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2068
        current => previous%next ! make list point to next node of head
2069
        DEALLOCATE(previous) ! deallocate current head node
2070
        previous => current  ! make current point to new head
2071
     END DO
2072

    
2073
     previous => Gauss_comment ! make current point to head of list
2074
     DO
2075
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2076
        current => previous%next ! make list point to next node of head
2077
        DEALLOCATE(previous) ! deallocate current head node
2078
        previous => current  ! make current point to new head
2079
     END DO
2080

    
2081
     previous => Gauss_end ! make current point to head of list
2082
     DO
2083
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2084
        current => previous%next ! make list point to next node of head
2085
        DEALLOCATE(previous) ! deallocate current head node
2086
        previous => current  ! make current point to new head
2087
     END DO
2088

    
2089

    
2090
  END IF
2091

    
2092
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2093
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2094
  IF (COORD.EQ.'BAKER') THEN
2095
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2096
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2097
  END IF
2098

    
2099
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2100

    
2101
  Close(IOIN)
2102
  Close(IOOUT)
2103
  IF (AnaGeom) Close(IODAT)
2104

    
2105
END PROGRAM PathOpt