Statistiques
| Révision :

root / src / Path.f90 @ 8

Historique | Voir | Annoter | Télécharger (89,18 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 namelis.
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
PROGRAM PathOpt
209

    
210
  use VarTypes
211
  use Path_module
212
  use Io_module
213

    
214
  IMPLICIT NONE
215

    
216
  CHARACTER(132) :: FileIn, FileOut, Line,Line2
217
  CHARACTER(32) :: Input,poscar
218

    
219
  INTEGER(KINT) :: I,J,K,IGeom,iat, IGeom0, IGeomF
220
  INTEGER(KINT) :: IOpt
221
  INTEGER(KINT) :: Idx, LineL,LenLine
222
  INTEGER(KINT) :: N3at, NFree, Nfr
223
  ! Temporary integer counter
224
  INTEGER(KINT) :: NTmp
225

    
226
  INTEGER(KINT) :: List(2*MaxFroz)
227

    
228
  REAL(KREAL) :: E,FactStep,B(3) !E=Energy
229
  REAL(KREAL) :: Alpha,sa,ca,norm,s,NormStep,NormGrad
230
  REAL(KREAL) :: EReac,EProd,HP,HEAT ! HEAT=E
231
  REAL(KREAL), ALLOCATABLE :: GeomTmp(:),GradTmp(:),ZeroVec(:) !NCoord
232
  REAL(KREAL), ALLOCATABLE :: GradOld(:,:) ! (NGeomF,NCoord)
233
  REAL(KREAL), ALLOCATABLE :: GeomTmp2(:),GradTmp2(:) ! NCoord
234
  REAL(KREAL), ALLOCATABLE :: HessTmp(:,:)
235
  REAL(KREAL), ALLOCATABLE :: Hprim(:,:) !(NPrim,NPrim)
236
  REAL(KREAL), ALLOCATABLE :: HprimU(:,:) !(NPrim,3*Nat-6))
237
  LOGICAL, SAVE :: allocation_flag
238

    
239
  ! Flag to see if frozen atoms are frozen (see below)
240
  LOGICAL :: FChkFrozen
241

    
242
  ! Energies for all points along the path at the previous iteration
243
  REAL(KREAL), ALLOCATABLE ::  EPathold(:) ! NGeomF, where it is deallocated: Prakash
244
  ! Maximum step allowed for this geometry
245
  REAL(KREAL), ALLOCATABLE ::  MaxStep(:) ! NGeomF, where it is deallocated: Prakash 
246

    
247
! these are used to read temporary coordinates
248
  REAL(KREAL) :: xtmp,ytmp,ztmp
249

    
250
  LOGICAL :: FFrozen,FCart
251

    
252
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
253
  !
254
  ! To analyse the frozen atoms
255
  !
256
  REAL(KREAL), ALLOCATABLE :: x1(:),y1(:),z1(:) ! nat
257
  REAL(KREAL), ALLOCATABLE :: x2(:),y2(:),z2(:) ! nat
258
  REAL(KREAL) :: MRot(3,3),Norm1,Norm2,Dist
259
  LOGICAL, ALLOCATABLE :: ListAtFroz(:)
260
  INTEGER(KINT) :: Iat1, Iat2, Iat3
261

    
262
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
263
  !
264
  ! For VASP
265
  !
266
  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
267
  INTEGER(KINT) :: NbType, NbTypeUser
268

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

    
271
  LOGICAL :: TChk
272
  LOGICAL :: Debug, Fini,Flag_Opt_Geom
273

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

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

    
284

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

    
291
     SUBROUTINE Read_Geom(Input)
292
       CHARACTER(32), INTENT(IN) :: input
293
     END SUBROUTINE READ_GEOM
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
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)') "OpenPath v1.4 (c) PFL/PD 2007-2011"
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 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(*,*) "MaxCyc: maximum number of iterations for the path optimization"
410
        WRITE(*,*) "Smax: Maximum length of a step during path optimization"
411
        WRITE(*,*) "SThresh: Step Threshold to consider that the path is stationary"
412
        WRITE(*,*) "GThresh: Gradient Threshold to consider that the path is stationary, only orthogonal part is taken"
413
        WRITE(*,*) "FTan: We moving the path, this gives the proportion of the displacement tangent to the path"
414
        WRITE(*,*) "      that is kept. FTan=1. corresponds to the full displacement, "
415
        WRITE(*,*) "      whereas FTan=0. gives a displacement orthogonal to the path."
416
        WRITE(*,*) "IReparam: The path is not reparameterised at each iteration. This gives the frequency of reparameterization."
417
        WRITE(*,*) "ISpline: By default, a linear interpolation is used to generate the path."
418
        WRITE(*,*) "         This option indicates the first step where spline interpolation is used."
419
        WRITE(*,*) "Boxtol: Real between 0. and 1. When doing periodic calculations, "
420
        WRITE(*,*) "        it might happen that an atom moves out of the unit cell."
421
        WRITE(*,*) "        Path detects this by comparing the displacement to boxtol:"
422
        WRITE(*,*) "        if an atom moves by more than Boxtol, then it is moved back into the unit cell. Default value: 0.5"
423
        WRITE(*,*) ""
424
        WRITE(*,*) "Arrays:"
425
        WRITE(*,*) "Rcov: Array containing the covalent radii of the first 80 elements."
426
        WRITE(*,*) "      You can modify it using, rcov(6)=0.8."
427
        WRITE(*,*) "Mass: Array containing the atomic mass of the first 80 elements."
428
        WRITE(*,*) "AtTypes: Name of the different atoms used in a VASP calculations."
429
        WRITe(*,*) "If not given, Path will read the POTCAR file."
430
        WRITE(*,*) ""
431
        WRITE(*,*) "Flags:"
432
        WRITE(*,*) "MW:  Flag. True if one wants to work in Mass Weighted coordinates. Default=.TRUE."
433
        WRITE(*,*) "Renum: Flag. True if one wants to reoder the atoms in the initial order. default is .TRUE. most of the time."
434
        WRITE(*,*) "OptProd: True if one wants to optimize the geometry of the products."
435
        WRITE(*,*) "OptReac: True if one wants to optimize the geometry of the reactants."
436
        WRITE(*,*) "CalcEReac: if TRUE the reactants energy will be computed. Default is FALSE. Not Compatible with RunMode=Para"
437
        WRITE(*,*) "CalcEProd: if TRUE the products energy will be computed. Default is FALSE. Not compatible with RunMode=Para"
438
        WRITE(*,*) "PathOnly:TRUE if one wants to generate the initial path, and stops." 
439
        WRITE(*,*) "Hinv: if True, then Hessian inversed is used."
440
        WRITE(*,*) "IniHup: if True, then Hessian inverse is extrapolated using the initial path calculations."
441
        WRITE(*,*) "HupNeighbour: if True, then Hessian inverse is extrapolated using the neighbouring points of the path."
442
        WRITE(*,*) "FFrozen: True if one wants to freeze the positions of some atoms." 
443
        WRITE(*,*) "         If True, a &frozenlist namelist containing the list of frozen atoms must be given."
444
        WRITE(*,*) "          If VASP is used, and frozen is not given"
445
        WRITE(*,*) " here, the program will use the F flags of the input geometry"
446
        WRITE(*,*) "FCart:  True if one wants to describe some atoms using cartesian coordinates. "
447
        WRITE(*,*) "        *** Only used in 'mixed' calculations. ***"
448
        WRITE(*,*) "      If True, a &cartlist namelist containing the list of cart atoms must be given."
449
        WRITE(*,*) "      By default, only frozen atoms are described in cartesian coordinates."
450
        WRITE(*,*) ""
451
        WRITE(*,*) "DynMaxStep: if TRUE, the maximum allowed step is updated at each step, for each geometry."
452
        WRITE(*,*) "        If energy goes up, Smax=Smax*0.8, if not Smax=Smax*1.2. "
453
        WRITE(*,*) "       It is ensured that the dynamical Smax is within [0.5*SMax_0,2*Smax_0]"
454
        WRITE(*,*) "Autocart: True if you want to let the program choosing the cartesian atoms."
455
        WRITE(*,*) "VMD: TRUE if you want to use VMD to look at the Path. Used only for VASP for now"
456
        WRITE(*,*) "WriteVASP: TRUE if you want to print the images coordinates in POSCAR files."
457
        WRITE(*,*) "See also the POSCAR option. This can be used only if prog or input=VASP."
458
        WRITE(*,*) ""
459
        STOP
460
     END IF
461
     OPEN(IOIN,FILE=FileIn,STATUS='UNKNOWN')
462
     IOOUT=6
463
  CASE (0)
464
     IOIN=5
465
     IOOUT=6
466
  END SELECT
467

    
468

    
469

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

    
524
  ! Number of point used to compute the length of the path,
525
  ! and to reparameterize the path
526
  NMaxPtPath=-1 
527
  FrozTol=1e-4
528

    
529
  ! Read the namelist
530
  read (IOIN,path)
531

    
532
  debug=valid("Main")
533
  FCalcHess=valid("CalcHess")
534
  PathName=AdjustL(PathName)
535
  Coord=AdjustL(Coord)
536
  CALL UpCase(coord)
537

    
538
  Step_method=AdjustL(Step_method)
539
  CALL UpCase(Step_method)
540

    
541
  Prog=AdjustL(Prog)
542
  CALL UpCase(prog)
543

    
544
  Input=AdjustL(Input)
545
  CALL UpCase(input)
546

    
547
  Runmode=AdjustL(Runmode)
548
  Call UpCase(Runmode)
549
  IF (Runmode(1:4).EQ.'PARA')    Runmode="PARA"
550

    
551
  if ((RunMode.EQ."PARA").AND.(CalcEReac.OR.CalcEProd)) THEN
552
     WRITE(*,*) "CalcEReac and CalcEProd are not compatible (for now) with RunMode='PARA'"
553
     WRITE(*,*) "Setting CalcERead and CalcEProd to FALSE"
554
     CalcEReac=.FALSE.
555
     CalcEProd=.FALSE.
556
  END IF
557

    
558
  If (IReparamT<=0) IReparamT=IReparam
559

    
560
! We take care of HesUpd flag
561
! this is needed because some users might still use it
562
  if (HesUpd.NE."EMPTY") THEN
563
! We are pityless:
564
     WRITE(IOOUT,*) "HesUpd is deprecated. Use Hupdate instead"
565
     WRITE(*,*) "HesUpd is deprecated. Use Hupdate instead"
566
     STOP
567
  END IF
568

    
569

    
570
  IF (COORD.EQ.'CART') THEN
571
     Renum=.FALSE.
572
     IGeomRef=1
573
  END IF
574

    
575
  if (Prog.EQ.'EMPTY') THEN
576
     Prog='GAUSSIAN'
577
     If (ProgExe=="EMPTY") ProgExe="g09"
578
  END IF
579
  if (Prog.EQ.'G03') THEN
580
     Prog='GAUSSIAN'
581
     If (ProgExe=="EMPTY") ProgExe="g03"
582
  END IF
583
  if (Prog.EQ.'G09') THEN
584
     Prog='GAUSSIAN'
585
     If (ProgExe=="EMPTY") ProgExe="g09"
586
  END IF
587

    
588
  if (Prog.EQ.'TM') Prog="TURBOMOLE"
589
  if (Prog.EQ.'TEST2') Prog="CHAMFRE"
590

    
591

    
592
  Select Case (Prog)
593
    CASE ("EXT")
594
     if (CalcName=="EMPTY") CalcName="Ext"
595
     if (ISuffix=="EMPTY")  ISuffix="in"
596
     if (OSuffix=="EMPTY")  OSuffix="out"
597
    CASE ('MOPAC')
598
     if (CalcName=="EMPTY") CalcName="Path"
599
     if (ISuffix=="EMPTY")  ISuffix="mop"
600
     if (OSuffix=="EMPTY")  OSuffix="out"
601
    CASE ("GAUSSIAN")
602
     if (CalcName=="EMPTY") CalcName="Path"
603
     if (ISuffix=="EMPTY")  ISuffix="com"
604
     if (OSuffix=="EMPTY")  OSuffix="log"
605
    CASE DEFAULT
606
     if (CalcName=="EMPTY") CalcName="Path"
607
     if (ISuffix=="EMPTY")  ISuffix="in"
608
     if (OSuffix=="EMPTY")  OSuffix="out"
609
  END SELECT 
610

    
611
  if (Input.EQ.'EMPTY') THEN
612
     Select Case (Prog)
613
       CASE ("VASP")
614
          Input=Prog
615
       CASE ("CHAMFRE")
616
          Input="VASP"
617
       CASE DEFAULT
618
          Input='XYZ'
619
     END SELECT
620
    WRITE(*,*) "Input has been set to the default: ",INPUT
621
 END IF
622

    
623
  WriteVASP=WriteVASP.AND.((Prog=="VASP").OR.(Input=="VASP"))
624

    
625
  IF ((RunMode=="PARA").AND.((Prog/="GAUSSIAN").AND.(Prog/="VASP"))) THEN
626
     WRITE(*,*) "Program=",TRIM(Prog)," not yet supported for Para mode."
627
     STOP
628
  END IF
629

    
630
  if (OptGeom.GE.1) THEN
631
     HUpdate="BFGS"
632
     IF (IGeomRef.LE.0) THEN
633
        IGeomRef=OptGeom
634
     ELSE
635
        IF (IGeomRef.NE.OptGeom) THEN
636
           WRITE(IOOUT,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
637
           WRITE(IOOUT,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
638

    
639
           WRITE(*,*) "STOP: IGeomRef and OptGeom both set... but to different values !"
640
           WRITE(*,*) "Check it : IGeomRef=",IGeomRef," OptGeom=",OptGeom
641
           STOP
642
        END IF
643
     END IF
644
  END IF
645

    
646
  IF (ProgExe.EQ.'EMPTY') THEN
647
     SELECT CASE (PROG)
648
     CASE ("GAUSSIAN")
649
        ProgExe="g03"
650
     CASE ("MOPAC")
651
        ProgExe="mopac"
652
     CASE ("EXT")
653
        ProgExe="./Ext.exe"
654
     CASE ("VASP")
655
        ProgExe='VASP'
656
     CASE ("TURBOMOLE")
657
        ProgExe='jobex -c 1 -keep'
658
     CASE ("TEST")
659
        ProgExe=""
660
     CASE ("CHAMFRE")
661
        ProgExe=""
662
     CASE DEFAULT
663
        WRITE(*,*) "Prog=",Adjustl(trim(Prog))," not recognized. STOP"
664
        STOP
665
     END SELECT
666
  END IF
667

    
668
  if (FCart.AND.FFrozen) THEN
669
     ! We have to be carreful because user might use any order
670
     ! to give the two lists CartList and FrozenList
671
     READ(IOIN,'(A)') Line
672
     Call Upcase(Line)
673
     if (Index(Line,"CART").NE.0) THEN
674
        ALLOCATE(Cart(Nat))
675
        BackSpace(IOIN)
676
        List=0
677
        READ(IOIN,cartlist)
678
        IF (COORD.NE.'CART') THEN
679
           Cart=List(1:Nat)
680
           if (debug) WRITE(*,*) "DBG Path L512 Cart",Cart
681
        ELSE
682
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
683
           WRITE(*,*) "I will discard the cartlist"
684
           Cart=0
685
        END IF
686

    
687
        ALLOCATE(Frozen(Nat))
688
        BACKSPACE(IOIN)
689
        List=0
690
        READ(IOIN,Frozenlist)
691
        Frozen=List(1:Nat)
692
        if (debug) WRITE(*,*) "DBG Path L517 Frozen",Frozen
693
     ELSE
694
        ALLOCATE(Frozen(Nat))
695
        BACKSPACE(IOIN)
696
        List=0
697
        READ(IOIN,Frozenlist)
698
        Frozen=List(1:Nat)
699
        if (debug) WRITE(*,*) "DBG Path L523 Frozen",Frozen
700
        ALLOCATE(Cart(Nat))
701
        BackSpace(IOIN)
702
        List=0
703
        READ(IOIN,cartlist)
704
        IF (COORD.NE.'CART') THEN
705
           Cart=List(1:Nat)
706
           if (debug) WRITE(*,*) "DBG Path L548 Cart",Cart
707
        ELSE
708
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
709
           WRITE(*,*) "I will discard the cartlist"
710
           Cart=0
711
        END IF
712
     END IF
713
  ELSE
714
     if (FCart) THEN
715
        ALLOCATE(Cart(Nat))
716
        List=0
717
        READ(IOIN,cartlist)
718
        IF (COORD.NE.'CART') THEN
719
           Cart=List(1:Nat)
720
           if (debug) WRITE(*,*) "DBG Path L561 Cart",Cart
721
        ELSE
722
           WRITE(*,*) "FCart=T is not allowed when Coord='CART'"
723
           WRITE(*,*) "I will discard the cartlist"
724
           Cart=0
725
        END IF
726
     ELSE
727
        IF ((PROG=="VASP").OR.(AutoCart)) THEN
728
           ALLOCATE(Cart(Nat))
729
        ELSE
730
           ALLOCATE(Cart(1))
731
        END IF
732
        Cart=0
733
     END IF
734

    
735
     if (FFrozen) THEN
736
        ALLOCATE(Frozen(Nat))
737
        List=0
738
        READ(IOIN,Frozenlist)
739
        Frozen=List(1:Nat)
740
        WRITE(*,*) Frozen
741
        if (debug) WRITE(*,*) "DBG Path L549 Frozen",Frozen
742
     ELSE
743
        IF (PROG=="VASP") THEN
744
           ALLOCATE(Frozen(Nat))
745
        ELSE
746
           ALLOCATE(Frozen(1))
747
        END IF
748
        Frozen=0
749
     END IF
750
  END IF
751

    
752

    
753
  If (FCart) Call Expand(Cart,NCart,Nat)
754
  IF (FFrozen) Call Expand(Frozen,NFroz,Nat)
755

    
756
  if (debug) WRITE(*,*) "DBG Main L609, Ncart,Cart",NCart,Cart
757
  if (debug) WRITE(*,*) "DBG Main L610, NFroz,Frozen", NFroz,Frozen
758

    
759
  if (NMaxPtPath.EQ.-1) then
760
     NMaxPtPath=min(50*NGeomF,2000)
761
     NMaxPtPath=max(20*NGeomF,NMaxPtPath)
762
  end if
763

    
764
  !  NMaxPtPath=10000
765

    
766
  ! We ensure that FTan is in [0.:1.]
767
  FTan=min(abs(FTan),1.d0)
768

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

    
783
  FTan=1.-FTan
784

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

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

    
796
  ! We read the initial geometries
797
  Call Read_Geom(input)
798

    
799
  ! We convert atome names into atom mass number
800
  DO I=1,NAt
801
     !     WRITE(*,*) 'I,AtName',I,AtName(I)
802
     Atome(I)=ConvertNumAt(AtName(I))
803
  END DO
804

    
805

    
806
  ! PFL 23 Sep 2008 ->
807
  ! We check that frozen atoms are really frozen, ie that their coordinates do not change
808
  ! between the first geometry and the others.
809
  IF (NFroz.GT.0) THEN
810
     ALLOCATE(ListAtFroz(Nat),x1(Nat),y1(nat),z1(nat))
811
     ALLOCATE(x2(nat),y2(nat),z2(nat)) 
812
     ListAtFroz=.FALSE.
813

    
814
     x1(1:Nat)=XyzgeomI(1,1,1:Nat)
815
     y1(1:Nat)=XyzgeomI(1,2,1:Nat)
816
     z1(1:Nat)=XyzgeomI(1,3,1:Nat)
817

    
818
     SELECT CASE (NFroz)
819
        ! We have Frozen atoms
820
        ! PFL 28 Oct 2008 ->
821
        ! It might happen that the geometries are translated/rotated.
822
        ! So if we have 3 (or more) frozen atoms, we re-orient the geometries, based on the first
823
        ! three frozen atoms.
824

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

    
899
     DO I=2,NGeomI
900
        ! First, we align the Ith geometry on the first one, using only
901
        ! the positions of the "frozen" atoms. (see above: it might be that
902
        ! ListAtFroz contains non frozen atoms usd to align the geometries
903
        x2(1:Nat)=XyzgeomI(I,1,1:Nat)
904
        y2(1:Nat)=XyzgeomI(I,2,1:Nat)
905
        z2(1:Nat)=XyzgeomI(I,3,1:Nat)
906
        Call Alignpartial(Nat,x1,y1,z1,x2,y2,z2,ListAtFroz,MRot)
907

    
908

    
909
        FChkFrozen=.FALSE.
910
        DO J=1,NFroz
911
           IAt=Frozen(J)
912
           IF ((abs(X1(Iat)-X2(Iat)).GT.FrozTol).OR. &
913
                (abs(y1(Iat)-y2(Iat)).GT.FrozTol).OR.  &
914
                (abs(z1(Iat)-z2(Iat)).GT.FrozTol))                     THEN
915
              IF (.NOT.FChkFrozen) WRITE(*,*) "*** WARNING ****"
916
              FChkFrozen=.TRUE.
917
              WRITE(*,*) "Atom ",Iat," is set to Frozen but its coordinates change"
918
              WRITE(*,*) "from geometry 1 to geometry ",I
919
           END IF
920
        END DO
921
     END DO
922
     IF (FChkFrozen) THEN
923
        WRITE(*,*) "Some frozen atoms are not frozen... check this and play again"
924
        STOP
925
     END IF
926
  END IF
927

    
928

    
929
  N3at=3*Nat
930
  NFree=N3at-6   ! ATTENTION: THIS IS NOT VALID FOR ALL TYPES OF COORDINATES.
931

    
932
  ! if prog=gaussian, we have to read a full input
933
  ! to mimic it later
934
  if (prog=='GAUSSIAN') THEN
935
     ! We read the Gaussian input file
936
     ! First, the root
937
     IF (DEBUG) WRITE(*,*) "Reading Gauss Root"
938
     ALLOCATE(Gauss_Root)
939
     NULLIFY(Gauss_Root%next)
940
     Current => Gauss_root
941
     LineL=1
942
     DO WHILE (LineL.NE.0)
943
        READ(IOIN,'(A)') Line
944
        Line=AdjustL(Line)
945
        LineL=len(Trim(Line))
946
        Idx=INDEX(Line,"chk")
947
        IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN
948
           current%Line=TRIM(Line)
949
           ALLOCATE(current%next)
950
           Current => Current%next
951
           Nullify(Current%next)
952
        END IF
953
     END DO
954

    
955
!     Current => Gauss_root
956
!     DO WHILE (ASSOCIATED(Current%next))
957
!        WRITE(*,'(1X,A)') Trim(current%line)
958
!        Current => current%next
959
!     END DO
960

    
961
     ! Now the comment... 
962
     IF (DEBUG) WRITE(*,*) "Reading Gauss Comment"
963
     ALLOCATE(Gauss_Comment)
964
     NuLLIFY(Gauss_Comment%Next)
965
     Current => Gauss_comment
966
     LineL=1
967
     DO WHILE (LineL.NE.0)
968
        READ(IOIN,'(A)') Line
969
        Line=AdjustL(Line)
970
        LineL=len(Trim(Line))
971
        IF (LineL.NE.0) THEN
972
           current%Line=TRIM(Line)
973
           ALLOCATE(current%next)
974
           Current => Current%next
975
           Nullify(Current%next)
976
        END IF
977
     END DO
978

    
979
 !    Current => Gauss_comment
980
 !    DO WHILE (ASSOCIATED(Current%next))
981
 !       WRITE(*,'(1X,A)') Trim(current%line)
982
 !       Current => current%next
983
 !    END DO
984

    
985
     ! Now the charge
986
     IF (DEBUG) WRITE(*,*) "Reading Gauss Charge"
987
     READ(IOIN,'(A)') Gauss_Charge
988
!     WRITE(*,*) Gauss_charge
989
     ! We now read the Paste part...
990
     ALLOCATE(Gauss_Paste(NAt))
991
     LenLine=1
992
     Iat=0
993
     Gauss_paste=" "
994
     DO While (LenLine.GT.0)
995
        READ(IOIN,'(A)') Line
996
        Line=AdjustL(Line)
997
        LenLine=Len_TRIM(Line)
998
        IF (LenLine.GT.0) Iat=Iat+1
999
        Idx=Index(Line,'.',BACK=.TRUE.)
1000
        !          WRITE(*,*) Idx,TRIM(Line)
1001
        !          WRITE(*,*) Idx,TRIM(Line(Idx+1:))
1002
        Line=ADJUSTL(Line(Idx+1:))
1003
        Idx=Index(Line," ")
1004
        LineL=LEN_TRIM(Line(Idx:))
1005
        !          WRITE(*,*) LineL,TRIM(Line(Idx:))
1006
        IF (LineL.GT.0) THEN 
1007
           Gauss_paste(Iat)=Line(Idx:)
1008
        END IF
1009
     END DO
1010
     IF (Iat.NE.Nat) THEN
1011
        WRITE(*,*) "I found ", Iat," lines for the geometry instead of ",Nat
1012
        WRITE(*,*) "ERROR. STOP"
1013
        WRITE(IOOUT,*) "I found ", Iat," lines for the geometry instead of ",Nat
1014
        WRITE(IOOUT,*) "ERROR. STOP"
1015
        STOP
1016
     END IF
1017
     ! We now read the last part
1018
     IF (DEBUG) WRITE(*,*) "Reading Gauss End"
1019
     !     READ(IOIN,'(A)') Line
1020
     ALLOCATE(Gauss_End)
1021
     NuLLIFY(Gauss_End%Next)
1022
     Current => Gauss_End
1023
     LineL=1
1024
     DO WHILE (1.EQ.1)
1025
        READ(IOIN,'(A)',END=999) Line
1026
        Line=AdjustL(Line)
1027
        LineL=len(Trim(Line))
1028
        current%Line=TRIM(Line)
1029
        ALLOCATE(current%next)
1030
        Current => Current%next
1031
        Nullify(Current%next)
1032
     END DO
1033
999  CONTINUE
1034

    
1035
     !     ! Write the gaussian input file for testing purposes
1036
     !     Current => Gauss_root
1037
     !     DO WHILE (ASSOCIATED(Current%next))
1038
     !        WRITE(*,'(1X,A)') Trim(current%line)
1039
     !        Current => current%next
1040
     !     END DO
1041

    
1042
     !        WRITE(*,*) 
1043
     !        WRITE(*,*) 'Comment original'
1044

    
1045
     !        Current => Gauss_comment
1046
     !        DO WHILE (ASSOCIATED(Current%next))
1047
     !           WRITE(*,'(1X,A)') Trim(current%line)
1048
     !           Current => current%next
1049
     !        END DO
1050

    
1051
     !        WRITE(*,*) 
1052
     !        WRITE(*,*) Trim(Gauss_charge)
1053

    
1054
     !        DO I=1,Nat
1055
     !           WRITE(*,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(I)),XyzGeomI(1:3,I,1), TRIM(Gauss_Paste(I))
1056
     !        END DO
1057

    
1058
     !        WRITE(*,*) 
1059
     !        Current => Gauss_End
1060
     !        DO WHILE (ASSOCIATED(Current%next))
1061
     !           WRITE(*,'(1X,A)') Trim(current%line)
1062
     !           Current => current%next
1063
     !        END DO
1064

    
1065
     !        WRITE(*,*) 
1066

    
1067
  END IF
1068

    
1069
! if prog=mopac, we have to read a full input to reproduce it later
1070
! The structure is:
1071
! A MOPAC data set normally consists of one line of keywords, two lines of user-defined text, then the coordinates
1072
! Then  a blank line or a line of 0.
1073
! then the symmetry description.
1074
! comment lines start with * and can be anywhere !!!
1075

    
1076
  if (prog=='MOPAC') THEN
1077
     ! We read the Gaussian input file
1078
     ! First, the root
1079
     IF (DEBUG) WRITE(*,*) "Reading Mopac input"
1080
     ALLOCATE(Mopac_Root)
1081
     NULLIFY(Mopac_Root%next)
1082
     ALLOCATE(Mopac_Comment)
1083
     NULLIFY(Mopac_Comment%next)
1084
     ALLOCATE(Mopac_End)
1085
     NuLLIFY(Mopac_End%Next)
1086
     Current => Mopac_root
1087
     CurCom => Mopac_Comment
1088
     LineL=1
1089
     NTmp=0
1090
     DO WHILE (NTmp.LT.3)
1091
        READ(IOIN,'(A)') Line
1092
        Line=AdjustL(Line)
1093
        LineL=len(Trim(Line))
1094
        IF (Line(1:1)/="*") THEN
1095
           IF (NTmp==0) THEN
1096
              Line2=Line
1097
              Call UpCase(Line2)
1098
              Idx=Index(Line2,'GRADIENTS')
1099
              If (Idx==0) Line=TRIM(Line) // " GRADIENTS"
1100
              Idx=Index(Line2,'1SCF')
1101
              If (Idx==0) Line=TRIM(Line) // " 1SCF"
1102
           END IF
1103
           current%Line=TRIM(Line)
1104
           ALLOCATE(current%next)
1105
           Current => Current%next
1106
           Nullify(Current%next)
1107
           NTmp=NTmp+1
1108
        ELSE
1109
           CurCom%Line=TRIM(LINE)
1110
           ALLOCATE(CurCom%Next)
1111
           CurCom => CurCom%Next
1112
           NULLIFY(CurCom%Next)
1113
        END IF
1114
     END DO
1115

    
1116
!     Current => Mopac_root
1117
!     DO WHILE (ASSOCIATED(Current%next))
1118
!        WRITE(*,'(1X,A)') Trim(current%line)
1119
!        Current => current%next
1120
!     END DO
1121

    
1122
! Now the geometry... that we just skip :)
1123
     IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry"
1124
     Mopac_EndGeom=""
1125
     LineL=1
1126
     DO WHILE (LineL.NE.0)
1127
        READ(IOIN,'(A)',END=989) Line
1128
        Line=AdjustL(Line)
1129
        LineL=len(Trim(Line))
1130
! The last line might be either blank or filled with 0
1131
        IF (Line(1:1)=="0") THEN
1132
           LineL=0
1133
           Mopac_EndGeom=Trim(Line)
1134
        END IF
1135
        IF (Line(1:1)=="*") THEN
1136
           CurCom%Line=TRIM(LINE)
1137
           ALLOCATE(CurCom%Next)
1138
           CurCom => CurCom%Next
1139
           NULLIFY(CurCom%Next)
1140
        END IF
1141
     END DO
1142

    
1143
! If we are here, there might be something else to read: Mopac_end
1144

    
1145
     ! We now read the last part
1146
     IF (DEBUG) WRITE(*,*) "Reading Mopac End"
1147
     !     READ(IOIN,'(A)') Line
1148
     Current => Mopac_End
1149
     LineL=1
1150
     DO WHILE (1.EQ.1)
1151
        READ(IOIN,'(A)',END=989) Line
1152
        Line=AdjustL(Line)
1153
        LineL=len(Trim(Line))
1154
        IF (Line(1:1)/="*") THEN
1155
           current%Line=TRIM(Line)
1156
           ALLOCATE(current%next)
1157
           Current => Current%next
1158
           Nullify(Current%next)
1159
           NTmp=NTmp+1
1160
        ELSE
1161
           CurCom%Line=TRIM(LINE)
1162
           ALLOCATE(CurCom%Next)
1163
           CurCom => CurCom%Next
1164
           NULLIFY(CurCom%Next)
1165
        END IF
1166
     END DO
1167
989  CONTINUE
1168

    
1169
  END IF ! IF (PROG="MOPAC")
1170

    
1171
  if ((Prog=="VASP").AND.(input/="VASP")) THEN
1172
     ! Input was not Vasp, so many parameters are missing like lattice 
1173
     ! constants...
1174
     ! we read them now !
1175
     ALLOCATE(FFF(3,nat))
1176
     ! First geometry is a bit special for VASP as we have to set
1177
     ! many things
1178
     IF (DEBUG) WRITE(*,*) "Reading Vasp Parameters"
1179
     READ(IOIN,'(A)') Vasp_Title
1180
     READ(IOIN,*) Vasp_param
1181

    
1182
     READ(IOIN,*) Lat_a
1183
     READ(IOIN,*) Lat_b
1184
     READ(IOIN,*) Lat_c
1185

    
1186
     Lat_a=Lat_a*Vasp_param
1187
     Lat_b=Lat_b*Vasp_param
1188
     Lat_c=Lat_c*Vasp_param
1189

    
1190
     ALLOCATE(NbAtType(nat))
1191
     READ(IOIN,'(A)') Vasp_types
1192
     ! First, how many different types ?
1193
     NbAtType=0
1194
     READ(Vasp_types,*,END=998) NbAtType
1195
998  CONTINUE
1196
     NbType=0
1197
     DO WHILE (NbAtType(NbType+1).NE.0)
1198
        NbType=NbType+1
1199
     END DO
1200

    
1201
     ! Do we know the atom types ?
1202
     IF (AtTypes(1).EQ.'  ') THEN
1203
        ! user has not provided atom types... we have to find them ourselves
1204
        ! by looking into the POTCAR file...
1205
        INQUIRE(File="POTCAR", EXIST=Tchk)
1206
        IF (.NOT.Tchk) THEN
1207
           WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found"
1208
           STOP
1209
        END IF
1210
        OPEN(IOTMP,File="POTCAR")
1211
        DO I=1,NbType
1212
           Line='Empty'
1213
           DO WHILE (Line(1:2).NE.'US')
1214
              READ(IOTMP,'(A)') Line
1215
              Line=AdjustL(Line)
1216
           END DO
1217
           Line=adjustl(Line(3:))
1218
           AtTypes(I)=Line(1:2)
1219
        END DO
1220
        if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType)
1221
        CLOSE(IOTMP)
1222

    
1223
     ELSE  !AtTypes(1).EQ.'  '
1224
        ! user has provided atom types
1225
        NbTypeUser=0
1226
        DO WHILE (AtTypes(NbTypeUser+1).NE.'  ')
1227
           NbTypeUser=NbTypeUser+1
1228
        END DO
1229
        IF (NbType.NE.NbTypeUser) THEN
1230
           WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes"
1231
           STOP
1232
        END IF
1233
     END IF
1234

    
1235
     IAt=1
1236
     DO I=1,NbType
1237
        DO J=1,NbAtType(I)
1238
           AtName(Iat)=AtTypes(I)
1239
           Iat=Iat+1
1240
        END DO
1241
     END DO
1242
     DEALLOCATE(NbAtType)
1243

    
1244
     NbTypes=NbType
1245

    
1246
     READ(IOIN,'(A)') Vasp_comment
1247
     READ(IOIN,'(A)') Vasp_direct
1248
     V_direct=Adjustl(Vasp_direct)
1249
     Call UpCase(V_direct)
1250
! PFL 2011 Mar 8 ->
1251
! We have to read the FFF flags :
1252
     DO I=1,Nat
1253
        READ(IOIN,*) Xtmp,ytmp,ztmp,FFF(1:3,I)
1254
        DO J=1,3
1255
           FFF(J,I)=AdjustL(FFF(J,I))
1256
           CALL Upcase(FFF(J,I))
1257
        END DO
1258
     END DO
1259
! <- PFL 2011 Mar 8 
1260

    
1261
  END IF
1262

    
1263
  ! In the case of VASP there is always the problem  of moving from one side
1264
  ! of the box to the other...
1265
! PFL 2011 Mar 14 ->
1266
! In fact, we have to make this check as soon as either
1267
! Input and/or Prog = VASP
1268
! PFL 2010 Dec 2 ->
1269
! Here we deal with input and not prog in fact 
1270
!  if (prog.EQ.'VASP') THEN
1271
  if ((Input.EQ.'VASP').OR.(Prog.EQ.'VASP')) THEN
1272
! <- PFL 2010 Dec 2
1273
! <- PFL 2011 Mar 14
1274
     Renum=.TRUE.
1275

    
1276
     ! V_direct has been set in Read_geom
1277
     IF (V_direct(1:6).EQ.'DIRECT') THEN
1278
        Latr(1:3,1)=Lat_a
1279
        Latr(1:3,2)=Lat_b
1280
        Latr(1:3,3)=Lat_c
1281
        B=1.
1282
        CALL Gaussj(Latr,3,3,B,1,1)
1283
     ELSE
1284
        Latr=0.
1285
        Latr(1,1)=1.d0
1286
        Latr(2,2)=1.d0
1287
        Latr(3,3)=1.d0
1288
     END IF
1289

    
1290
     ! Actualization of Frozen using the FFFF... 
1291
     ! Frozen has the advantage ie if given, it imposes _ALL_ the FFF flags.
1292
     IF (Frozen(1).NE.0) THEN
1293
        WRITE(IOOUT,*) "Frozen is given. Flags of the given POSCAR are overriden"
1294
        FFF='T'
1295

    
1296
        NFroz=0
1297
        DO WHILE (Frozen(NFroz+1).NE.0)
1298
           NFroz=NFroz+1
1299
           FFF(1:3,Frozen(NFroz))='F'
1300
        END DO
1301
     ELSE
1302
        WRITE(IOOUT,*) "Frozen not given : using  Flags of the given POSCAR"
1303
        NFroz=0
1304
        Frozen=0
1305
        DO I=1,Nat
1306
           IF ((FFF(1,I).EQ.'F').OR.(FFF(2,I).EQ.'F').OR.(FFF(3,I).EQ.'F')) THEN
1307
              FFF(1:3,I)='F'
1308
              NFroz=NFroz+1
1309
              Frozen(NFroz)=I
1310
           END IF
1311
        END DO
1312
        WRITE(IOOUT,*) "Frozen atoms:",Frozen(1:NFroz)
1313
     END IF
1314

    
1315

    
1316
  END IF
1317

    
1318
  IF (Prog=="VASP") THEN
1319
     IF (Vmd) THEN
1320
        if (debug) WRITE(*,*) "DBG MAIN, L803, VMD=T,NbTypes,AtTypes",NbTypes,AtTypes(1:NbTypes)
1321
        Line=""
1322
        DO I=1,NbTypes
1323
           Line=TRIM(Line) // " " // TRIM(AdjustL(AtTypes(I))) 
1324
        END DO
1325
        Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title))
1326
        if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title)
1327
     END IF
1328
     Call CheckPeriodicBound
1329
  END IF
1330

    
1331
  IF (COORD=="MIXED") Call TestCart(AutoCart)
1332

    
1333
  SELECT CASE(NFroz)
1334
  CASE (2)
1335
     IntFroz=1
1336
  CASE (3)
1337
     IntFroz=3
1338
  CASE (4:)
1339
     IntFroz=(NFroz-2)*3  !(NFroz-3)*3+3
1340
  END SELECT
1341

    
1342
  if (debug) WRITE(*,'(1X,A,A,5I5)') "DBG Path, L825, Coord, NCart, NCoord, N3at, NFroz: "&
1343
       ,TRIM(COORD),Ncart,NCoord,N3at,NFroz
1344
  if (debug) WRITE(*,*) "DBG Path, L827, Cart",Cart
1345

    
1346
  if (((NCart+NFroz).EQ.0).AND.(Coord=="MIXED")) THEN
1347
     WRITE(IOOUT,*) "Coord=MIXED but Cart list empty: taking Coord=ZMAT."
1348
     COORD="ZMAT"
1349
  END IF
1350

    
1351
  IF ((Coord=="MIXED").AND.(NFroz.NE.0)) IntFroz=3*NFroz
1352

    
1353
  SELECT CASE(COORD)
1354
  CASE ('ZMAT')
1355
     NCoord=Nfree
1356
     ALLOCATE(dzdc(3,nat,3,nat))
1357
  CASE ('MIXED')
1358
     SELECT CASE (NCart+NFroz)
1359
     CASE (1) 
1360
        NCoord=N3At-3
1361
     CASE (2)
1362
        NCoord=N3At-1
1363
     CASE (3:)
1364
        NCoord=N3At
1365
     END SELECT
1366
     ALLOCATE(dzdc(3,nat,3,nat))
1367
  CASE ('BAKER')
1368
     Symmetry_elimination=0
1369
     NCoord=3*Nat-6-Symmetry_elimination !NFree = 3*Nat-6
1370
     ALLOCATE(BMat_BakerT(3*Nat,NCoord))
1371
     ALLOCATE(BTransInv(NCoord,3*Nat))
1372
     ALLOCATE(BTransInv_local(NCoord,3*Nat))
1373
  CASE ('HYBRID')
1374
     NCoord=N3at
1375
  CASE ('CART')
1376
     NCoord=N3at
1377
  END SELECT
1378

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

    
1381
  ALLOCATE(IntCoordI(NGeomI,NCoord),IntCoordF(NGeomF,NCoord))
1382
  ALLOCATE(XyzGeomF(NGeomF,3,Nat),XyzTangent(NGeomF,3*Nat))
1383
  ALLOCATE(Vfree(NCoord,NCoord),IntTangent(NGeomF,NCoord))
1384
  ALLOCATE(Energies(NgeomF),ZeroVec(NCoord)) ! Energies is declared in  Path_module.f90
1385
  ALLOCATE(Grad(NGeomF,NCoord),GeomTmp(NCoord),GradTmp(NCoord))
1386
  ALLOCATE(GeomOld_all(NGeomF,NCoord),GradOld(NGeomF,NCoord))
1387
  ALLOCATE(Hess(NGeomF,NCoord,NCoord),SGeom(NGeomF)) ! SGeom is declared in  Path_module.f90
1388
  ALLOCATE(EPathOld(NGeomF),MaxStep(NGeomF))
1389

    
1390
  ZeroVec=0.d0
1391
  DO I =1, NGeomF
1392
     IntTangent(I,:)=0.d0
1393
  END DO
1394

    
1395
  if (debug) THEN
1396
     Print *, 'L886, IntTangent(1,:)='
1397
     Print *, IntTangent(1,:)
1398
  END IF
1399

    
1400
  if (.NOT.OptReac) Energies(1)=Ereac
1401
  if (.NOT.OptProd) Energies(NGeomF)=EProd
1402
  MaxStep=SMax
1403

    
1404
  ! We initialize the mass array
1405
  if (MW) THEN
1406
     WRITE(*,*) 'Working in MW coordinates'
1407
     DO I=1,Nat
1408
        massAt(I)=Mass(Atome(I))
1409
     END DO
1410
  ELSE
1411
     DO I=1,Nat
1412
        MassAt(I)=1.d0
1413
     END DO
1414
  END IF
1415

    
1416
  WRITE(*,*) "Prog=",TRIM(PROG)
1417

    
1418
  ALLOCATE(FrozAtoms(Nat))
1419

    
1420
  !  IF (debug) WRITe(*,*) "DBG Main L940 NFroz",NFroz
1421
  !  IF (debug) WRITe(*,*) "DBG Main L940 Frozen",Frozen
1422
  !  IF (debug) WRITe(*,*) "DBG Main L940 FrozAtoms",FrozAtoms
1423
  IF (NFroz.EQ.0) THEN
1424
     FrozAtoms=.TRUE.
1425
  ELSE
1426
     I=1
1427
     NFr=0
1428
     FrozAtoms=.False.
1429
     DO WHILE (Frozen(I).NE.0)
1430
        IF (Frozen(I).LT.0) THEN
1431
           DO J=Frozen(I-1),abs(Frozen(I))
1432
              IF (.NOT.FrozAtoms(J)) THEN
1433
                 NFr=NFr+1
1434
                 FrozAtoms(J)=.TRUE.
1435
              END IF
1436
           END DO
1437
        ELSEIF (.NOT.FrozAtoms(Frozen(I))) THEN
1438
           FrozAtoms(Frozen(I))=.TRUE.
1439
           NFr=NFr+1
1440
        END IF
1441
        I=I+1
1442
     END DO
1443
     IF (NFr.NE.NFroz) THEN
1444
        WRITE(*,*) "Pb in PATH: Number of Frozen atoms not good :-( ",NFroz,NFr
1445
        STOP
1446
     END IF
1447
  END IF
1448

    
1449
  if (debug) THEN
1450
     !Print *, 'L968, IntTangent(1,:)='
1451
     !Print *, IntTangent(1,:)
1452
  END IF
1453

    
1454
  ! We have read everything we needed... time to work :-)
1455
  IOpt=0
1456
  FirstTimePathCreate = .TRUE. ! Initialized.
1457
  Call PathCreate(IOpt)    ! IntTangent is also computed in PathCreate.
1458
  FirstTimePathCreate = .FALSE.
1459

    
1460
  if (debug) THEN
1461
     Print *, 'L980, IntTangent(1,:)='
1462
     Print *, IntTangent(1,:)
1463
  END IF
1464

    
1465
  ! PFL 30 Mar 2008 ->
1466
  ! The following is now allocated in Calc_Baker. This is more logical to me
1467
  !   IF (COORD .EQ. "BAKER") Then
1468
  !      ALLOCATE(XprimitiveF(NgeomF,NbCoords))
1469
  !      ALLOCATE(UMatF(NGeomF,NbCoords,NCoord))
1470
  !      ALLOCATE(BTransInvF(NGeomF,NCoord,3*Nat))
1471
  !   END IF
1472
  ! <- PFL 30 mar 2008
1473

    
1474
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1475
  !
1476
  ! For Debugging purposes
1477
  !
1478
  ! OptGeom is the index of the geometry which is to be optimized.
1479
  IF (Optgeom.GT.0) THEN
1480
     Flag_Opt_Geom=.TRUE.
1481
     SELECT CASE(Coord)
1482
     CASE ('CART','HYBRID')
1483
!!! CAUTION : PFL 29.AUG.2008 ->
1484
        ! XyzGeomI/F stored in (NGeomI/F,3,Nat) and NOT (NGeomI/F,Nat,3)
1485
        ! This should be made  consistent !!!!!!!!!!!!!!!
1486
        GeomTmp=Reshape(reshape(XyzGeomF(OptGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/))
1487
        !           GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1488
!!! <- CAUTION : PFL 29.AUG.2008
1489
     CASE ('ZMAT','MIXED')
1490
        !GeomTmp=IntCoordF(OptGeom,:)
1491
        GeomTmp=IntCoordI(OptGeom,:)
1492
     CASE ('BAKER')
1493
        !GeomTmp=IntCoordF(OptGeom,:)
1494
        GeomTmp=IntCoordI(OptGeom,:)
1495
     CASE DEFAULT
1496
        WRITE(*,*) "Coord=",TRIM(COORD)," not recognized in Path L935.STOP"
1497
        STOP
1498
     END SELECT
1499

    
1500
     ! NCoord = NFree, Coord = Type of coordinate; e.g.: Baker
1501
     Flag_Opt_Geom = .TRUE.
1502
     Call Opt_geom(OptGeom,Nat,Ncoord,Coord,GeomTmp,E,Flag_Opt_Geom,Step_method)
1503

    
1504
     WRITE(Line,"('Geometry ',I3,' optimized')") Optgeom
1505
     WRITE(*,*) "Before PrintGeom, L1059, Path.f90"
1506
     Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,IoOut,Atome,Order,OrderInv,IndZmat)
1507
     STOP
1508
  END IF ! IF (Optgeom.GT.0) THEN
1509

    
1510
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1511
  ! End of GeomOpt
1512
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1513

    
1514
  IF (PathOnly) THEN
1515
     WRITE(*,*) "PathOnly=.T. , Stop here"
1516
     Call Write_path(-1)
1517
     if ((prog.EQ.'VASP').OR.WriteVasp) Call Write_vasp(poscar)
1518
     STOP
1519
  END IF
1520

    
1521
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1522

    
1523
  DEALLOCATE(XyzGeomI,IntCoordI)
1524
  NGeomI=NGeomF
1525
  ALLOCATE(XyzGeomI(NGeomF,3,Nat),IntCoordI(NGeomF,NCoord))
1526

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

    
1532
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1533
     DO IGeom=1,NGeomF
1534
        WRITE(*,*) "Before Calc_baker_allgeomf UMatF for IGeom=",IGeom
1535
        DO I=1,3*Nat-6
1536
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1537
        END DO
1538
     END DO
1539
  END IF
1540

    
1541
  ! To calculate B^T^-1 for all extrapolated final geometries:
1542
  ! PFL 18 Jan 2008 ->
1543
  ! Added a test for COORD=BAKER before the call to calc_baker_allgeomF
1544
  IF (COORD.EQ."BAKER")   Call Calc_baker_allGeomF()
1545
  ! <-- PFL 18 Jan 2008
1546

    
1547
  if (DEBUG.AND.(COORD.EQ."BAKER")) THEN
1548
     DO IGeom=1,NGeomF
1549
        WRITE(*,*) "After Calc_baker_allgeomf UMatF for IGeom=",IGeom
1550
        DO I=1,3*Nat-6
1551
           WRITE(*,'(50(1X,F12.8))') UMatF(IGeom,:,I)
1552
        END DO
1553
     END DO
1554
  END IF
1555

    
1556
  IOpt=0
1557
  Call EgradPath(IOpt,Flag_Opt_Geom)
1558

    
1559
  if (debug)  WRITE(*,*) "Calling Write_path, L1002, Path.f90, IOpt=",IOpt
1560
  Call Write_path(IOpt)
1561
  if ((debug.AND.(prog.EQ.'VASP')).OR.WriteVasp)  Call Write_vasp(poscar)
1562

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

    
1567
  ! By default, Hess=Id
1568
  Hess=0.
1569
  IF(Coord .EQ. "ZMAT") Then
1570
     ! We use the fact that we know that approximate good hessian values
1571
     ! are 0.5 for bonds, 0.1 for valence angle and 0.02 for dihedral angles
1572
     DO IGeom=1,NGeomF
1573
        Hess(IGeom,1,1)=0.5d0
1574
        Hess(IGeom,2,2)=0.5d0
1575
        Hess(IGeom,3,3)=0.1d0
1576
        DO J=1,Nat-3
1577
           Hess(IGeom,3*J+1,3*J+1)=0.5d0
1578
           Hess(IGeom,3*J+2,3*J+2)=0.1d0
1579
           Hess(IGeom,3*J+3,3*J+3)=0.02d0
1580
        END DO
1581
        IF (HInv) THEN
1582
           DO I=1,NCoord
1583
              Hess(IGeom,I,I)=1.d0/Hess(IGeom,I,I)
1584
           END DO
1585
        END IF
1586
     END DO
1587
  ELSE
1588
     DO I=1,NCoord
1589
        DO J=1,NGeomF
1590
           Hess(J,I,I)=1.d0
1591
        END DO
1592
     END DO
1593
  END IF
1594

    
1595
  IF (COORD .EQ. "BAKER") THEN
1596
     ! UMat(NPrim,NCoord)
1597
     ALLOCATE(Hprim(NPrim,NPrim),HprimU(NPrim,3*Nat-6))
1598
     Hprim=0.d0
1599
     ScanCoord=>Coordinate
1600
     I=0
1601
     DO WHILE (Associated(ScanCoord%next))
1602
        I=I+1
1603
        SELECT CASE (ScanCoord%Type)
1604
        CASE ('BOND')
1605
           !DO J=1,NGeomF ! why all geometries?? Should be only 1st and NGeomF??
1606
           Hprim(I,I) = 0.5d0
1607
           !END DO
1608
        CASE ('ANGLE')
1609
           Hprim(I,I) = 0.2d0
1610
        CASE ('DIHEDRAL')
1611
           Hprim(I,I) = 0.1d0
1612
        END SELECT
1613
        ScanCoord => ScanCoord%next
1614
     END DO
1615

    
1616
     ! HprimU = H^prim U:
1617
     HprimU=0.d0
1618
     DO I=1, 3*Nat-6
1619
        DO J=1,NPrim
1620
           HprimU(:,I) = HprimU(:,I) + Hprim(:,J)*UMat(J,I)
1621
        END DO
1622
     END DO
1623

    
1624
     ! Hess = U^T HprimU = U^T H^prim U:
1625
     Hess=0.d0
1626
     DO I=1, 3*Nat-6
1627
        DO J=1,NPrim
1628
           ! UMat^T is needed below.
1629
           Hess(1,:,I) = Hess(1,:,I) + UMat(J,:)*HprimU(J,I)
1630
        END DO
1631
     END DO
1632
     DO K=2,NGeomF
1633
        Hess(K,:,:)=Hess(1,:,:)
1634
     END DO
1635
     DEALLOCATE(Hprim,HprimU)
1636
  end if !  ends IF COORD=BAKER
1637
  ALLOCATE(GradTmp2(NCoord),GeomTmp2(NCoord))
1638
  ALLOCATE(HessTmp(NCoord,NCoord))
1639

    
1640
  IF (IniHUp) THEN
1641
     IF (FCalcHess) THEN
1642
        ! The numerical calculation of the Hessian has been switched on
1643
        DO IGeom=1,NGeomF
1644
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1645
              GeomTmp=IntCoordF(IGeom,:)
1646
           ELSE
1647
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
1648
           END IF
1649
           Call CalcHess(NCoord,GeomTmp,Hess(IGeom,:,:),IGeom,Iopt)
1650
        END DO
1651
     ELSE
1652
        ! We generate a forward hessian and a backward one and we mix them.
1653
        ! First, the forward way:
1654
        DO IGeom=2,NGeomF-1
1655
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1656
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
1657
           ELSE
1658
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
1659
           END IF
1660
           GradTmp=Grad(IGeom-1,:)-Grad(IGeom,:)
1661
           HessTmp=Hess(IGeom-1,:,:)
1662
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1663
           Hess(IGeom,:,:)=HessTmp
1664
        END DO
1665

    
1666
        ! Now, the backward way:
1667
        HessTmp=Hess(NGeomF,:,:)
1668
        DO IGeom=NGeomF-1,2,-1
1669
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1670
              GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
1671
           ELSE
1672
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
1673
           END IF
1674
           GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
1675
           Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
1676
           alpha=0.5d0*Pi*(IGeom-1.)/(NGeomF-1.)
1677
           ca=cos(alpha)
1678
           sa=sin(alpha)
1679
           Hess(IGeom,:,:)=ca*Hess(IGeom,:,:)+sa*HessTmp(:,:)
1680
        END DO
1681
     END IF ! matches IF (FCalcHess)
1682
  END IF ! matches IF (IniHUp) THEN
1683

    
1684
  ! For debug purposes, we diagonalize the Hessian matrices
1685
  IF (debug) THEN
1686
     !WRITE(*,*) "DBG Main, L1104, - EigenSystem for Hessian matrices"
1687
     DO I=1,NGeomF
1688
        WRITE(*,*) "DBG Main - Point ",I
1689
        Call Jacobi(Hess(I,:,:),NCoord,GradTmp2,HessTmp,NCoord)
1690
        DO J=1,NCoord
1691
           WRITE(*,'(1X,I3,1X,F10.6," : ",20(1X,F8.4))') J,GradTmp2(J),HessTmp(J,1:min(20,NCoord))
1692
        END DO
1693
     END DO
1694
  END IF ! matches IF (debug) THEN
1695

    
1696
  ! we have the hessian matrices and the gradients, we can start the optimizations !
1697
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1698
  !
1699
  ! Main LOOP : Optimization of the Path
1700
  !
1701
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1702
  if (debug)   Print *, 'NGeomF=', NGeomF
1703
  allocation_flag = .TRUE.
1704

    
1705
  Fini=.FALSE.
1706

    
1707
  DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT. Fini))
1708
     if (debug) Print *, 'IOpt at the begining of the loop, L1128=', IOpt
1709
     IOpt=IOpt+1
1710

    
1711
     SELECT CASE (COORD)
1712
     CASE ('ZMAT','MIXED','BAKER')
1713
        GeomOld_all=IntCoordF
1714
     CASE ('CART','HYBRID')
1715
        GeomOld_all=Reshape(XyzGeomF,(/NGeomF,NCoord/))
1716
     CASE DEFAULT
1717
        WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1149.STOP'
1718
        STOP
1719
     END SELECT
1720

    
1721
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1722
          valid("printtangent")) THEN
1723
        WRITE(*,*) "Visualization of tangents"
1724
        Idx=min(12,NCoord)
1725
        DO I=1,NGeomF
1726
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-0.5d0*IntTangent(I,1:Idx)
1727
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+0.5d0*IntTangent(I,1:Idx)
1728
           WRITE(*,*) 
1729
        END DO
1730
        WRITE(*,*) "END of tangents"
1731
     END IF
1732

    
1733
     IF (((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER').OR.(COORD.EQ.'MIXED')).AND. &
1734
          valid("printgrad")) THEN
1735
        WRITE(*,*) "Visualization of gradients"
1736
        Idx=min(12,NCoord)
1737
        DO I=1,NGeomF
1738
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)-1.d0*Grad(I,1:Idx)
1739
           WRITE(*,'(12(1X,F12.5))') IntCoordF(I,1:Idx)+1.d0*Grad(I,1:Idx)
1740
           WRITe(*,*) 
1741
        END DO
1742
        WRITE(*,*) "END of gradients"
1743
     END IF
1744

    
1745
     Fini=.TRUE.
1746
     IF (OptReac) THEN !OptReac: True if one wants to optimize the geometry of the reactants. Default is FALSE.
1747
        IGeom=1
1748
        if (debug) WRITE(*,*) '**************************************'
1749
        if (debug) WRITE(*,*) 'Optimizing reactant'
1750
        if (debug) WRITE(*,*) '**************************************'
1751
        SELECT CASE (COORD)
1752
        CASE ('ZMAT','MIXED','BAKER')
1753
           SELECT CASE (Step_method)
1754
           CASE ('RFO')
1755
              !GradTmp(NCoord) becomes Step in Step_RFO_all and has INTENT(OUT).
1756
              Call Step_RFO_all(NCoord,GradTmp,1,IntCoordF(1,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1757
              Print *, GradTmp
1758
           CASE ('GDIIS')
1759
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1760
              Call Step_DIIS_all(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),HP,HEAT,&
1761
                   Hess(1,:,:),NCoord,allocation_flag,ZeroVec)
1762
              Print *, 'OptReac, Steps from GDIIS, GeomNew - IntCoordF(1,:)='
1763
              Print *, GradTmp
1764
           CASE ('GEDIIS')
1765
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT).:
1766
              ! Energies are updated in EgradPath.f90
1767
              Call Step_GEDIIS_All(NGeomF,1,GradTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),&
1768
                   NCoord,allocation_flag,ZeroVec)
1769
              !Call Step_GEDIIS(GeomTmp,IntCoordF(1,:),Grad(1,:),Energies(1),Hess(1,:,:),NCoord,allocation_flag)
1770
              !GradTmp = GeomTmp - IntCoordF(1,:)
1771
              !Print *, 'Old Geometry:'
1772
              !Print *, IntCoordF(1,:)
1773
              Print *, 'OptReac, Steps from GEDIIS, GradTmp = GeomTmp - IntCoordF(1,:)='
1774
              Print *, GradTmp
1775
              !Print *, 'GeomTmp:'
1776
              !Print *, GeomTmp
1777
              GeomTmp(:) = GradTmp(:) + IntCoordF(1,:)
1778
           CASE DEFAULT
1779
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1194. Stop"
1780
              STOP
1781
           END SELECT
1782
        CASE ('HYBRID','CART')
1783
           Call Step_RFO_all(NCoord,GradTmp,1,XyzGeomF(1,:,:),Grad(1,:),Hess(1,:,:),ZeroVec)
1784
        CASE DEFAULT
1785
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1200. STOP'
1786
           STOP
1787
        END SELECT
1788

    
1789
        IF (debug) THEN
1790
           WRITE(Line,*) 'DBG Path, L1227,', TRIM(Step_method), 'Step for IOpt=',IOpt, ', IGeom=1'
1791
           Call PrintGeom(Line,Nat,NCoord,GeomTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1792
        END IF
1793

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

    
1807
        GradTmp=GradTmp*FactStep
1808

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

    
1829
        IF (debug) THEN
1830
           !WRITE(Line,*) 'DBG Path, L1244, Step for IOpt=',IOpt, ', IGeom=1'
1831
           !Call PrintGeom(Line,Nat,NCoord,GradTmp,Coord,6,Atome,Order,OrderInv,IndZmat)
1832
        END IF
1833

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

    
1846
        SELECT CASE (COORD)
1847
        CASE ('ZMAT','MIXED','BAKER')
1848
           IntcoordF(1,:)=IntcoordF(1,:)+GradTmp(:) !IGeom=1
1849
        CASE ('HYBRID','CART')
1850
           XyzGeomF(1,:,:)=XyzGeomF(1,:,:)+reshape(GradTmp(:),(/3,Nat/))
1851
        CASE DEFAULT
1852
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1201.STOP'
1853
           STOP
1854
        END SELECT
1855

    
1856
        IF (debug) THEN
1857
           WRITE(*,*) "DBG Main, L1271, IOpt=",IOpt, ', IGeom=1'
1858
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1859
              WRITE(*,'(12(1X,F8.3))') IntCoordF(1,1:min(12,NFree))
1860
           ELSE
1861
              DO Iat=1,Nat
1862

    
1863
!!!!!!!!!! There was an OrderInv here... should I put it back ?
1864
                 ! It was with indzmat. IndZmat cannot appear here...
1865
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1866
                      XyzGeomF(IGeom,1:3,Iat)
1867
              END DO
1868
           END IF
1869
        END IF ! matches IF (debug) THEN
1870

    
1871
        IF (.NOT.OptReac) WRITE(*,*) "Reactants optimized"
1872
     ELSE ! Optreac
1873
        IF (debug) WRITE(*,*) "Reactants already optimized, Fini=",Fini
1874
     END IF ! matches IF (OptReac) THEN
1875

    
1876
     IF (OptProd) THEN !OptProd: True if one wants to optimize the geometry of the products. Default is FALSE.
1877
        IGeom=NGeomF
1878
        if (debug) WRITE(*,*) '******************'
1879
        if (debug) WRITE(*,*) 'Optimizing product'
1880
        if (debug) WRITE(*,*) '******************'
1881
        SELECT CASE (COORD)
1882
        CASE ('ZMAT','MIXED','BAKER')
1883
           Print *, 'Optimizing product'
1884
           SELECT CASE (Step_method)
1885
           CASE ('RFO')
1886
              !GradTmp(NCoord)) becomes "Step" in Step_RFO_all and has INTENT(OUT).
1887
              Call Step_RFO_all(NCoord,GradTmp,NGeomF,IntCoordF(NGeomF,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1888
              Print *, GradTmp
1889
           CASE ('GDIIS')
1890
              ! GradTmp becomes "step" below and has INTENT(OUT):
1891
              Call Step_DIIS_all(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),&
1892
                   HP,HEAT,Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1893
              Print *, 'OptProd, Steps from GDIIS, GeomNew - IntCoordF(NGeomF,:)='
1894
              Print *, GradTmp
1895
           CASE ('GEDIIS')
1896
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
1897
              Call Step_GEDIIS_All(NGeomF,NGeomF,GradTmp,IntCoordF(NGeomF,:),Grad(NGeomF,:),Energies(NGeomF),&
1898
                   Hess(NGeomF,:,:),NCoord,allocation_flag,ZeroVec)
1899
              Print *, 'OptProd, Steps from GEDIIS, GeomNew - IntCoordF(NGeomF,:)='
1900
              Print *, GradTmp
1901
           CASE DEFAULT
1902
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1309. Stop"
1903
              STOP
1904
           END SELECT
1905
        CASE ('HYBRID','CART')
1906
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(NGeomF,:,:),Grad(NGeomF,:),Hess(NGeomF,:,:),ZeroVec)
1907
        CASE DEFAULT
1908
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1460.STOP'
1909
           STOP
1910
        END SELECT
1911

    
1912
        ! we have the Newton+RFO step, we will check that the displacement is not greater than Smax
1913
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
1914
        FactStep=min(1.d0,MaxStep(IGeom)/NormStep)
1915
        Fini=(Fini.AND.(NormStep.LE.SThresh))
1916
        OptProd=.NOT.(NormStep.LE.SThresh)
1917
        NormGrad=sqrt(DOT_Product(reshape(Grad(NGeomF,:),(/Nfree/)),reshape(Grad(NGeomF,:),(/NFree/))))
1918
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
1919
        OptProd=(OptProd.OR.(NormGrad.GT.GThresh))
1920
        IF (debug) WRITE(*,*) "NormGrad, GThresh, OptProd=",NormGrad, GThresh, OptProd
1921

    
1922
        GradTmp=GradTmp*FactStep
1923

    
1924
        ! we take care that frozen atoms do not move
1925
        IF (NFroz.NE.0) THEN
1926
           SELECT CASE (COORD)
1927
           CASE ('ZMAT','MIXED','BAKER')
1928
              GradTmp(1:IntFroz)=0.d0
1929
           CASE ('CART','HYBRID')
1930
              DO I=1,Nat
1931
                 IF (FrozAtoms(I)) THEN
1932
                    Iat=Order(I)
1933
                    GradTmp(3*Iat-2:3*Iat)=0.d0
1934
                 END IF
1935
              END DO
1936
           CASE DEFAULT
1937
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1045.STOP'
1938
              STOP
1939
           END SELECT
1940
        END IF ! matches IF (NFroz.NE.0) THEN
1941

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

    
1954
        SELECT CASE (COORD)
1955
        CASE ('ZMAT','MIXED','BAKER')
1956
           IntcoordF(NGeomF,:)=IntcoordF(NGeomF,:)+GradTmp(:) ! IGeom is NGeomF.
1957
        CASE ('HYBRID','CART')
1958
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
1959
        CASE DEFAULT
1960
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1050.STOP'
1961
           STOP
1962
        END SELECT
1963

    
1964
        IF (debug) THEN
1965
           WRITE(*,*) "DBG Main, L1376, IGeom=, IOpt=",IGeom,IOpt
1966
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
1967
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
1968
           ELSE
1969
              DO Iat=1,Nat
1970
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
1971
                      XyzGeomF(IGeom,1:3,Iat)
1972
              END DO
1973
           END IF
1974
        END IF
1975
     ELSE ! Optprod
1976
        IF (debug) WRITE(*,*) "Product already optimized, Fini=",Fini
1977
     END IF !matches IF (OptProd) THEN 
1978

    
1979
     IF (debug) WRITE(*,*) "Fini before L1491 path:",Fini
1980

    
1981
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1982
     !
1983
     !  Optimizing other geometries
1984
     !
1985
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1986

    
1987

    
1988

    
1989
     DO IGeom=2,NGeomF-1 ! matches at L1556
1990
        if (debug) WRITE(*,*) '****************************'
1991
        if (debug) WRITE(*,*) 'All other geometries, IGeom=',IGeom
1992
        if (debug) WRITE(*,*) '****************************'
1993

    
1994
        SELECT CASE (COORD)
1995
        CASE ('ZMAT','BAKER','MIXED')
1996
           GradTmp2=reshape(IntTangent(IGeom,:),(/NCoord/))
1997
           ! PFL 6 Apr 2008 ->
1998
           ! Special case, if FTan=0. we do not consider Tangent at all.
1999
           ! To DO: add the full treatment of FTan
2000
           if (FTan==0.) GradTmp2=ZeroVec
2001
           ! <- PFL 6 Apr 2008
2002
           if (debug) THEN
2003
              Print *, 'L1500, IntTangent(',IGeom,',:)='
2004
              Print *, IntTangent(IGeom,:)
2005
           END IF
2006
           !GradTmp(NCoord)) is actually "Step" in Step_RFO_all and has INTENT(OUT).
2007
           SELECT CASE (Step_method)
2008
           CASE ('RFO')
2009
              Call Step_RFO_all(NCoord,GradTmp,IGeom,IntCoordF(IGeom,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
2010
           CASE ('GDIIS')
2011
              !The following has no effect at IOpt==1
2012
              !Print *, 'Before call of Step_diis_all, All other geometries, IGeom=', IGeom
2013
              Call Step_DIIS_all(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),HP,HEAT,&        
2014
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
2015
              Print *, 'All others, Steps from GDIIS, GeomNew - IntCoordF(IGeom,:)='
2016
              Print *, GradTmp
2017
           CASE ('GEDIIS')
2018
              ! GradTmp(NCoord) becomes "step" below and has INTENT(OUT):
2019
              Call Step_GEDIIS_All(NGeomF,IGeom,GradTmp,IntCoordF(IGeom,:),Grad(IGeom,:),Energies(IGeom),&
2020
                   Hess(IGeom,:,:),NCoord,allocation_flag,GradTmp2)
2021
              Print *, 'All others, Steps from GEDIIS, GeomNew - IntCoordF(IGeom,:)='
2022
              Print *, GradTmp
2023
           CASE DEFAULT
2024
              WRITE (*,*) "Step_method=",TRIM(Step_method)," not recognized, Path.f90, L1430. Stop"
2025
              STOP
2026
           END SELECT
2027
        CASE ('HYBRID','CART')
2028
           ! XyzTangent is implicitely in Nat,3... but all the rest is 3,Nat
2029
           ! so we change it
2030
           DO I=1,Nat
2031
              DO J=1,3
2032
                 GradTmp2(J+(I-1)*3)=XyzTangent(IGeom,(J-1)*Nat+I)
2033
              END DO
2034
           END DO
2035
           !           GradTmp2=XyzTangent(IGeom,:)
2036
           ! PFL 6 Apr 2008 ->
2037
           ! Special case, if FTan=0. we do not consider Tangent at all.
2038
           ! To DO: add the full treatment of FTan
2039
           if (FTan==0.) GradTmp2=ZeroVec
2040
           ! <- PFL 6 Apr 2008
2041

    
2042
           Call Step_RFO_all(NCoord,GradTmp,IGeom,XyzGeomF(IGeom,:,:),Grad(IGeom,:),Hess(IGeom,:,:),GradTmp2)
2043
        CASE DEFAULT
2044
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1083.STOP'
2045
           STOP
2046
        END SELECT
2047

    
2048
        ! we take care that frozen atoms do not move
2049
        IF (NFroz.NE.0) THEN
2050
           SELECT CASE (COORD)
2051
           CASE ('ZMAT','MIXED','BAKER')
2052
              IF (debug) THEN 
2053
                 WRITE(*,*) 'Step computed. Coord=',Coord
2054
                 WRITE(*,*) 'Zeroing components for frozen atoms. Intfroz=', IntFroz
2055
              END IF
2056
              GradTmp(1:IntFroz)=0.d0
2057
           CASE ('CART','HYBRID')
2058
              if (debug) THEN
2059
                 WRITE(*,*) "DBG Path Zeroing components L1771. Step before zeroing"
2060
                 DO I=1,Nat
2061
                    WRITe(*,*) I,GradTmp(3*I-2:3*I)
2062
                 END DO
2063
              END IF
2064
              DO I=1,Nat
2065
                 IF (FrozAtoms(I)) THEN
2066
                    if (debug) THEN
2067
                       write(*,*) 'Step Computed. Coord=',Coord
2068
                       WRITE(*,*) 'Atom',I,' is frozen. Zeroing',Order(I)
2069
                    END IF
2070
                    Iat=Order(I)
2071
                    GradTmp(3*Iat-2:3*Iat)=0.d0
2072
                 END IF
2073
              END DO
2074

    
2075
                 if (debug) THEN
2076
                    WRITE(*,*) "DBG Path Zeroing components L1785. Step After zeroing"
2077
                    DO I=1,Nat
2078
                       WRITe(*,*) I,GradTmp(3*I-2:3*I)
2079
                    END DO
2080
                 END IF
2081

    
2082
           CASE DEFAULT
2083
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1338.STOP'
2084
              STOP
2085
           END SELECT
2086
        END IF !matches IF (NFroz.NE.0) THEN
2087

    
2088
        IF (debug) THEN
2089
           SELECT CASE (COORD)
2090
           CASE ('ZMAT','MIXED','BAKER')
2091
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',IntTangent(IGeom,1:min(12,NFree))
2092
              WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
2093
           CASE ('CART','HYBRID')
2094
              WRITE(*,'(A,12(1X,F8.3))') 'Step tan:',XyzTangent(IGeom,1:min(12,NFree))
2095
               WRITE(*,'(A,12(1X,F8.3))') 'Ortho Step:',GradTmp(1:min(12,NFree))
2096
           CASE DEFAULT
2097
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1588.STOP'
2098
              STOP
2099
           END SELECT
2100
        END IF ! matches if (debug) THEN
2101

    
2102
        ! we check that the maximum displacement is not greater than Smax
2103
        If (debug) WRITE(*,*) "Fini before test:",Fini
2104
        NormStep=sqrt(DOT_Product(GradTmp,GradTmp))      
2105
        FactStep=min(1.d0,maxStep(Igeom)/NormStep)
2106
        Fini=(Fini.AND.(NormStep.LE.SThresh))
2107
        IF (debug) WRITE(*,*) "IGeom,NormStep, SThresh,Fini",IGeom,NormStep, SThresh, Fini
2108

    
2109
        GradTmp=GraDTmp*FactStep
2110

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

    
2114
        IF (debug) THEN
2115
           WRITE(*,*) "DBG Main, L1479, IGeom=, IOpt=",IGeom,IOpt
2116
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2117
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
2118
           ELSE
2119
              DO Iat=1,Nat
2120
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
2121
                      XyzGeomF(IGeom,1:3,Iat)
2122
              END DO
2123
           END IF
2124
        END IF ! matches if (debug) THEN
2125

    
2126
        SELECT CASE (COORD)
2127
        CASE ('ZMAT')
2128
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2129
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2130
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2131
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2132
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2133
                OrderInv(indzmat(2,2)),IntcoordF(IGeom,1)
2134
           WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2135
                OrderInv(indzmat(3,2)),IntcoordF(IGeom,2), OrderInv(indzmat(3,3)),IntcoordF(IGeom,3)*180./Pi
2136
           Idx=4
2137
           DO Iat=4,Nat
2138
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2139
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2140
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2141
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2142
              Idx=Idx+3
2143
           END DO
2144

    
2145
        CASE ('BAKER')
2146
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2147
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2148
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2149
           WRITE(IOOUT,*) "COORD=BAKER, printing IntCoordF is not meaningfull."
2150

    
2151
        CASE ('MIXED')
2152
           IntcoordF(IGeom,:)=IntcoordF(IGeom,:)+GradTmp(:)
2153
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2154
           WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2155
           DO Iat=1,NCart
2156
              WRITE(IOOUT,'(1X,A5,3(1X,F13.8))') Nom(Atome(OrderInv(indzmat(1,1)))),IntcoordF(IGeom,3*Iat-2:3*Iat)
2157
           END DO
2158
           Idx=3*NCart+1
2159
           IF (NCart.EQ.1) THEN
2160
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2161
                   OrderInv(indzmat(2,2)),IntcoordF(IGeom,4)
2162
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2163
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,5), OrderInv(indzmat(3,3)),IntcoordF(IGeom,6)*180./Pi
2164

    
2165
              Idx=7
2166
           END IF
2167
           IF (NCart.EQ.2) THEN
2168
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2169
                   OrderInv(indzmat(3,2)),IntcoordF(IGeom,7), OrderInv(indzmat(3,3)),IntcoordF(IGeom,8)*180./Pi
2170
              Idx=9
2171
           END IF
2172

    
2173
           
2174
           DO Iat=max(NCart+1,4),Nat
2175
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2176
                   OrderInv(indzmat(iat,2)),IntcoordF(IGeom,Idx), &
2177
                   OrderInv(indzmat(iat,3)),IntcoordF(IGeom,Idx+1)*180./pi, &
2178
                   OrderInv(indzmat(iat,4)),IntcoordF(IGeom,Idx+2)*180/Pi
2179
              Idx=Idx+3
2180
           END DO
2181
           if (debug) Print *, 'New Geom=IntcoordF(',IGeom,',:)=', IntcoordF(IGeom,:)
2182

    
2183
        CASE ('HYBRID','CART')
2184
           XyzGeomF(IGeom,:,:)=XyzGeomF(IGeom,:,:)+reshape(GradTmp(:),(/3,Nat/))
2185
        CASE DEFAULT
2186
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1537.STOP'
2187
           STOP
2188
        END SELECT
2189

    
2190
        IF (debug) THEN
2191
           WRITE(*,*) "DBG Main, L1516, IGeom=, IOpt=",IGeom,IOpt
2192
           IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2193
              WRITE(*,'(12(1X,F8.3))') IntCoordF(IGeom,1:min(12,NFree))
2194
           ELSE
2195
              DO Iat=1,Nat
2196
                 WRITE(*,'(1X,A5,3(1X,F11.6))') Nom(Atome(Iat)), &
2197
                      XyzGeomF(IGeom,1:3,Iat)
2198
              END DO
2199
           END IF
2200
        END IF ! matches if (debug) THEN
2201

    
2202
        ! We project out the tangential part of the gradient to know if we are done
2203
        GradTmp=Grad(IGeom,:)
2204
        NormGrad=sqrt(dot_product(GradTmp,GradTmp))
2205
        if (debug) WRITE(*,*) "Norm Grad tot=",NormGrad
2206
        SELECT CASE (COORD)
2207
        CASE ('ZMAT','MIXED','BAKER')
2208
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2209
           Norm=DOT_PRODUCT(IntTangent(IGeom,:),IntTangent(IGeom,:))
2210
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2211
           Ca=S/(sqrt(Norm)*NormGrad)
2212
           S=DOT_PRODUCT(GradTmp,IntTangent(IGeom,:))
2213
           GradTmp=GradTmp-S/Norm*IntTangent(IGeom,:)
2214
        CASE ('CART','HYBRID')
2215
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2216
           Norm=DOT_PRODUCT(XyzTangent(IGeom,:),XyzTangent(IGeom,:))
2217
           Ca=S/(sqrt(Norm)*NormGrad)
2218
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2219
           S=DOT_PRODUCT(GradTmp,XyzTangent(IGeom,:))
2220
           GradTmp=GradTmp-S/Norm*XyzTangent(IGeom,:)
2221
        CASE DEFAULT
2222
           WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1190.STOP'
2223
           STOP
2224
        END SELECT
2225
        IF (debug) WRITE(*,'(1X,A,3(1X,F10.6))') "cos and angle between grad and tangent",ca, acos(ca), acos(ca)*180./pi
2226
        NormGrad=sqrt(DOT_Product(GradTmp,GradTmp))
2227
        Fini=(Fini.AND.(NormGrad.LE.GThresh))
2228
        IF (debug) WRITE(*,*) "NormGrad_perp, GThresh, Fini=",NormGrad, GThresh, Fini
2229

    
2230
     END DO ! matches DO IGeom=2,NGeomF-1
2231
     ! we save the old gradients
2232
     GradOld=Grad
2233
     EPathOld=Energies
2234

    
2235
     ! Shall we re-parameterize the path ?
2236
     ! PFL 2007/Apr/10 ->
2237
     ! We call PathCreate in all cases, it will take care of the 
2238
     ! reparameterization, as well as calculating the tangents.
2239
     !     IF (MOD(IOpt,IReparam).EQ.0) THEN
2240
     XyzGeomI=XyzGeomF
2241
     IntCoordI=IntCoordF
2242

    
2243
     Call PathCreate(IOpt)
2244
     !     END IF
2245
     ! <- PFL 2007/Apr/10
2246

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

    
2250
        ! We have the new path, we calculate its energy and gradient
2251

    
2252
        Call EgradPath(IOpt,Flag_Opt_Geom)
2253
        !IF(IOPT .EQ. 10) Then
2254
        !         Print *, 'Stopping at my checkpoint.'
2255
        !           STOP !This is my temporary checkpoint.
2256
        !ENDIF
2257

    
2258
        ! PFL 31 Mar 2008 ->
2259
        ! Taken from v3.99 but we added a flag...
2260
        ! Added in v3.99 : MaxStep is modified according to the evolution of energies
2261
        ! if Energies(IGeom)<=EPathOld then MaxStep is increased by 1.1 
2262
        ! else it is decreased by 0.8
2263

    
2264
        if (DynMaxStep) THEN
2265
           if (debug) WRITE(*,*) "Dynamically updating the Maximum step"
2266
           if (OptReac) THEN
2267
              If (Energies(1)<EPathOld(1)) Then
2268
                 MaxStep(1)=min(1.1*MaxStep(1),2.*SMax)
2269
              ELSE
2270
                 MaxStep(1)=max(0.8*MaxStep(1),SMax/2.)
2271
              END IF
2272
           END IF
2273

    
2274
           if (OptProd) THEN
2275
              If (Energies(NGeomF)<EPathOld(NGeomF)) Then
2276
                 MaxStep(NGeomF)=min(1.1*MaxStep(NGeomF),2.*SMax)
2277
              ELSE
2278
                 MaxStep(NGeomF)=max(0.8*MaxStep(NGeomF),SMax/2.)
2279
              END IF
2280
           END IF
2281

    
2282
           DO IGeom=2,NGeomF-1
2283
              If (Energies(IGeom)<EPathOld(IGeom)) Then
2284
                 MaxStep(IGeom)=min(1.1*MaxStep(IGeom),2.*SMax)
2285
              ELSE
2286
                 MaxStep(IGeom)=max(0.8*MaxStep(IGeom),SMax/2.)
2287
              END IF
2288
           END DO
2289
           if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2290
        END IF ! test on DynMaxStep
2291
        if (debug) WRITE(*,*) "Maximum step",MaxStep(1:NgeomF)
2292
        ! <- PFL 31 MAr 2008
2293
        ! Also XyzGeomF is updated in EgradPath for Baker Case.
2294
        !  It should get updated for other cases too.
2295

    
2296
        DO IGeom=1,NGeomF
2297
           SELECT CASE (COORD)
2298
           CASE('ZMAT')
2299
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2300
              GeomTmp=IntCoordF(IGeom,:)
2301
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(1,1))))
2302
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(2,1)))), &
2303
                   OrderInv(indzmat(2,2)),Geomtmp(1)
2304
              WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(3,1)))), &
2305
                   OrderInv(indzmat(3,2)),Geomtmp(2), &
2306
                   OrderInv(indzmat(3,3)),Geomtmp(3)*180./Pi
2307
              Idx=4
2308
              DO Iat=4,Nat
2309
                 WRITE(IOOUT,'(1X,A5,3(1X,I4,1X,F11.6))') Nom(Atome(OrderInv(indzmat(iat,1)))), &
2310
                      OrderInv(indzmat(iat,2)),Geomtmp(Idx), &
2311
                      OrderInv(indzmat(iat,3)),Geomtmp(Idx+1)*180./pi, &
2312
                      OrderInv(indzmat(iat,4)),Geomtmp(Idx+2)*180/Pi
2313
                 Idx=Idx+3
2314
              END DO
2315
           CASE('BAKER')
2316
              !WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2317
              !GeomTmp=IntCoordF(IGeom,:)
2318
           CASE('CART','HYBRID','MIXED')
2319
              WRITE(IOOUT,'(1X,I5)') Nat
2320
              WRITE(IOOUT,"('Geometry ',I3,'/',I3,' for iteration ',I3)") IGeom,NgeomF,IOpt
2321
              GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2322
              DO I=1,Nat
2323
                 Iat=I
2324
                 If (renum) Iat=OrderInv(I)
2325
                 WRITE(IOOUT,'(1X,A5,3(1X,F11.6))') Nom(Atome(iat)), Geomtmp(3*Iat-2:3*Iat)
2326
              END DO
2327
           CASE DEFAULT
2328
              WRITE(*,*) 'Coord=',TRIM(COORD),' not recognized. Path L1356.STOP'
2329
              STOP
2330
           END SELECT
2331
        END DO ! matches DO IGeom=1,NGeomF
2332

    
2333
        Call Write_path(IOpt)
2334

    
2335
        ! We update the Hessian matrices
2336
        IF (debug) WRITE(*,*) "Updating Hessian matrices",NCoord
2337
        ! First using the displacement from the old path
2338
        IGeom0=2
2339
        IGeomF=NGeomF-1
2340
        IF (OptReac) IGeom0=1
2341
        IF (OptProd) IGeomF=NGeomF
2342

    
2343
        IF (FCalcHess) THEN
2344
           DO IGeom=IGeom0,IGeomF
2345
              SELECT CASE (COORD)
2346
              CASE ('ZMAT','MIXED','BAKER')
2347
                 GeomTmp2=IntCoordF(IGeom,:)
2348
              CASE ('CART','HYBRID')
2349
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2350
              CASE DEFAULT
2351
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2352
                 STOP
2353
              END SELECT
2354
              Call CalcHess(NCoord,GeomTmp2,Hess(IGeom,:,:),IGeom,Iopt)
2355
           END DO
2356
        ELSE
2357
           DO IGeom=IGeom0,IGeomF
2358
              GeomTmp=GeomOld_all(IGeom,:)
2359
              SELECT CASE (COORD)
2360
              CASE ('ZMAT','MIXED','BAKER')
2361
                 GeomTmp2=IntCoordF(IGeom,:)
2362
              CASE ('CART','HYBRID')
2363
                 GeomTmp2=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))
2364
              CASE DEFAULT
2365
                 WRITE(*,*) "Coord=",TRIM(COORD)," not recognized. PATH L1267.STOP"
2366
                 STOP
2367
              END SELECT
2368
              GeomTmp=GeomTmp2-GeomTmp
2369
              GradTmp=Grad(IGeom,:)-GradOld(IGeom,:)
2370
              HessTmp=Hess(IGeom,:,:)
2371
              Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2372
              Hess(IGeom,:,:)=HessTmp
2373
           END DO ! matches DO IGeom=IGeom0,IGeomF
2374

    
2375
           ! Second using the neighbour points
2376
           IF (HupNeighbour) THEN
2377
              ! Only one point for end points.
2378
              IF (OptReac) THEN
2379
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2380
                    GeomTmp=IntCoordF(1,:)-IntCoordF(2,:)
2381
                 ELSE
2382
                    GeomTmp=Reshape(XyzGeomF(1,:,:),(/3*Nat/))-Reshape(XyzGeomF(2,:,:),(/3*Nat/))
2383
                 END IF
2384
                 GradTmp=Grad(1,:)-Grad(2,:)
2385
                 HessTmp=Hess(1,:,:)
2386
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)
2387
                 Hess(1,:,:)=HessTmp
2388
              END IF
2389

    
2390
              IF (OptProd) THEN
2391
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2392
                    GeomTmp=IntCoordF(NGeomF,:)-IntCoordF(NGeomF-1,:)
2393
                 ELSE
2394
                    GeomTmp=Reshape(XyzGeomF(NGeomF,:,:),(/3*Nat/))-Reshape(XyzGeomF(NGeomF-1,:,:),(/3*Nat/))
2395
                 END IF
2396
                 GradTmp=Grad(NGeomF,:)-Grad(NGeomF-1,:)
2397
                 HessTmp=Hess(NGeomF,:,:)
2398
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2399
                 Hess(NGeomF,:,:)=HessTmp
2400
              END IF
2401

    
2402
              ! Two points for the rest of the path
2403
              DO IGeom=2,NGeomF-1
2404
                 IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'BAKER')) THEN
2405
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom+1,:)
2406
                 ELSE
2407
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom+1,:,:),(/3*Nat/))
2408
                 END IF
2409
                 GradTmp=Grad(IGeom,:)-Grad(IGeom+1,:)
2410
                 HessTmp=Hess(IGeom,:,:)
2411
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2412

    
2413
                 IF ((COORD.EQ.'ZMAT') .OR. (COORD.EQ.'BAKER')) THEN
2414
                    GeomTmp=IntCoordF(IGeom,:)-IntCoordF(IGeom-1,:)
2415
                 ELSE
2416
                    GeomTmp=Reshape(XyzGeomF(IGeom,:,:),(/3*Nat/))-Reshape(XyzGeomF(IGeom-1,:,:),(/3*Nat/))
2417
                 END IF
2418
                 GradTmp=Grad(IGeom,:)-Grad(IGeom-1,:)
2419

    
2420
                 Call Hupdate_all(NCoord,GeomTmp,GradTmp,HessTmp)        
2421
                 Hess(IGeom,:,:)=HessTmp
2422
              END DO ! matches DO IGeom=2,NGeomF-1
2423
           END IF !matches IF (HupNeighbour) THEN
2424
        END IF ! matches IF (FCalcHess)
2425
     END IF ! If (not.fini.).AND.(IOpt.LE.MaxCyc)
2426
  END DO ! matches DO WHILE ((IOpt.LT.MaxCyc).AND.(.NOT.Fini))
2427

    
2428
  IF (PROG=="GAUSSIAN") THEN
2429
     DEALLOCATE(Gauss_paste)
2430
     DEALLOCATE(Gauss_root,Gauss_comment,Gauss_end,current)
2431
  END IF
2432
  DEALLOCATE(XyzGeomF, IntCoordF)
2433
  DEALLOCATE(XyzGeomI, IntCoordI)
2434
  DEALLOCATE(XyzTangent,IntTangent)
2435
  DEALLOCATE(AtName,IndZmat,Order,Atome,MassAt)
2436
  DEALLOCATE(GeomOld_all)
2437

    
2438
  if (PROG=="GAUSSIAN") THEN
2439
     ! We de-allocate the Gauss_XX lists
2440
     ! transverse the list and deallocate each node
2441
     previous => Gauss_root ! make current point to head of list
2442
     DO
2443
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2444
        current => previous%next ! make list point to next node of head
2445
        DEALLOCATE(previous) ! deallocate current head node
2446
        previous => current  ! make current point to new head
2447
     END DO
2448

    
2449
     previous => Gauss_comment ! make current point to head of list
2450
     DO
2451
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2452
        current => previous%next ! make list point to next node of head
2453
        DEALLOCATE(previous) ! deallocate current head node
2454
        previous => current  ! make current point to new head
2455
     END DO
2456

    
2457
     previous => Gauss_end ! make current point to head of list
2458
     DO
2459
        IF (.NOT. ASSOCIATED(previous)) EXIT ! exit if null pointer
2460
        current => previous%next ! make list point to next node of head
2461
        DEALLOCATE(previous) ! deallocate current head node
2462
        previous => current  ! make current point to new head
2463
     END DO
2464

    
2465
     DEALLOCATE(current)
2466
  END IF
2467

    
2468
  DEALLOCATE(GeomTmp, Hess, GradTmp)
2469
  IF (COORD.EQ.'ZMAT')  DEALLOCATE(dzdc)
2470
  IF (COORD.EQ.'BAKER') THEN
2471
     DEALLOCATE(BMat_BakerT,BTransInv,BTransInv_local,Xprimitive_t)
2472
     DEALLOCATE(XprimitiveF,UMatF,BTransInvF,UMat,UMat_local)
2473
  END IF
2474

    
2475
  WRITE(IOOUT,*) "Exiting Path, IOpt=",IOpt
2476

    
2477
END PROGRAM PathOpt